
#####################################################################################################
#####################################################################################################
# This sequence/project begins with a script that uses Claude's API to obtain data about a list of 
# fragrances (perfumes, colognes, etc.). This method (using Claude) proved to be too expensive and 
# was replaced mid-project, but I include it here for full transparency. Claude's API has advanced 
# since then (January 2026), to the point where the script below might not even work if used today. 
# One of the unexpected benefits of using Claude first (for awhile) is that the new API I settled 
# on (Grok, below) showed systematic differences with Claude in its assessments of fragrances. This 
# gave rise to the creation of a correction script (shown below), which is arguably the centerpiece 
# of the manuscript for which this .r document is a supplement.
#####################################################################################################
#####################################################################################################

library(xml2)
library(rvest)
library(psych)
library(dotenv)
library(tidyverse)
library(data.table)
library(httr)
library(jsonlite)
library(h2o)
library(tidyverse)
library(dplyr)
library(readr)

# Configuration
load_dot_env("claude_api_key.env")
API_KEY <- Sys.getenv("ANTHROPIC_API_KEY")
API_URL <- "https://api.anthropic.com/v1/messages"

# Haiku 4.5: Input $0.80/MTok, Output $4/MTok
# Web search: $10/1,000 searches (same for all models)
MODEL <- "claude-haiku-4-5-20251001"

BATCH_SIZE <- 4
OUTPUT_FILE <- "fragrance_data_raw.csv"
CHECKPOINT_FILE <- "fragrances_checkpoint_claude.csv"
ERROR_LOG_FILE <- "fragrances_errors_claude.log"

# System prompt to strictly enforce CSV-only output
SYSTEM_PROMPT <- "You are a data extraction machine. You output ONLY valid CSV data.
You NEVER output explanations, commentary, markdown, apologies, or any text
that is not part of the CSV. You NEVER say things like 'I'll search for...' or
'Based on my research...' — you just output the CSV rows.
If you cannot find information about a fragrance, use 2 (Unknown) for all its values.
Every fragrance name in column 1 MUST be wrapped in double quotes.
Every row must have exactly the right number of comma-separated integer values (0-4)."

# ---------- Error logging ----------
log_error <- function(msg) {
  timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S")
  line <- sprintf("[%s] %s\n", timestamp, msg)
  cat(line)
  cat(line, file = ERROR_LOG_FILE, append = TRUE)
}

# Read the questions
questions_text <- readLines("TEMP_frag_Qs.txt", warn = FALSE)
questions <- questions_text[questions_text != ""]

# Read the fragrances
fragrances_df <- read.csv("fragrances.csv", stringsAsFactors = FALSE)
fragrances <- fragrances_df$fragrance

# Initialize results matrix
n_fragrances <- length(fragrances)
n_questions <- length(questions)
results_matrix <- matrix(2, nrow = n_fragrances, ncol = n_questions)
rownames(results_matrix) <- fragrances
colnames(results_matrix) <- paste0("Q", 1:n_questions)

# Track cost
total_input_tokens <- 0
total_output_tokens <- 0
total_searches <- 0

# Function to create the prompt for a batch of fragrances
create_prompt <- function(fragrance_batch, questions) {
  questions_formatted <- paste0(1:length(questions), ". ", questions, collapse = "\n")

  # Quote fragrance names to avoid CSV ambiguity
  fragrances_list <- paste0('- "', fragrance_batch, '"', collapse = "\n")

  # Escape any % in user data to prevent sprintf errors
  questions_safe <- gsub("%", "%%", questions_formatted)
  fragrances_safe <- gsub("%", "%%", fragrances_list)

  prompt <- sprintf('Search the web for information about these %d fragrances, then answer the questions below.

Fragrances:
%s

Questions (answer each with 0-4):
%s

Scale: 0=Definitely No, 1=Probably No, 2=Unknown, 3=Probably Yes, 4=Definitely Yes

Output ONLY a CSV table. First column "Fragrance" (values in double quotes), then Q1 through Q%d.
Exactly %d data rows. No other text. Example format:
"Fragrance",Q1,Q2,Q3
"Example Scent",4,0,2

Begin:',
    length(fragrance_batch), fragrances_safe, questions_safe,
    length(questions), length(fragrance_batch))

  return(prompt)
}

# ---------- API call ----------
call_claude_api <- function(prompt, api_key, max_retries = 6) {
  headers <- c(
    "x-api-key" = api_key,
    "anthropic-version" = "2023-06-01",
    "Content-Type" = "application/json"
  )

  body <- list(
    model = MODEL,
    max_tokens = 4096,
    system = SYSTEM_PROMPT,
    messages = list(
      list(role = "user", content = prompt)
    ),
    tools = list(
      list(
        type = "web_search_20250305",
        name = "web_search",
        max_uses = 5
      )
    ),
    temperature = 0.0  # Zero temp for maximum format compliance
  )

  for (attempt in 1:max_retries) {
    tryCatch({
      response <- POST(
        url = API_URL,
        add_headers(.headers = headers),
        body = toJSON(body, auto_unbox = TRUE),
        encode = "json",
        timeout(120)
      )

      status <- status_code(response)
      raw_body <- content(response, "text", encoding = "UTF-8")

      if (status == 200) {
        parsed <- fromJSON(raw_body, simplifyVector = FALSE)

        # Track token usage
        if (!is.null(parsed$usage)) {
          input_tok <- parsed$usage$input_tokens %||% 0
          output_tok <- parsed$usage$output_tokens %||% 0
          search_count <- 0
          if (!is.null(parsed$usage$server_tool_use)) {
            search_count <- parsed$usage$server_tool_use$web_search_requests %||% 0
          }
          total_input_tokens <<- total_input_tokens + input_tok
          total_output_tokens <<- total_output_tokens + output_tok
          total_searches <<- total_searches + search_count
          cat(sprintf("  Tokens: %d in / %d out | Searches: %d\n",
                      input_tok, output_tok, search_count))
        }

        # Handle pause_turn
        stop_reason <- parsed$stop_reason %||% ""
        all_content <- parsed$content

        while (stop_reason == "pause_turn") {
          cat("  Claude paused turn, sending continuation...\n")
          cont_body <- list(
            model = MODEL,
            max_tokens = 4096,
            system = SYSTEM_PROMPT,
            messages = list(
              list(role = "user", content = prompt),
              list(role = "assistant", content = all_content),
              list(role = "user", content = "Continue the CSV output.")
            ),
            tools = list(
              list(
                type = "web_search_20250305",
                name = "web_search",
                max_uses = 5
              )
            ),
            temperature = 0.0
          )

          cont_response <- POST(
            url = API_URL,
            add_headers(.headers = headers),
            body = toJSON(cont_body, auto_unbox = TRUE),
            encode = "json",
            timeout(120)
          )

          if (status_code(cont_response) != 200) {
            log_error(sprintf("Continuation failed with status %d", status_code(cont_response)))
            break
          }

          cont_parsed <- fromJSON(content(cont_response, "text", encoding = "UTF-8"),
                                  simplifyVector = FALSE)

          if (!is.null(cont_parsed$usage)) {
            total_input_tokens <<- total_input_tokens + (cont_parsed$usage$input_tokens %||% 0)
            total_output_tokens <<- total_output_tokens + (cont_parsed$usage$output_tokens %||% 0)
            if (!is.null(cont_parsed$usage$server_tool_use)) {
              total_searches <<- total_searches + (cont_parsed$usage$server_tool_use$web_search_requests %||% 0)
            }
          }

          all_content <- c(all_content, cont_parsed$content)
          stop_reason <- cont_parsed$stop_reason %||% ""
        }

        # Extract ALL text blocks
        text_parts <- c()
        for (block in all_content) {
          if (!is.null(block$type) && block$type == "text") {
            text_parts <- c(text_parts, block$text)
          }
        }

        result_text <- paste(text_parts, collapse = "\n")

        if (nchar(trimws(result_text)) == 0) {
          log_error("API returned 200 but no text blocks found")
          log_error(sprintf("Content block types: %s",
                            paste(sapply(all_content, function(b) b$type %||% "NULL"), collapse = ", ")))
          if (attempt < max_retries) {
            Sys.sleep(2^attempt)
            next
          }
        }

        return(result_text)

      } else if (status %in% c(429, 529)) {
        wait_time <- 2^attempt * 10
        label <- ifelse(status == 429, "Rate limited", "API overloaded")
        cat(sprintf("%s (HTTP %d). Waiting %d seconds...\n", label, status, wait_time))
        Sys.sleep(wait_time)
      } else {
        log_error(sprintf("HTTP %d: %s", status, substr(raw_body, 1, 500)))
        if (attempt < max_retries) {
          Sys.sleep(2^attempt)
        } else {
          stop(sprintf("API error after %d attempts: HTTP %d", max_retries, status))
        }
      }
    }, error = function(e) {
      log_error(sprintf("Attempt %d error: %s", attempt, e$message))
      if (attempt == max_retries) {
        stop(sprintf("Failed after %d attempts: %s", max_retries, e$message))
      }
      Sys.sleep(2^attempt)
    })
  }
}

# ---------- Parse response ----------
parse_response <- function(response_text, fragrance_batch) {
  if (is.null(response_text) || nchar(trimws(response_text)) == 0) {
    log_error("Empty response text received")
    return(matrix(2, nrow = length(fragrance_batch), ncol = n_questions))
  }

  tryCatch({
    # Strip markdown code fences
    cleaned <- gsub("```csv\\s*\n?", "", response_text)
    cleaned <- gsub("```\\s*$", "", cleaned)
    cleaned <- gsub("```", "", cleaned)
    cleaned <- trimws(cleaned)

    lines <- strsplit(cleaned, "\n")[[1]]
    lines <- trimws(lines)
    lines <- lines[lines != ""]  # Remove blank lines

    # Find the header line (contains "Fragrance" or "Q1")
    header_idx <- grep('^"?Fragrance"?,', lines, ignore.case = TRUE)
    if (length(header_idx) == 0) {
      header_idx <- grep("Q1", lines)
    }

    if (length(header_idx) > 0) {
      start <- header_idx[1]
      end <- min(start + length(fragrance_batch), length(lines))
      csv_lines <- lines[start:end]
      cleaned <- paste(csv_lines, collapse = "\n")
      cat(sprintf("  Extracted CSV from lines %d-%d of %d total lines\n",
                  start, end, length(lines)))
    } else {
      # Last resort: try to find lines that look like CSV data rows
      # (start with a quote and contain mostly commas and digits)
      data_lines <- grep('^".*"(,[0-4]){5,}', lines)
      if (length(data_lines) > 0) {
        cat(sprintf("  Found %d data-like lines by pattern matching\n", length(data_lines)))
        # Construct a header
        n_cols <- length(strsplit(lines[data_lines[1]], ",")[[1]]) - 1
        header <- paste0('"Fragrance",', paste0("Q", 1:n_cols, collapse = ","))
        csv_lines <- c(header, lines[data_lines])
        cleaned <- paste(csv_lines, collapse = "\n")
      } else {
        log_error(sprintf("No CSV header or data rows found. First 300 chars: %s",
                          substr(cleaned, 1, 300)))
      }
    }

    # Parse with read.csv — the quoting should handle commas in fragrance names
    response_df <- read.csv(text = cleaned, header = TRUE, stringsAsFactors = FALSE,
                            quote = '"', check.names = FALSE)

    # Extract numeric columns (skip the Fragrance name column)
    # Find which column is the fragrance name (usually first)
    name_col <- which(grepl("fragrance", names(response_df), ignore.case = TRUE))[1]
    if (is.na(name_col)) name_col <- 1

    numeric_cols <- response_df[, -name_col, drop = FALSE]

    # Force all columns to numeric
    result_matrix <- matrix(
      as.numeric(as.matrix(numeric_cols)),
      nrow = nrow(numeric_cols),
      ncol = ncol(numeric_cols)
    )

    # Validate values are in range 0-4
    result_matrix[is.na(result_matrix)] <- 2
    result_matrix[result_matrix < 0 | result_matrix > 4] <- 2

    # Validate row count
    if (nrow(result_matrix) != length(fragrance_batch)) {
      log_error(sprintf("Row mismatch: got %d, expected %d", nrow(result_matrix), length(fragrance_batch)))
      if (nrow(result_matrix) < length(fragrance_batch)) {
        padding <- matrix(2, nrow = length(fragrance_batch) - nrow(result_matrix),
                          ncol = ncol(result_matrix))
        result_matrix <- rbind(result_matrix, padding)
      } else {
        result_matrix <- result_matrix[1:length(fragrance_batch), , drop = FALSE]
      }
    }

    # Validate column count
    if (ncol(result_matrix) != n_questions) {
      log_error(sprintf("Column mismatch: got %d, expected %d", ncol(result_matrix), n_questions))
      if (ncol(result_matrix) > n_questions) {
        result_matrix <- result_matrix[, 1:n_questions, drop = FALSE]
      } else {
        # Pad with 2s
        padding <- matrix(2, nrow = nrow(result_matrix),
                          ncol = n_questions - ncol(result_matrix))
        result_matrix <- cbind(result_matrix, padding)
      }
    }

    cat(sprintf("  Successfully parsed %dx%d matrix\n", nrow(result_matrix), ncol(result_matrix)))
    return(result_matrix)

  }, error = function(e) {
    log_error(sprintf("Parse error: %s", e$message))
    log_error(sprintf("Raw response (first 500 chars): %s", substr(response_text, 1, 500)))
    return(matrix(2, nrow = length(fragrance_batch), ncol = n_questions))
  })
}

# Function to save checkpoint
save_checkpoint <- function(results_matrix, last_processed_index) {
  checkpoint_data <- as.data.frame(results_matrix[1:last_processed_index, , drop = FALSE])
  # Quote fragrance names to handle special characters
  checkpoint_data <- cbind(Fragrance = rownames(results_matrix)[1:last_processed_index],
                           checkpoint_data)
  write.csv(checkpoint_data, CHECKPOINT_FILE, row.names = FALSE)
  cat(sprintf("Checkpoint saved at fragrance %d\n", last_processed_index))
}

# Check for existing checkpoint
start_index <- 1
if (file.exists(CHECKPOINT_FILE)) {
  checkpoint_data <- read.csv(CHECKPOINT_FILE, stringsAsFactors = FALSE)
  if (nrow(checkpoint_data) > 0) {
    cat("Found checkpoint file. Loading previous progress...\n")
    start_index <- nrow(checkpoint_data) + 1

    for (i in 1:nrow(checkpoint_data)) {
      fragrance_name <- checkpoint_data$Fragrance[i]
      fragrance_idx <- which(fragrances == fragrance_name)
      if (length(fragrance_idx) > 0) {
        results_matrix[fragrance_idx, ] <- as.numeric(checkpoint_data[i, -1])
      }
    }
    cat(sprintf("Resuming from fragrance %d\n", start_index))
  } else {
    cat("Checkpoint file exists but is empty, starting fresh.\n")
  }
}

# ========== DRY RUN: Test a single batch first ==========
cat("\n===== DRY RUN: Testing with first batch =====\n")
test_batch <- fragrances[1:min(BATCH_SIZE, n_fragrances)]
test_prompt <- create_prompt(test_batch, questions)
cat(sprintf("Prompt length: %d characters\n", nchar(test_prompt)))

test_response <- call_claude_api(test_prompt, API_KEY)
cat(sprintf("\nRaw response (first 1000 chars):\n%s\n", substr(test_response, 1, 1000)))
cat("---\n")

test_parsed <- parse_response(test_response, test_batch)
cat(sprintf("Parsed matrix: %d x %d\n", nrow(test_parsed), ncol(test_parsed)))
cat("Values:\n")
print(test_parsed)

# Check if it's all 2s (failure)
if (all(test_parsed == 2)) {
  cat("\n*** WARNING: Dry run returned all 2s -- parsing likely failed! ***\n")
  cat("*** Check the raw response above and the error log. ***\n")
  cat("*** Fix the issue before running the full batch. ***\n")
  est <- total_input_tokens / 1e6 * 0.80 + total_output_tokens / 1e6 * 4.00 + total_searches * 0.01
  cat(sprintf("\nEstimated cost of dry run: $%.4f\n", est))
  stop("Dry run failed. Aborting to prevent wasting money.")
}

cat("\n===== Dry run passed! Proceeding with full run =====\n\n")
# ==========================================================

# Main processing loop
cat(sprintf("Processing %d fragrances in batches of %d...\n", n_fragrances, BATCH_SIZE))
for (i in seq(start_index, n_fragrances, by = BATCH_SIZE)) {
  end_idx <- min(i + BATCH_SIZE - 1, n_fragrances)
  fragrance_batch <- fragrances[i:end_idx]

  cat(sprintf("\nProcessing fragrances %d-%d of %d...\n", i, end_idx, n_fragrances))

  prompt <- create_prompt(fragrance_batch, questions)

  max_retries <- 4
  attempt <- 1
  batch_results <- parse_response(call_claude_api(prompt, API_KEY), fragrance_batch)
  while (attempt < max_retries &&
         (ncol(batch_results) != n_questions ||
          nrow(batch_results) != length(fragrance_batch) ||
          all(batch_results == 2))) {
    cat(sprintf("Incomplete/failed response (attempt %d/%d), retrying...\n", attempt, max_retries))
    Sys.sleep(5 * attempt)
    batch_results <- parse_response(call_claude_api(prompt, API_KEY), fragrance_batch)
    attempt <- attempt + 1
  }
  if (ncol(batch_results) != n_questions || nrow(batch_results) != length(fragrance_batch)) {
    log_error(sprintf("Batch %d-%d failed after retries, filling with unknowns", i, end_idx))
    batch_results <- matrix(2, nrow = length(fragrance_batch), ncol = n_questions)
  }

  results_matrix[i:end_idx, ] <- batch_results

  # Save checkpoint every 15 fragrances
  if ((end_idx %% 15) == 0 || end_idx == n_fragrances) {
    save_checkpoint(results_matrix, end_idx)
  }

  if (end_idx < n_fragrances) {
    Sys.sleep(2)
  }

  progress_pct <- round(100 * end_idx / n_fragrances, 1)
  est_cost <- total_input_tokens / 1e6 * 0.80 +
              total_output_tokens / 1e6 * 4.00 +
              total_searches * 0.01
  cat(sprintf("Progress: %.1f%% | Est. cost so far: $%.2f | Searches: %d\n",
              progress_pct, est_cost, total_searches))
}

# Save final results
cat("\nSaving final results...\n")
results_df <- as.data.frame(results_matrix)
results_df <- cbind(Fragrance = rownames(results_matrix), results_df)
write.csv(results_df, OUTPUT_FILE, row.names = FALSE)
cat(sprintf("\nAnalysis complete! Results saved to %s\n", OUTPUT_FILE))
cat(sprintf("Matrix dimensions: %d fragrances x %d questions\n", n_fragrances, n_questions))

if (file.exists(CHECKPOINT_FILE)) {
  file.remove(CHECKPOINT_FILE)
  cat("Checkpoint file removed.\n")
}

# Final cost summary
cat(sprintf("\n===== Cost Summary =====\n"))
cat(sprintf("Total input tokens:  %d\n", total_input_tokens))
cat(sprintf("Total output tokens: %d\n", total_output_tokens))
cat(sprintf("Total web searches:  %d\n", total_searches))
est_cost <- total_input_tokens / 1e6 * 0.80 +
            total_output_tokens / 1e6 * 4.00 +
            total_searches * 0.01
cat(sprintf("Estimated total cost: $%.2f\n", est_cost))

# Summary statistics
cat("\nResponse distribution:\n")
response_counts <- table(as.vector(results_matrix))
for (val in 0:4) {
  if (val %in% names(response_counts)) {
    pct <- round(100 * response_counts[as.character(val)] / (n_fragrances * n_questions), 2)
    cat(sprintf(" %d: %d (%.2f%%)\n", val, response_counts[as.character(val)], pct))
  }
}

#####################################################################################################
#####################################################################################################
# This concludes the use of Claude's API. It was used to collect data on approximately ~15,000 
# fragrances and 60 items. I then switched to the Grok API, which is used in the script below.
#####################################################################################################
#####################################################################################################

library(psych)
library(httr)
library(jsonlite)
library(dplyr)
library(readr)
library(dotenv)
library(tidyverse)
library(h2o)

# This script processes fragrance names (perfumes, colognes, etc.) and answers
# questions about each one, using the xAI BATCH API (asynchronous) rather than
# one synchronous call per batch. It mirrors the FashionDive batch workflow:
#   build JSONL -> upload -> create batch -> poll -> paginate results ->
#   robust parse (name match + column alignment) -> checkpoint -> retry JSONL.

# ====================== CONFIGURATION ======================
load_dot_env("grok_api_key.env")
API_KEY <- Sys.getenv("XAI_API_KEY")
BASE_URL <- "https://api.x.ai/v1"

# Model + reasoning effort (verified against docs.x.ai, July 2026).
#   * grok-4-1-fast-reasoning was RETIRED 2026-05-15. Requests to that slug now
#     silently redirect to grok-4.3 at LOW reasoning effort and are billed at
#     grok-4.3 rates ($1.25/1M input, $2.50/1M output) -- so keeping the old slug
#     saves nothing and hides the effort setting. We call grok-4.3 explicitly.
#   * grok-4.5 is NOT supported by the Batch API and would be rejected.
MODEL <- "grok-4.3"

# Reasoning effort: "low" / "medium" / "high" (set to "" to omit and take the
# model default). These characterization questions lean on judgment + web search,
# so we default to "high" for accuracy. NOTE: on a run of this size, high effort
# materially increases reasoning-token cost -- drop to "medium" or "low" if you
# want to trade some accuracy for cost/speed.
REASONING_EFFORT <- "medium"

BATCH_SIZE <- 5 # Number of fragrances to process per API call
OUTPUT_FILE <- "fragrance_data_raw_forPub.csv"
CHECKPOINT_FILE <- "checkpoint_1.csv"
BATCH_JSONL_FILE <- "batch_engine.jsonl"
BATCH_ID_FILE <- "batch_id.txt"

# To force a brand new batch (ignoring any saved batch ID), set this to TRUE.
# Otherwise, the script auto-resumes if BATCH_ID_FILE exists with a valid batch.
FORCE_NEW_BATCH <- FALSE

# Set to TRUE when you have a pre-built JSONL file (e.g., a retry file) that you
# want to upload as-is. When FALSE (default), the script generates a fresh JSONL
# from the full fragrance list.
SKIP_JSONL_GENERATION <- FALSE

# ====================== LOAD DATA ======================
questions <- readLines("TEMP_frag_Qs.txt", warn = FALSE)
questions <- questions[questions != ""]

fragrances_df <- read.csv("new_fragrances.csv", stringsAsFactors = FALSE, fileEncoding = "UTF-8")
fragrances <- fragrances_df$fragrance

# Clean only problematic characters, preserve valid UTF-8 (Cyrillic, accents, etc.)
sanitize_for_json <- function(x) {
  x <- gsub('"', "'", x, fixed = TRUE)
  x <- gsub('`', "'", x, fixed = TRUE)
  x <- iconv(x, from = "UTF-8", to = "UTF-8", sub = "")
  x <- gsub("\\\\", "", x)
  x
}
fragrances <- sapply(fragrances, sanitize_for_json, USE.NAMES = FALSE)

n_fragrances <- length(fragrances)
n_questions <- length(questions)
results_matrix <- matrix(2, nrow = n_fragrances, ncol = n_questions) # 2 = "Unknown"
rownames(results_matrix) <- fragrances
colnames(results_matrix) <- paste0("Q", 1:n_questions)

# ====================== CREATE PROMPT ======================
create_prompt <- function(fragrance_batch, questions) {
  questions_formatted <- paste0(1:length(questions), ". ", questions, collapse = "\n")
  fragrances_list <- paste0("- ", fragrance_batch, collapse = "\n")

  sprintf("I need you to answer questions about the following %d fragrances (colognes, perfumes, etc.):
%s
For each fragrance, please answer these %d questions using this scale:
0 = Definitely No
1 = Probably No
2 = Unknown (or not applicable)
3 = Probably Yes
4 = Definitely Yes
Questions:
%s
Please respond with ONLY a CSV format where each row represents a fragrance (in the same order as listed above) and each column represents the answer to the corresponding question (Q1 through Q%d). Include a header row with the fragrance names in the first column and Q1, Q2, etc. for the question columns.
Important:
- Many fragrances are quite obscure and the questions are detailed, so you will almost certainly need to search the internet thoroughly for accurate and complete information about each fragrance before attempting to answer the questions
- Provide exactly %d rows of data (one per fragrance)
- Each row must have EXACTLY %d numeric values (one per question)
- NEVER leave any cell empty; if you can't find the answer to a question even after effortful and meticulous search, use 2
- No additional text or explanation, just the CSV data",
    length(fragrance_batch), fragrances_list, n_questions, questions_formatted,
    n_questions, length(fragrance_batch), n_questions)
}

# Build one JSONL request object for a batch. Adds the reasoning-effort field
# only when REASONING_EFFORT is a non-empty string, so setting it to "" cleanly
# omits the field and falls back to the model default.
make_request <- function(custom_id, prompt) {
  body <- list(
    model = MODEL,
    input = list(list(role = "user", content = prompt)),
    tools = list(list(type = "web_search"))
  )
  if (!is.null(REASONING_EFFORT) && !is.na(REASONING_EFFORT) && nzchar(REASONING_EFFORT)) {
    body$reasoning <- list(effort = REASONING_EFFORT)   # -> {"reasoning":{"effort":"..."}}
  }
  list(custom_id = custom_id, method = "POST", url = "/v1/responses", body = body)
}


# ====================== PARSE & CHECKPOINT ======================

# Normalize a fragrance name for fuzzy matching: lowercase, strip any
# "(...)" parenthetical, drop punctuation and whitespace. This lets
# "Aventus (Creed)" match a returned "Aventus", etc.
normalize_fragrance_name <- function(x) {
  x <- as.character(x)
  x <- tolower(x)
  x <- gsub("\\(.*?\\)", "", x)        # remove parenthetical
  x <- gsub("https?://", "", x)        # strip protocol if present
  x <- gsub("www\\.", "", x)           # strip www.
  x <- gsub("[^a-z0-9]", "", x)        # keep only alphanumerics
  trimws(x)
}

# Match a single returned name against the normalized batch names.
# Returns the index of the unique best match, or NA if none / ambiguous.
match_fragrance_name <- function(returned_norm, batch_norm) {
  if (is.na(returned_norm) || nchar(returned_norm) == 0) return(NA_integer_)

  # 1. Exact normalized match
  exact <- which(batch_norm == returned_norm)
  if (length(exact) == 1) return(exact)
  if (length(exact) > 1) return(NA_integer_)  # ambiguous

  # 2. Substring match (returned name contained in a batch name or vice versa).
  #    Require the shorter string be at least 4 chars to avoid spurious hits.
  contains <- which(sapply(batch_norm, function(b) {
    if (nchar(b) == 0 || nchar(returned_norm) == 0) return(FALSE)
    shorter <- min(nchar(b), nchar(returned_norm))
    shorter >= 4 && (grepl(returned_norm, b, fixed = TRUE) || grepl(b, returned_norm, fixed = TRUE))
  }))
  if (length(contains) == 1) return(contains)
  if (length(contains) > 1) return(NA_integer_)  # ambiguous

  NA_integer_
}

# Returns a list: $mat (the matrix), $status, $detail.
# Statuses:
#   "ok"            - clean parse, every batch fragrance matched by name
#   "col_aligned"   - Q columns realigned by header; rows still matched by name
#   "row_matched"   - rows placed by name match, but some fragrances had no returned row (left as 2s)
#   "row_unmatched" - one or more returned rows could not be matched to a batch fragrance
#   "col_shift"     - Q columns could not be aligned by header (data discarded upstream)
#   "error"         - parse failed entirely
parse_response <- function(response_text, fragrance_batch) {
  status <- "ok"
  detail <- ""

  tryCatch({
    response_text <- gsub("```csv\\s*\n|```", "", response_text)
    response_text <- trimws(response_text)
    df <- read.csv(text = response_text, header = TRUE, stringsAsFactors = FALSE, check.names = FALSE)

    col_names <- colnames(df)
    data_cols <- col_names[-1]
    q_numbers <- suppressWarnings(as.integer(gsub("[^0-9]", "", data_cols)))

    # Returned fragrance names live in the first column. Keep them for matching.
    returned_names <- as.character(df[[1]])

    # Numeric answer matrix (everything except the name column).
    numeric_df <- df[, -1, drop = FALSE]
    numeric_df[] <- lapply(numeric_df, function(col) as.numeric(trimws(as.character(col))))
    raw_mat <- as.matrix(numeric_df)

    n_batch <- length(fragrance_batch)

    # ---------- Step 1: build a per-returned-row answer matrix aligned to n_questions ----------
    # (Column alignment first, so each returned row has a clean Q1..Qn vector.)
    aligned <- matrix(2, nrow = nrow(raw_mat), ncol = n_questions)

    if (all(!is.na(q_numbers)) && all(q_numbers >= 1) && all(q_numbers <= n_questions)) {
      for (ci in seq_along(q_numbers)) {
        aligned[, q_numbers[ci]] <- raw_mat[, ci]
      }
      if (ncol(raw_mat) != n_questions) {
        missing_qs <- setdiff(1:n_questions, q_numbers)
        status <- "col_aligned"
        detail <- paste0(detail, sprintf("Aligned by header. Missing: Q%s. ", paste(missing_qs, collapse = ",")))
      }
    } else if (ncol(raw_mat) == n_questions) {
      aligned <- raw_mat
    } else {
      # Columns can't be trusted — signal col_shift; caller discards the batch.
      detail <- paste0(detail, sprintf("Got %d cols, expected %d. Headers: %s. COLUMN ALIGNMENT FAILED. ",
                                        ncol(raw_mat), n_questions, paste(head(data_cols, 5), collapse = ",")))
      return(list(mat = matrix(2, nrow = n_batch, ncol = n_questions),
                  status = "col_shift", detail = detail))
    }
    aligned[is.na(aligned)] <- 2

    # ---------- Step 2: place each returned row into the correct fragrance slot by NAME ----------
    mat <- matrix(2, nrow = n_batch, ncol = n_questions)
    batch_norm    <- normalize_fragrance_name(fragrance_batch)
    returned_norm <- normalize_fragrance_name(returned_names)

    filled <- rep(FALSE, n_batch)   # which batch slots got real data
    unmatched_rows <- integer(0)    # returned rows we couldn't place

    for (r in seq_len(nrow(aligned))) {
      idx <- match_fragrance_name(returned_norm[r], batch_norm)
      if (!is.na(idx) && !filled[idx]) {
        mat[idx, ] <- aligned[r, ]
        filled[idx] <- TRUE
      } else {
        unmatched_rows <- c(unmatched_rows, r)
      }
    }

    n_matched   <- sum(filled)
    n_missing   <- sum(!filled)
    n_unmatched <- length(unmatched_rows)

    # ---------- Step 3: positional fallback ONLY when it is unambiguously safe ----------
    # If nothing matched by name AND the row count is exactly right, fall back to
    # positional assignment (covers the case where Grok returned numeric IDs or
    # otherwise unrecognizable names but in the correct order).
    if (n_matched == 0 && nrow(aligned) == n_batch) {
      mat <- aligned
      status <- if (status == "ok") "row_fix" else status
      detail <- paste0(detail, "No names matched; used positional fallback (row count exact). ")
    } else {
      # Escalate status based on coverage problems.
      if (n_unmatched > 0) {
        status <- "row_unmatched"
        detail <- paste0(detail, sprintf("%d returned row(s) unmatched to any batch fragrance. ", n_unmatched))
      }
      if (n_missing > 0) {
        # Don't downgrade a more severe status.
        if (!status %in% c("row_unmatched", "col_aligned", "col_shift", "error")) {
          status <- "row_matched"
        }
        missing_fragrances <- fragrance_batch[!filled]
        detail <- paste0(detail, sprintf("%d fragrance(s) had no matching row (left as 2s): %s. ",
                                          n_missing, paste(head(missing_fragrances, 5), collapse = "; ")))
      }
    }

    mat[is.na(mat)] <- 2

    list(mat = mat, status = status, detail = detail)
  }, error = function(e) {
    list(mat = matrix(2, nrow = length(fragrance_batch), ncol = n_questions),
         status = "error",
         detail = e$message)
  })
}

save_checkpoint <- function() {
  df <- as.data.frame(results_matrix)
  df <- cbind(Fragrance = rownames(results_matrix), df)
  write.csv(df, CHECKPOINT_FILE, row.names = FALSE)
}

# Load checkpoint
if (file.exists(CHECKPOINT_FILE)) {
  cat("Loading checkpoint...\n")
  cp <- read.csv(CHECKPOINT_FILE, stringsAsFactors = FALSE)
  for (i in seq_len(nrow(cp))) {
    idx <- which(fragrances == cp$Fragrance[i])
    if (length(idx) > 0) results_matrix[idx, ] <- as.numeric(cp[i, -1])
  }
  cat(sprintf("Loaded checkpoint with %d fragrances.\n", nrow(cp)))
}

# ====================== BATCH CREATION OR RESUME ======================
cat("=== Starting xAI Batch API run (FragranceDive) ===\n")

resume_mode <- FALSE
if (!FORCE_NEW_BATCH && file.exists(BATCH_ID_FILE)) {
  saved_batch_id <- trimws(readLines(BATCH_ID_FILE, warn = FALSE)[1])
  if (nchar(saved_batch_id) > 0) {
    cat(sprintf("Found existing batch ID: %s\n", saved_batch_id))
    cat("Checking batch status...\n")

    check_resp <- GET(
      paste0(BASE_URL, "/batches/", saved_batch_id),
      add_headers(Authorization = paste("Bearer", API_KEY))
    )

    if (status_code(check_resp) == 200) {
      check_st <- content(check_resp, "parsed")
      cancel_time <- check_st$cancel_time
      is_cancelled <- is.character(cancel_time) && nchar(cancel_time) > 0

      if (is_cancelled) {
        cancel_msg <- check_st$cancel_by_xai_message
        cat(sprintf("Batch was cancelled: %s\n", if (is.character(cancel_msg)) cancel_msg else "unknown reason"))

        n_completed <- (check_st$state$num_success %||% 0) + (check_st$state$num_error %||% 0)
        if (n_completed > 0) {
          cat(sprintf("Batch has %d completed results. Will retrieve them.\n", n_completed))
          batch_id <- saved_batch_id
          resume_mode <- TRUE
        } else {
          cat("No completed results. Will create a new batch.\n")
        }
      } else {
        batch_id <- saved_batch_id
        resume_mode <- TRUE
        cat(sprintf("Resuming batch. Current state: %d/%d complete (%d pending)\n",
                    (check_st$state$num_success %||% 0) + (check_st$state$num_error %||% 0),
                    check_st$state$num_requests %||% 0,
                    check_st$state$num_pending %||% 0))
      }
    } else {
      cat("Could not retrieve saved batch. Will create a new one.\n")
    }
  }
}

if (!resume_mode) {
  # ---------- Generate JSONL (or use existing file) ----------
  if (!SKIP_JSONL_GENERATION) {
    cat("Generating batch requests...\n")
    con <- file(BATCH_JSONL_FILE, "w")
    for (i in seq(1, n_fragrances, by = BATCH_SIZE)) {
      end_i <- min(i + BATCH_SIZE - 1, n_fragrances)
      batch_fragrances <- fragrances[i:end_i]
      prompt <- create_prompt(batch_fragrances, questions)
      custom_id <- paste0("req_", i, "_", end_i)
      req <- make_request(custom_id, prompt)
      writeLines(toJSON(req, auto_unbox = TRUE), con)
    }
    close(con)
    n_requests <- length(readLines(BATCH_JSONL_FILE))
    cat(sprintf("Created %s with %d requests\n", BATCH_JSONL_FILE, n_requests))
  } else {
    cat(sprintf("Using existing JSONL file: %s\n", BATCH_JSONL_FILE))
    if (!file.exists(BATCH_JSONL_FILE)) stop("BATCH_JSONL_FILE not found: ", BATCH_JSONL_FILE)
    n_requests <- length(readLines(BATCH_JSONL_FILE))
    cat(sprintf("  File contains %d requests\n", n_requests))
  }

  # ---------- Validate JSONL ----------
  cat("Validating JSONL...\n")
  jsonl_lines <- readLines(BATCH_JSONL_FILE, warn = FALSE)
  jsonl_lines <- jsonl_lines[jsonl_lines != ""]
  for (li in seq_along(jsonl_lines)) {
    tryCatch({
      parsed <- fromJSON(jsonl_lines[li], simplifyVector = FALSE)
      required <- c("custom_id", "method", "url", "body")
      missing_fields <- setdiff(required, names(parsed))
      if (length(missing_fields) > 0) {
        cat(sprintf("  WARNING line %d: missing fields: %s\n", li, paste(missing_fields, collapse = ", ")))
      }
      if (li <= 2) {
        cat(sprintf("  Line %d: custom_id=%s, url=%s, model=%s\n",
                    li, parsed$custom_id, parsed$url, parsed$body$model))
      }
    }, error = function(e) {
      cat(sprintf("  ERROR line %d: %s\n", li, e$message))
      stop("JSONL validation failed.")
    })
  }
  cat("JSONL validation passed.\n")

  # ---------- Upload JSONL ----------
  cat("Uploading JSONL file...\n")
  upload_resp <- POST(
    paste0(BASE_URL, "/files"),
    add_headers(Authorization = paste("Bearer", API_KEY)),
    body = list(file = upload_file(BATCH_JSONL_FILE)),
    encode = "multipart"
  )
  if (status_code(upload_resp) != 200) {
    cat("File upload failed:\n")
    print(content(upload_resp, "text"))
    stop("File upload failed.")
  }
  file_id <- content(upload_resp, "parsed")$id
  cat(sprintf("File uploaded. ID: %s\n", file_id))

  # ---------- Create batch ----------
  cat("Creating batch...\n")
  batch_body <- list(
    name = paste0("FragranceDive_", format(Sys.time(), "%Y%m%d_%H%M")),
    input_file_id = file_id
  )
  resp <- POST(
    paste0(BASE_URL, "/batches"),
    add_headers(Authorization = paste("Bearer", API_KEY), `Content-Type` = "application/json"),
    body = toJSON(batch_body, auto_unbox = TRUE)
  )
  if (status_code(resp) != 200) {
    cat("Batch creation failed:\n")
    print(content(resp, "text"))
    stop("Batch creation failed.")
  }
  batch_id <- content(resp, "parsed")$batch_id
  writeLines(batch_id, BATCH_ID_FILE)
  cat(sprintf("Batch created! ID: %s\n", batch_id))
}

# ====================== POLL FOR COMPLETION ======================
cat("\nPolling for completion...\n")
cat("  Safe to Ctrl+C and re-run later — the script will resume from the saved batch ID.\n\n")

zero_request_polls <- 0
MAX_ZERO_POLLS <- 20
poll_interval <- 30
batch_done <- FALSE

for (attempt in 1:2400) {
  status_resp <- GET(
    paste0(BASE_URL, "/batches/", batch_id),
    add_headers(Authorization = paste("Bearer", API_KEY))
  )
  st <- content(status_resp, "parsed")

  if (attempt == 1) {
    cat("--- First poll ---\n")
    cat(toJSON(st, auto_unbox = TRUE, pretty = TRUE), "\n---\n\n")
  }

  # Check cancellation
  cancel_time <- st$cancel_time
  if (is.character(cancel_time) && nchar(cancel_time) > 0) {
    cancel_msg <- st$cancel_by_xai_message
    cat(sprintf("\n*** BATCH CANCELLED at %s ***\n", cancel_time))
    if (is.character(cancel_msg)) cat(sprintf("Reason: %s\n", cancel_msg))
    n_completed <- (st$state$num_success %||% 0)
    if (n_completed > 0) {
      cat(sprintf("Will retrieve %d completed results.\n", n_completed))
    } else {
      cat("No completed results available.\n")
      stop("Batch cancelled with no results.")
    }
    batch_done <- TRUE
    break
  }

  num_pending  <- st$state$num_pending  %||% 0
  num_success  <- st$state$num_success  %||% 0
  num_error    <- st$state$num_error    %||% 0
  num_requests <- st$state$num_requests %||% 0

  elapsed_min <- round((attempt - 1) * poll_interval / 60, 1)
  cat(sprintf("   [%s] (~%.0f min) %d/%d complete (%d success, %d error, %d pending)\n",
              format(Sys.time(), "%H:%M:%S"), elapsed_min,
              num_success + num_error, num_requests,
              num_success, num_error, num_pending))

  if (num_requests > 0 && num_pending == 0) {
    cat("Batch processing complete!\n")
    batch_done <- TRUE
    break
  }

  if (num_requests == 0) {
    zero_request_polls <- zero_request_polls + 1
    if (zero_request_polls == 5) {
      cat("  NOTE: num_requests still 0 after ~2.5 min.\n")
      cat(toJSON(st, auto_unbox = TRUE, pretty = TRUE), "\n")
    }
    if (zero_request_polls >= MAX_ZERO_POLLS) {
      cat(toJSON(st, auto_unbox = TRUE, pretty = TRUE), "\n")
      stop("Batch stalled with 0 requests. JSONL likely failed validation.")
    }
  } else {
    zero_request_polls <- 0
    if (poll_interval < 60) {
      poll_interval <- 60
      cat("  Requests ingested. Switching to 60s poll interval.\n")
    }
  }

  Sys.sleep(poll_interval)
}

if (!batch_done) {
  cat("Polling loop exhausted. Check batch status in xAI Console and re-run to resume.\n")
  stop("Polling timeout.")
}

# ====================== RETRIEVE RESULTS ======================
cat("\nRetrieving results...\n")
all_results <- list()
pagination_token <- NULL
page_num <- 0

repeat {
  results_url <- paste0(BASE_URL, "/batches/", batch_id, "/results?limit=100")
  if (!is.null(pagination_token)) {
    results_url <- paste0(results_url, "&pagination_token=", URLencode(pagination_token, reserved = TRUE))
  }

  results_resp <- GET(results_url, add_headers(Authorization = paste("Bearer", API_KEY)))

  if (status_code(results_resp) != 200) {
    cat(sprintf("Results retrieval failed on page %d:\n", page_num + 1))
    print(content(results_resp, "text"))
    if (length(all_results) > 0) {
      cat("Proceeding with partial results.\n")
      break
    }
    stop("Could not retrieve any results.")
  }

  page <- content(results_resp, "parsed")
  page_num <- page_num + 1

  n_page_results <- length(page$results %||% list())
  if (n_page_results > 0) {
    all_results <- c(all_results, page$results)
  }

  if (page_num %% 10 == 0 || n_page_results > 0) {
    cat(sprintf("  Page %d: +%d results (total: %d)\n", page_num, n_page_results, length(all_results)))
  }

  pagination_token <- page$pagination_token
  if (is.null(pagination_token)) break
}

cat(sprintf("Retrieved %d results across %d pages.\n\n", length(all_results), page_num))

# ====================== PARSE RESULTS ======================
cat("Parsing results...\n")
n_success <- 0
n_fail <- 0
n_web_search_used <- 0
n_web_search_missing <- 0
checkpoint_every <- 100

# Brief structure check on first result
if (length(all_results) > 0) {
  first_br <- all_results[[1]]$batch_result$response
  if (!is.null(first_br$chat_get_completion)) {
    n_ch <- length(first_br$chat_get_completion$choices)
    last_fr <- first_br$chat_get_completion$choices[[n_ch]]$finish_reason %||% "?"
    last_len <- nchar(first_br$chat_get_completion$choices[[n_ch]]$message$content %||% "")
    cat(sprintf("  First result: %d choices, last finish_reason='%s', content length=%d\n\n", n_ch, last_fr, last_len))
  }
}

has_web_search <- function(x) {
  if (is.character(x)) return(grepl("web_search", x, fixed = TRUE))
  if (is.list(x)) return(any(sapply(x, has_web_search)))
  FALSE
}

# Track problematic requests for potential retry
problem_log <- data.frame(
  custom_id = character(),
  start_i = integer(),
  end_i = integer(),
  status = character(),
  detail = character(),
  stringsAsFactors = FALSE
)

for (ri in seq_along(all_results)) {
  result <- all_results[[ri]]

  custom_id <- result$batch_request_id %||% result$custom_id %||% NULL

  if (is.null(custom_id)) { n_fail <- n_fail + 1; next }

  id_parts <- strsplit(gsub("req_", "", custom_id), "_")[[1]]
  if (length(id_parts) < 2) { n_fail <- n_fail + 1; next }
  start_i <- as.integer(id_parts[1])
  end_i   <- as.integer(id_parts[2])
  batch_fragrances <- fragrances[start_i:end_i]

  # Web search check
  if (has_web_search(result)) {
    n_web_search_used <- n_web_search_used + 1
  } else {
    n_web_search_missing <- n_web_search_missing + 1
  }

  # Extract response text
  resp_text <- NULL
  tryCatch({
    br <- result$batch_result$response

    # Strategy 1: chat_get_completion — find choice with finish_reason="stop"
    if (!is.null(br$chat_get_completion)) {
      choices <- br$chat_get_completion$choices
      if (is.list(choices)) {
        for (ch in choices) {
          fr <- ch$finish_reason %||% ""
          ct <- ch$message$content %||% ""
          if (fr == "stop" && is.character(ct) && nchar(ct) > 0) {
            resp_text <- ct
          }
        }
        if (is.null(resp_text) && length(choices) > 0) {
          ct <- choices[[length(choices)]]$message$content %||% ""
          if (is.character(ct) && nchar(ct) > 0) resp_text <- ct
        }
      }
    }

    # Strategy 2: responses API — output_text blocks
    if (is.null(resp_text)) {
      output_list <- br$responses$output %||% br$output %||% NULL
      if (is.list(output_list)) {
        for (item in output_list) {
          if (identical(item$type, "output_text") && !is.null(item$text)) {
            resp_text <- item$text
            break
          }
        }
      }
    }

    if (is.null(resp_text) && n_fail < 5) {
      raw_json <- toJSON(br, auto_unbox = TRUE)
      cat(sprintf("  Debug [%s]: No text. Keys: %s. First 300 chars: %s\n",
                  custom_id, paste(names(br), collapse=","), substr(as.character(raw_json), 1, 300)))
    }
  }, error = function(e) {
    cat(sprintf("  Error [%s]: %s\n", custom_id, e$message))
  })

  if (!is.null(resp_text) && nchar(resp_text) > 0) {
    tryCatch({
      parsed <- parse_response(resp_text, batch_fragrances)

      # col_shift now returns all-2s directly from parse_response (columns couldn't
      # be aligned by header). This block is a defensive backstop in case that
      # changes — ensures misaligned data is never silently used as valid.
      if (parsed$status == "col_shift") {
        parsed$mat <- matrix(2, nrow = length(batch_fragrances), ncol = n_questions)
        parsed$detail <- paste0(parsed$detail, "DATA DISCARDED (replaced with 2s). ")
      }

      results_matrix[start_i:end_i, ] <- parsed$mat
      n_success <- n_success + 1

      # Log any non-ok results
      if (parsed$status != "ok") {
        problem_log <- rbind(problem_log, data.frame(
          custom_id = custom_id, start_i = start_i, end_i = end_i,
          status = parsed$status, detail = parsed$detail,
          stringsAsFactors = FALSE
        ))
      }
    }, error = function(e) {
      cat(sprintf("  Parse/assign error [%s] (fragrances %d-%d): %s\n", custom_id, start_i, end_i, e$message))
      problem_log <<- rbind(problem_log, data.frame(
        custom_id = custom_id, start_i = start_i, end_i = end_i,
        status = "assign_error", detail = e$message,
        stringsAsFactors = FALSE
      ))
      n_fail <<- n_fail + 1
    })
  } else {
    n_fail <- n_fail + 1
    problem_log <- rbind(problem_log, data.frame(
      custom_id = custom_id, start_i = start_i, end_i = end_i,
      status = "no_text", detail = "No response text extracted",
      stringsAsFactors = FALSE
    ))
  }

  # Incremental checkpoint
  if (ri %% checkpoint_every == 0) {
    save_checkpoint()
    cat(sprintf("  [%d/%d] %d success, %d fail — checkpoint saved\n",
                ri, length(all_results), n_success, n_fail))
  }
}

# Final checkpoint
save_checkpoint()

cat(sprintf("\nParsing complete: %d success, %d fail out of %d results.\n",
            n_success, n_fail, length(all_results)))
cat(sprintf("Web search: %d/%d used it, %d did not.\n",
            n_web_search_used, n_success + n_fail, n_web_search_missing))

# ====================== PROBLEM REPORT & RETRY ======================

# row_unmatched = at least one returned row couldn't be tied to a batch fragrance, so
#   placement is uncertain — worth a retry. row_matched is benign (some fragrances just
#   weren't returned and were left as 2s), so it is NOT retried by default; add it
#   here if you'd rather re-query partial batches too.
retry_statuses <- c("col_shift", "row_unmatched", "error", "assign_error", "no_text")

if (nrow(problem_log) > 0) {
  PROBLEM_LOG_FILE <- gsub("\\.csv$", "_problems.csv", OUTPUT_FILE)
  write.csv(problem_log, PROBLEM_LOG_FILE, row.names = FALSE)
  cat(sprintf("\nProblem log saved to %s (%d entries)\n", PROBLEM_LOG_FILE, nrow(problem_log)))

  cat("\nProblem summary:\n")
  for (s in unique(problem_log$status)) {
    n <- sum(problem_log$status == s)
    cat(sprintf("  %s: %d requests\n", s, n))
  }

  needs_retry <- problem_log[problem_log$status %in% retry_statuses, ]

  if (nrow(needs_retry) > 0) {
    RETRY_JSONL_FILE <- gsub("\\.jsonl$", "_retry.jsonl", BATCH_JSONL_FILE)
    cat(sprintf("\n%d requests need retry (statuses: %s)\n",
                nrow(needs_retry), paste(unique(needs_retry$status), collapse=", ")))
    cat(sprintf("Generating retry JSONL: %s\n", RETRY_JSONL_FILE))

    con <- file(RETRY_JSONL_FILE, "w")
    for (row_i in seq_len(nrow(needs_retry))) {
      si <- needs_retry$start_i[row_i]
      ei <- needs_retry$end_i[row_i]
      batch_fragrances <- fragrances[si:ei]
      prompt <- create_prompt(batch_fragrances, questions)
      custom_id <- paste0("req_", si, "_", ei)
      req <- make_request(custom_id, prompt)
      writeLines(toJSON(req, auto_unbox = TRUE), con)
    }
    close(con)

    n_retry <- length(readLines(RETRY_JSONL_FILE))
    cat(sprintf("Retry JSONL created with %d requests.\n", n_retry))
    cat("To run the retry batch:\n")
    cat("  1. Set BATCH_JSONL_FILE to the retry file path\n")
    cat("  2. Set SKIP_JSONL_GENERATION <- TRUE\n")
    cat("  3. Delete or rename BATCH_ID_FILE so a new batch is created\n")
    cat("  4. Re-run the script\n")
    cat("  (The checkpoint will preserve all previously successful results.)\n")
  } else {
    cat("\nAll problems were minor (col_aligned or row_fix) — no retry needed.\n")
  }
} else {
  cat("\nNo problems detected — all results parsed cleanly!\n")
}

# ====================== FINAL SAVE ======================
results_df <- cbind(Fragrance = rownames(results_matrix), as.data.frame(results_matrix))
write.csv(results_df, OUTPUT_FILE, row.names = FALSE)

cat(sprintf("\nDone! Results in %s\n", OUTPUT_FILE))
cat(sprintf("Batch ID: %s (saved in %s)\n", batch_id, BATCH_ID_FILE))


#####################################################################################################
#####################################################################################################
# This concludes the fragrance data collection phase.
# The next phase runs a factor analysis on the data collected above.
#####################################################################################################
#####################################################################################################

library(psych)

fnames <- c("Resinous/Smoky","Animalic","Expensive","Intense/Powerful","Earthy",
            "Gourmand/Sweet","Masculine","Aquatic/Ozonic","Spicy/Woody",
            "Carnation/Aldehydic","Minty","Fruity","Vintage/Waxy","US_Mainstream")

getTheta <- function(mod, items) {
  sc <- scoreIrt(mod, items = items, cut = 0.3)
  sc <- sc[, grep("theta", colnames(sc))]
  colnames(sc) <- fnames
  sc
}

# Apply fixed reference-set parameters; clamp = FALSE lets you see values beyond the bounds
postProcess <- function(th, lo, hi, clamp = TRUE) {
  out <- th
  for (j in seq_along(lo)) {
    v <- th[, j]
    if (clamp) v <- pmin(pmax(v, lo[j]), hi[j])
    out[, j] <- round(-4 + 8 * (v - lo[j]) / (hi[j] - lo[j]), 2)
  }
  out
}

## ---- Reference (large) set ----
x <- read.csv("fragrance_data_raw_forPub.csv")
x[, c(25, 72)] <- 4 - x[, c(25, 72)]
items <- x[, -1]
items[items == 3] <- 2
items[items == 4] <- 3

mod    <- irt.fa(items, 14, rotate = "oblimin")
th_ref <- getTheta(mod, items)

# Estimate post-processing parameters ONCE, from the reference set
trim <- 0.001
lo <- apply(th_ref, 2, quantile, probs = trim,     na.rm = TRUE)
hi <- apply(th_ref, 2, quantile, probs = 1 - trim, na.rm = TRUE)
write.csv(data.frame(Factor = fnames, lo = lo, hi = hi),
          "fragrance_scaling_params.csv", row.names = FALSE)

ref_out <- data.frame(Fragrance = x[, 1], postProcess(th_ref, lo, hi), check.names = FALSE)
write.csv(ref_out, "fragrance_data_forPub.csv", na = "", row.names = FALSE)

## ---- Corrected 50 ----
y <- read.csv("famous_fragrance50_data_raw_corrected.csv")
y[, c(25, 72)] <- 4 - y[, c(25, 72)]
items2 <- y[, -1]
items2[items2 == 2] <- NA  
items2[items2 == 3] <- 2
items2[items2 == 4] <- 3

th_new  <- getTheta(mod, items2)
new_out <- data.frame(Fragrance = y[, 1], postProcess(th_new, lo, hi), check.names = FALSE)
write.csv(new_out, "famous_fragrance50_factor_scores_corrected.csv", na = "", row.names = FALSE)

#####################################################################################################
#####################################################################################################
# This concludes the factor score calculation phase.
# The next phase plots some of those factor scores for 50 well-known fragrances. These fragrances 
# were obtained by asking ChatGPT 6 to list the 50 most well-known fragrances ever made, where 
# "well-known" should combine total sales, frequency of appearance in media/film, recognizability 
# of the container and brand, and some others.
#####################################################################################################
#####################################################################################################

# ---- 1. Settings -------------------------------------------------------------
input_file <- "famous_fragrance50_factor_scores_corrected.csv"
output_dir <- NULL   # NULL = create a folder alongside the CSV

# Choose the x- and y-axis columns by their COLUMN NUMBERS in the CSV.
# Column 1 is the fragrance-name column and cannot be used as a numeric axis.
# Each row below creates one graph. You can change these numbers, add rows,
# or remove rows; graph titles, axis labels, and filenames are generated
# automatically from the selected column names.
#
# NOTE: In the finalized .png images, the top panel is the original (using data from "famous_fragrance50_factor_scores.csv") and the 
# bottom panel is the corrected (using data from "famous_fragrance50_factor_scores_corrected.csv") 
#
# These defaults reproduce the three graphs from the previous script.
plot_axes <- data.frame(
  x_col = c(8L, 2L, 7L),
  y_col = c(6L, 3L, 13L)
)

auto_install_packages <- TRUE

# Landscape, print-oriented output size
plot_width_in  <- 17
plot_height_in <- 11
png_dpi        <- 300

# Labels are at least 3x larger than in the original script
label_size_mm    <- 6.0
label_wrap_width <- 22L
point_size_mm    <- 1.8

# Match the 50-label examples you attached (about 25% of the 200 points).
# All points are still shown; only some labels are omitted. Increase this value
# if you want to attempt more labels.
label_keep_fraction <- 0.4
always_keep_extremes <- FALSE

# Optional priority labels: enter DATA ROW NUMBERS from the CSV (header excluded).
# These rows are ALWAYS labeled, in bold dark red, in every plot for which they
# have valid x/y scores. They are added on top of the ordinary fraction-based
# labels rather than using up any of the ordinary label quota.
# Examples:
#   highlight_rows <- c(5L, 17L, 42L, 88L)
#   highlight_rows <- integer(0)   # no highlighted labels
highlight_rows <- c(15L, 45L, 3L, 19L, 30L)
highlight_label_color <- "darkred"

repel_seed       <- 20260907L
repel_force      <- 0.9
repel_force_pull <- 0.025
repel_max_time   <- 30
repel_max_iter   <- 200000L

# Smaller paddings let labels pack more tightly when needed.
repel_box_padding   <- 0.10
repel_point_padding <- 0.06
segment_size_mm     <- 0.22

# Actively USE the blank vertical space rather than merely reserving it.
# Selected labels are split into an upper and lower half according to the
# y-coordinate of their points, then given opposite starting nudges. This keeps
# labels associated with relatively high-y points above labels associated with
# relatively low-y points, while encouraging ggrepel to fill the large blank
# bands above +4 and below -4.
vertical_label_nudge <- 1.75
vertical_split_gap   <- 0.15

# White space available to receive displaced labels.
x_margin_fraction <- 0.10
y_margin_fraction <- 0.35

show_in_plots_pane <- FALSE

# ---- 2. Locate the input and prepare packages --------------------------------
if (!file.exists(input_file)) {
  if (interactive()) {
    message("Input CSV not found at the configured path. Please select it.")
    input_file <- file.choose()
  } else {
    stop("Input CSV not found: ", input_file,
         "\nEdit input_file or run the script from the CSV's folder.",
         call. = FALSE)
  }
}
input_file <- normalizePath(input_file, winslash = "/", mustWork = TRUE)

packages <- c("ggplot2", "ggrepel")
missing_packages <- packages[!vapply(
  packages, requireNamespace, logical(1), quietly = TRUE
)]
if (length(missing_packages)) {
  if (!auto_install_packages) {
    stop("Install the following packages first: ",
         paste(missing_packages, collapse = ", "), call. = FALSE)
  }
  install.packages(missing_packages, repos = "https://cloud.r-project.org")
}
for (pkg in packages) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    stop("Package could not be loaded: ", pkg, call. = FALSE)
  }
}
suppressPackageStartupMessages({
  library(ggplot2)
  library(ggrepel)
})

# ---- 3. Read and validate the data --------------------------------------------
scores <- utils::read.csv(
  input_file,
  check.names = FALSE,
  stringsAsFactors = FALSE,
  colClasses = "character",
  fileEncoding = "UTF-8-BOM",
  na.strings = c("", "NA", "N/A", "NaN")
)

# Show the column-number/name mapping in the console so it is easy to choose
# different axes without having to count columns manually.
message("\nCSV columns available for plotting:")
message(paste(sprintf("  %d = %s", seq_along(names(scores)), names(scores)),
              collapse = "\n"))

# Validate the requested axis-column numbers.
if (!is.data.frame(plot_axes) ||
    !all(c("x_col", "y_col") %in% names(plot_axes)) ||
    nrow(plot_axes) < 1L) {
  stop("plot_axes must be a data frame with x_col and y_col columns and at least one row.",
       call. = FALSE)
}

axis_numbers <- unlist(plot_axes[c("x_col", "y_col")], use.names = FALSE)
if (any(!is.finite(axis_numbers)) ||
    any(axis_numbers != as.integer(axis_numbers))) {
  stop("Every x_col and y_col value in plot_axes must be a whole-number CSV column index.",
       call. = FALSE)
}
axis_numbers <- as.integer(axis_numbers)

if (any(axis_numbers < 2L) || any(axis_numbers > ncol(scores))) {
  stop(
    "Axis column numbers must be between 2 and ", ncol(scores),
    " (column 1 contains fragrance names).\nAvailable columns:\n",
    paste(sprintf("  %d = %s", seq_along(names(scores)), names(scores)),
          collapse = "\n"),
    call. = FALSE
  )
}

# Make filesystem-safe but still human-readable names automatically.
make_filename_piece <- function(x) {
  x <- gsub("[^[:alnum:]]+", "_", x)
  x <- gsub("^_+|_+$", "", x)
  x
}

plot_specs <- data.frame(
  x = names(scores)[as.integer(plot_axes$x_col)],
  y = names(scores)[as.integer(plot_axes$y_col)],
  stringsAsFactors = FALSE
)
plot_specs$filename <- sprintf(
  "%02d_%s_vs_%s_print",
  seq_len(nrow(plot_specs)),
  make_filename_piece(plot_specs$x),
  make_filename_piece(plot_specs$y)
)

required_columns <- unique(c(plot_specs$x, plot_specs$y))

if (nrow(scores) == 0L || ncol(scores) < 2L) {
  stop("The input must contain fragrance names and score columns.", call. = FALSE)
}
if (anyDuplicated(names(scores))) {
  stop("The CSV contains duplicate column names; resolve these first.",
       call. = FALSE)
}
absent_columns <- setdiff(required_columns, names(scores))
if (length(absent_columns)) {
  stop("Required columns not found: ", paste(absent_columns, collapse = ", "),
       "\nAvailable columns: ", paste(names(scores), collapse = ", "),
       call. = FALSE)
}
fragrance_names <- scores[[1L]]
if (any(is.na(fragrance_names) | !nzchar(trimws(fragrance_names)))) {
  stop("Every input row must have a fragrance name in column 1.", call. = FALSE)
}

# Validate highlighted data-row numbers. Row 1 here means the first fragrance
# row after the CSV header, matching the Data_row values written to output files.
if (is.null(highlight_rows)) highlight_rows <- integer(0)
if (!is.numeric(highlight_rows) ||
    any(!is.finite(highlight_rows)) ||
    any(highlight_rows != as.integer(highlight_rows))) {
  stop("highlight_rows must contain whole-number CSV data-row indices (or integer(0)).",
       call. = FALSE)
}
highlight_rows <- sort(unique(as.integer(highlight_rows)))
if (any(highlight_rows < 1L | highlight_rows > nrow(scores))) {
  stop("highlight_rows contains row number(s) outside the valid range 1 to ",
       nrow(scores), ".", call. = FALSE)
}

for (column in required_columns) {
  original <- scores[[column]]
  numeric_values <- suppressWarnings(as.numeric(original))
  invalid <- !is.na(original) & nzchar(trimws(original)) & is.na(numeric_values)
  if (any(invalid)) {
    stop("Non-numeric values in '", column, "' at data row(s): ",
         paste(which(invalid), collapse = ", "), call. = FALSE)
  }
  scores[[column]] <- numeric_values
}

all_scores <- unlist(scores[required_columns], use.names = FALSE)
finite_scores <- all_scores[is.finite(all_scores)]
if (!length(finite_scores)) {
  stop("No finite scores were found in the requested columns.", call. = FALSE)
}

score_bound <- max(4, ceiling(max(abs(finite_scores))))
x_plot_limits <- c(-score_bound, score_bound) +
  c(-1, 1) * (2 * score_bound * x_margin_fraction)
y_plot_limits <- c(-score_bound, score_bound) +
  c(-1, 1) * (2 * score_bound * y_margin_fraction)

axis_breaks <- pretty(c(-score_bound, score_bound), n = 8)
axis_breaks <- axis_breaks[axis_breaks >= -score_bound & axis_breaks <= score_bound]

wrap_label <- function(x) {
  vapply(x, function(z) {
    paste(strwrap(z, width = label_wrap_width), collapse = "\n")
  }, character(1), USE.NAMES = FALSE)
}

# Rank points for label retention: prioritize isolated and/or extreme points.
rank_for_labeling <- function(df) {
  xy <- as.matrix(df[, c("x", "y")])
  n <- nrow(xy)
  if (n <= 1L) return(seq_len(n))

  d <- as.matrix(stats::dist(xy))
  diag(d) <- Inf
  nearest_neighbor <- apply(d, 1, min)
  local_density <- rowSums(d < 0.75, na.rm = TRUE)
  extremeness <- sqrt(df$x^2 + df$y^2)

  score <-
    0.55 * as.numeric(scale(nearest_neighbor)) +
    0.30 * as.numeric(scale(extremeness)) -
    0.15 * as.numeric(scale(local_density))

  score[is.na(score)] <- 0
  order(score, decreasing = TRUE)
}

if (is.null(output_dir)) {
  output_dir <- file.path(dirname(input_file), "Fragrance_factor_plots_print")
}
if (!dir.exists(output_dir)) dir.create(output_dir, recursive = TRUE)
if (!dir.exists(output_dir)) {
  stop("Could not create output directory: ", output_dir, call. = FALSE)
}
output_dir <- normalizePath(output_dir, winslash = "/", mustWork = TRUE)

pdf_device <- if (capabilities("cairo")) grDevices::cairo_pdf else grDevices::pdf

# ---- 4. Create and export the three plots ------------------------------------
plots <- list()
summary_rows <- list()
excluded_rows <- list()
label_omissions <- list()

for (i in seq_len(nrow(plot_specs))) {
  x_col <- plot_specs$x[i]
  y_col <- plot_specs$y[i]
  plot_id <- plot_specs$filename[i]

  valid <- is.finite(scores[[x_col]]) & is.finite(scores[[y_col]])
  omitted <- which(!valid)
  excluded_rows[[i]] <- data.frame(
    Plot = rep(plot_id, length(omitted)),
    Data_row = omitted,
    Fragrance = fragrance_names[omitted],
    Missing_or_nonfinite_score = vapply(omitted, function(j) {
      paste(c(x_col, y_col)[c(!is.finite(scores[[x_col]][j]),
                            !is.finite(scores[[y_col]][j]))], collapse = "; ")
    }, character(1)),
    stringsAsFactors = FALSE
  )

  plot_data <- data.frame(
    x = scores[[x_col]][valid],
    y = scores[[y_col]][valid],
    label = fragrance_names[valid],
    wrapped_label = wrap_label(fragrance_names[valid]),
    original_row = which(valid),
    stringsAsFactors = FALSE
  )

  # Priority/highlight rows are handled separately from the ordinary label quota.
  # This means, for example, 10% ordinary labels + 10 highlighted rows can yield
  # more than 10% labeled points overall.
  plot_data$highlighted <- plot_data$original_row %in% highlight_rows
  highlight_idx <- which(plot_data$highlighted)

  n_to_label <- max(1L, ceiling(label_keep_fraction * nrow(plot_data)))
  ranked <- rank_for_labeling(plot_data)

  # Highlighted rows do NOT consume slots in the fraction-based selection.
  ranked_regular <- ranked[!ranked %in% highlight_idx]
  keep_idx <- ranked_regular[seq_len(min(n_to_label, length(ranked_regular)))]

  if (always_keep_extremes && nrow(plot_data) >= 1L) {
    extreme_idx <- unique(c(
      which.min(plot_data$x), which.max(plot_data$x),
      which.min(plot_data$y), which.max(plot_data$y)
    ))
    # A highlighted extreme is already guaranteed a label, so only non-highlighted
    # extremes need to be forced into the ordinary fraction-based set.
    regular_extreme_idx <- setdiff(extreme_idx, highlight_idx)
    keep_idx <- unique(c(keep_idx, regular_extreme_idx))
    if (length(keep_idx) > n_to_label) {
      non_extreme <- setdiff(ranked_regular, regular_extreme_idx)
      keep_idx <- unique(c(
        regular_extreme_idx,
        non_extreme[seq_len(max(0, n_to_label - length(regular_extreme_idx)))]
      ))
    }
  }
  keep_idx <- sort(unique(keep_idx))

  plot_data$labeled <- FALSE
  plot_data$labeled[unique(c(keep_idx, highlight_idx))] <- TRUE
  plot_data$label_color <- ifelse(plot_data$highlighted,
                                  highlight_label_color, "black")
  plot_data$label_face <- ifelse(plot_data$highlighted, "bold", "plain")

  # Divide selected labels approximately 50/50 between the lower and upper
  # portions of the panel. The split follows point y-rank, so leader lines are
  # generally short and vertical ordering remains intuitive.
  selected_idx <- which(plot_data$labeled)
  selected_by_y <- selected_idx[order(plot_data$y[selected_idx],
                                      plot_data$x[selected_idx])]
  n_lower <- floor(length(selected_by_y) / 2)
  lower_label_idx <- if (n_lower > 0L) selected_by_y[seq_len(n_lower)] else integer(0)
  upper_label_idx <- setdiff(selected_idx, lower_label_idx)

  upper_label_data <- plot_data[upper_label_idx, , drop = FALSE]
  lower_label_data <- plot_data[lower_label_idx, , drop = FALSE]

  omitted_labels_df <- plot_data[!plot_data$labeled, c("original_row", "label")]
  names(omitted_labels_df) <- c("Data_row", "Fragrance")
  omitted_labels_df$Plot <- rep(plot_id, nrow(omitted_labels_df))
  omitted_labels_df <- omitted_labels_df[, c("Plot", "Data_row", "Fragrance")]
  label_omissions[[i]] <- omitted_labels_df

  n_unique_points <- nrow(unique(plot_data[c("x", "y")]))
  summary_rows[[i]] <- data.frame(
    Plot = plot_id,
    X_axis = x_col,
    Y_axis = y_col,
    Input_rows = nrow(scores),
    Plotted_rows = nrow(plot_data),
    Labeled_rows = sum(plot_data$labeled),
    Highlighted_rows = sum(plot_data$highlighted),
    Unlabeled_rows = sum(!plot_data$labeled),
    Excluded_rows = length(omitted),
    Distinct_coordinates = n_unique_points,
    stringsAsFactors = FALSE
  )

  if (length(omitted)) {
    warning(sprintf(
      "%s: omitted %d row(s) with missing/nonfinite scores: %s",
      plot_id, length(omitted), paste(fragrance_names[omitted], collapse = "; ")
    ), call. = FALSE, immediate. = TRUE)
  }
  if (!nrow(plot_data)) {
    warning("Skipping ", plot_id, ": no complete score pairs.", call. = FALSE)
    next
  }

  p <- ggplot(plot_data, aes(x = x, y = y)) +
    geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.4,
               alpha = 0.45) +
    geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.4,
               alpha = 0.45) +
    # Subtle identity/reference line: points on this line have equal x and y scores.
    geom_abline(slope = 1, intercept = 0, linetype = "dotted",
                linewidth = 0.35, alpha = 0.35) +
    geom_point(size = point_size_mm, alpha = 0.8) +
    # Two repel layers are deliberate: upper-half labels are initialized upward
    # and lower-half labels downward, which makes practical use of the otherwise
    # empty top and bottom of a landscape page.
    geom_text_repel(
      data = upper_label_data,
      aes(label = wrapped_label, colour = label_color, fontface = label_face),
      size = label_size_mm,
      family = "sans",
      show.legend = FALSE,
      lineheight = 0.92,
      seed = repel_seed + i,
      max.overlaps = Inf,
      max.time = repel_max_time,
      max.iter = repel_max_iter,
      force = repel_force,
      force_pull = repel_force_pull,
      nudge_y = vertical_label_nudge,
      box.padding = repel_box_padding,
      point.padding = repel_point_padding,
      point.size = point_size_mm,
      min.segment.length = 0,
      segment.size = segment_size_mm,
      segment.alpha = 0.5,
      direction = "both",
      ylim = c(vertical_split_gap, y_plot_limits[2])
    ) +
    geom_text_repel(
      data = lower_label_data,
      aes(label = wrapped_label, colour = label_color, fontface = label_face),
      size = label_size_mm,
      family = "sans",
      show.legend = FALSE,
      lineheight = 0.92,
      seed = repel_seed + 1000L + i,
      max.overlaps = Inf,
      max.time = repel_max_time,
      max.iter = repel_max_iter,
      force = repel_force,
      force_pull = repel_force_pull,
      nudge_y = -vertical_label_nudge,
      box.padding = repel_box_padding,
      point.padding = repel_point_padding,
      point.size = point_size_mm,
      min.segment.length = 0,
      segment.size = segment_size_mm,
      segment.alpha = 0.5,
      direction = "both",
      ylim = c(y_plot_limits[1], -vertical_split_gap)
    ) +
    scale_colour_identity(guide = "none") +
    scale_x_continuous(breaks = axis_breaks) +
    scale_y_continuous(breaks = axis_breaks) +
    coord_cartesian(xlim = x_plot_limits, ylim = y_plot_limits, expand = FALSE) +
    labs(
      title = paste(x_col, "and", y_col),
      subtitle = sprintf(
        "%d of %d fragrances plotted; %d labels shown (%.0f%% of plotted fragrances)",
        nrow(plot_data), nrow(scores), sum(plot_data$labeled),
        100 * sum(plot_data$labeled) / nrow(plot_data)
      ),
      x = x_col,
      y = y_col,
      caption = paste(
        "All complete cases are plotted as points.",
        "Selected large-print labels are displaced vertically to use the full landscape panel.",
        "Fragrances with identical score pairs share a point."
      )
    ) +
    theme_bw(base_size = 16, base_family = "sans") +
    theme(
      panel.grid.minor = element_blank(),
      panel.grid.major = element_line(linewidth = 0.22),
      plot.title = element_text(size = 24, face = "bold", hjust = 0),
      plot.subtitle = element_text(size = 15, margin = margin(b = 14)),
      axis.title = element_text(size = 19, face = "bold"),
      axis.text = element_text(size = 14),
      plot.caption = element_text(size = 11, hjust = 0, margin = margin(t = 12)),
      plot.margin = margin(18, 26, 18, 18)
    )

  plots[[plot_id]] <- p
  message(sprintf(
    "Exporting %s: %d plotted fragrances; %d labels retained (%d highlighted); %d unlabeled.",
    plot_id, nrow(plot_data), sum(plot_data$labeled),
    sum(plot_data$highlighted), sum(!plot_data$labeled)
  ))

  ggsave(filename = file.path(output_dir, paste0(plot_id, ".pdf")),
         plot = p, device = pdf_device,
         width = plot_width_in, height = plot_height_in, units = "in")
  ggsave(filename = file.path(output_dir, paste0(plot_id, ".png")),
         plot = p, device = "png", dpi = png_dpi,
         width = plot_width_in, height = plot_height_in, units = "in")

  if (show_in_plots_pane && interactive()) print(p)
}

# ---- 5. Save accounting tables ------------------------------------------------
plot_summary <- do.call(rbind, summary_rows)
exclusions <- do.call(rbind, excluded_rows)
omitted_labels <- do.call(rbind, label_omissions)

utils::write.csv(plot_summary, file.path(output_dir, "plot_summary.csv"),
                 row.names = FALSE, fileEncoding = "UTF-8")
utils::write.csv(exclusions, file.path(output_dir, "excluded_rows.csv"),
                 row.names = FALSE, fileEncoding = "UTF-8")
utils::write.csv(omitted_labels, file.path(output_dir, "label_omissions.csv"),
                 row.names = FALSE, fileEncoding = "UTF-8")
writeLines(capture.output(sessionInfo()), file.path(output_dir, "session_info.txt"))

message("\nFinished. All output is in:\n", output_dir)
print(plot_summary, row.names = FALSE)

#########################################################################################################
#########################################################################################################
# This concludes the factor score plotting phase.
# The next phase runs a systematic check of the values obtained from the Grok batch API runs above. This 
# is arguably the focus of the present manuscript. Because the project is self-funded, I can't afford 
# to run the correction script on the full data set of 35,449 fragrances. I therefore ran the correction 
# script on the 200 random fragrances and the 200 most well-known fragrances. The correction script was 
# set to complete 20 sweeps (of all columns), checking 5 outliers (or missing values) per column per 
# sweep. This is 100 rows checked per column - exactly 50% of the 200-row data set. These 20 sweeps cost 
# about $350 (so $700 for both the 200 random and 200 well-known). Extrapolating, running this on the 
# 35449-fragrance data set would cost about $62k. Realistically, though, one wouldn't run the correction 
# script to cover 50% of his/her data. A more realistic use would be perhaps 5%, resulting in much lower 
# cost. The exact script below runs it on the 200 random, but the 200 well-known were also run to 
# produce data analyzed further below.
#########################################################################################################
#########################################################################################################

# ==============================================================================
# FragranceDive_correction_script.R
#
# Iterative multivariate data-quality checking for 0-4 Likert-type annotations
# of fragrances (perfumes, colognes, etc.). Adapted from the PharmaDive compound
# QC script; consumes the matrix produced by the FragranceDive batch script.
#
# For each question/column, the script:
#
#   1. Optionally prioritizes genuinely missing/invalid cells in the target
#      column before doing any multivariate modeling.
#   2. If the target column still has review budget after missing cells are
#      selected, predicts that column from all other annotation columns using an
#      out-of-bag random-forest model.
#   3. Scores non-missing rows as multivariate outliers when the observed value
#      has low predicted probability and/or is ordinally far from the model
#      expectation.
#   4. Sends the selected missing cells/outliers for that one column/question to
#      a web-search-enabled LLM API.
#   5. Applies corrections only when the LLM returns parseable JSON and meets
#      the configured confidence threshold.
#   6. Persists output/checkpoint/review-log progress at the configured cadence
#      (default: end of each sweep) instead of after every LLM batch.
#
# Required packages:
#   install.packages(c("data.table", "ranger", "httr", "jsonlite", "dotenv"))
#
# Expected input CSV (the FragranceDive batch output):
#   - Column 1: fragrance name, e.g. "Aventus (Creed)"
#   - Columns 2..N: Q1, Q2, ... with integer values 0,1,2,3,4
#   Because that file uses generic Q1..Qn headers rather than question text,
#   QUESTIONS_FILE is effectively REQUIRED here (see config below).
#
# API keys:
#   - OpenAI:    OPENAI_API_KEY
#   - Anthropic: ANTHROPIC_API_KEY
#   - xAI/Grok:  XAI_API_KEY
#
# Strong recommendation:
#   Run first with DRY_RUN <- TRUE and N_SWEEPS <- 1, OUTLIERS_PER_COL <- 1.
#
# ==============================================================================

# ----------------------------- USER CONFIG ------------------------------------

INPUT_FILE <- "fragrance_data_raw_200.csv"
OUTPUT_FILE <- "fragrance_data_raw_200_corrected.csv"
CHECKPOINT_FILE <- "TEMP_fragrances200_checkpoint.rds"
REVIEW_LOG_FILE <- "TEMP_fragrances200_change_summary.csv"
ERROR_LOG_FILE <- "TEMP_fragrances200_error_log.log"
OUTLIER_LOG_FILE <- "TEMP_fragrances200_outlier_log.csv"

# Optional .env file(s). The first one is the file name used in your existing Grok scripts.
# The script loads every file in this vector that exists in the current working directory.
# Set to NA_character_ to skip dotenv loading entirely.
ENV_FILE <- c("grok_api_key.env", ".env", "claude_api_key.env", "openai_api_key.env")

ID_COL <- 1L                           # "Fragrance" column written by the batch script
FIRST_QUESTION_COL <- 2L
LAST_QUESTION_COL <- NA_integer_       # NA means use ncol(input)

# The FragranceDive batch output has generic Q1..Qn headers, so the actual
# question text has to come from the same file the batch script read. Without
# this the LLM would be asked to audit a column called "Q7" with no idea what
# Q7 means. The file must have exactly one question per non-empty line, in the
# same order used to build the matrix.
QUESTIONS_FILE <- "TEMP_frag_Qs.txt"

TARGET_COLUMNS_TO_RUN <- NULL          # NULL = all question columns; otherwise names or 1-based CSV indices

# Iteration controls
N_SWEEPS <- 20L
OUTLIERS_PER_COL <- 5L                 # total cells to review per target column per sweep
LLM_BATCH_SIZE <- 5L                    # fragrances per one LLM request, for one question only
MAX_COLUMNS_PER_SWEEP <- Inf           # useful for testing
SLEEP_BETWEEN_API_CALLS <- 1.7
SLEEP_BETWEEN_COLUMNS <- 0.0

# Persistence controls
# Controls how often the three main run-state files are written:
#   OUTPUT_FILE, CHECKPOINT_FILE, and REVIEW_LOG_FILE.
# Defaulting to "sweep" reduces disk churn while still giving one durable
# checkpoint/review-log update after each complete sweep.
# Options:
#   "batch"  = old behavior, persist after every LLM batch
#   "column" = persist after each target question/column
#   "sweep"  = persist after each full sweep (recommended default)
#   "final"  = persist only at the end of the run
SAVE_PROGRESS_EVERY <- "sweep"

# Missing-value controls
# Reported study inputs had all collection-stage code-2 (Unknown/not-applicable)
# responses converted to blanks/NA before statistical analysis. The switches below
# are retained for generic/raw or future correction runs; they do not give original
# code-2 responses an analytical role in the reported results.
# TREAT_2_AS_MISSING applies only when an input file still contains literal 2s.
# Any newly returned Unknown value is converted to a blank before statistical analysis.
PRIORITIZE_MISSING_VALUES <- TRUE
MISSING_VALUE_STRINGS <- c("", "NA", "N/A", "NULL", "MISSING", ".", "?")
TREAT_2_AS_MISSING <- FALSE
MIN_CONFIDENCE_TO_APPLY_MISSING <- 0.70  # set NA_real_ to reuse MIN_CONFIDENCE_TO_APPLY
APPLY_ANY_VALID_RECOMMENDATION_TO_MISSING <- FALSE
RETRY_UNFILLED_MISSING_VALUES <- FALSE    # permanently skip missing cells left as NA?

# ---- Fragrance-specific: 2s that were never actually researched ----
# The FragranceDive batch script initialises its whole matrix to 2 and only
# overwrites cells it successfully parsed. Two situations therefore leave
# "2 = Unknown" in the CSV without the model ever having looked at the item:
#   (a) a batch whose columns could not be aligned ("col_shift") had its data
#       deliberately discarded and replaced with all 2s;
#   (b) a batch that errored, returned no text, or whose rows could not be
#       matched by name left some fragrances untouched at the initial 2.
# Those are not "researched Unknown" — they are unanswered. Left alone they
# would never be prioritized, because a literal 2 looks like a valid answer.
#
# Point this at the batch script's problem CSV (it writes OUTPUT_FILE with
# "_problems.csv" substituted for ".csv"). Rows whose fragrance index falls in
# a listed start_i:end_i range AND whose current value is 2 get treated as
# missing, so they are filled first. Set to NA_character_ to disable, or if
# the file no longer exists.
PROBLEM_LOG_FILE <- "TEMP_fragrances200_problems.csv"

# Which problem statuses imply the cells were never really answered. "row_matched"
# is deliberately excluded by default: in that status the fragrances that WERE
# returned are genuine, and the ones that were not are usually a small subset.
# Add it if you would rather re-check those whole batches.
PROBLEM_LOG_UNANSWERED_STATUSES <- c(
  "col_shift", "error", "assign_error", "no_text", "row_unmatched"
)

# Predictor-side missingness in the random-forest outlier model.
#
# IMPORTANT BUGFIX vs the PharmaDive version. That script set
#     MISSING_FILL_VALUE_FOR_MODEL <- NA_integer_
#     RANGER_NA_ACTION <- "na.learn"
# and passed na.action to ranger::ranger(). ranger() has no na.action argument.
# It discards it with a "Unused arguments: na.action" warning, and then errors
# with "Missing data in columns: ..." the moment any predictor holds an NA.
# That error is swallowed by the tryCatch around select_outliers_for_column()
# in main(), which logs a line and returns zero candidates -- so as soon as the
# data had any gaps, outlier detection quietly stopped running altogether and
# every sweep only ever filled missing cells. Verified against ranger 0.16.0.
#
# ranger genuinely cannot take NA predictors, so the fix is to impute, and to
# guarantee that nothing NA ever reaches it:
#   MODEL_NA_STRATEGY "constant"    fill with MISSING_FILL_VALUE_FOR_MODEL
#                     "column_mode" fill with each column's most common value
# "constant" with a fill of 2 is a software-level placeholder for blank predictor
# cells. It is created only for model fitting, is not an observed Unknown response,
# and is not retained in the analysis data.
MODEL_NA_STRATEGY <- "constant"
MISSING_FILL_VALUE_FOR_MODEL <- 2L

# Model controls
MODEL_TYPE <- "classification"         # "classification" recommended; "regression" also supported
N_TREES <- 500L
MTRY <- NULL                           # NULL = floor(sqrt(p))
MIN_NODE_SIZE_CLASSIFICATION <- 10L
MIN_NODE_SIZE_REGRESSION <- 20L
MIN_TRAIN_CLASSES <- 2L
SEED <- 20260612L

# Outlier scoring controls for classification
EPS_PROB <- 1e-6
ORDINAL_GAP_WEIGHT <- 1.25
CLASS_GAP_WEIGHT <- 0.75
MIN_OUTLIER_SCORE <- -Inf              # set e.g. 2.0 to avoid weak candidates

# Correction controls
APPLY_LLM_CORRECTIONS <- TRUE
DRY_RUN <- FALSE                       # TRUE = call LLM and log but do not edit data
MIN_CONFIDENCE_TO_APPLY <- 0.70
MANUAL_APPROVAL <- FALSE               # TRUE = ask interactively before applying each change
SKIP_PREVIOUSLY_REVIEWED <- TRUE       # skip cells already parsed/checked in REVIEW_LOG_FILE
RECHECK_IF_VALUE_CHANGED_SINCE_REVIEW <- TRUE

# LLM provider. Supported: "openai", "anthropic", "xai"
# Defaulting to xAI/Grok because your existing scripts use XAI_API_KEY in grok_api_key.env.
PROVIDER <- "xai"

# ---- Reasoning effort (shared knob, provider-specific field) ----
# "low" / "medium" / "high", plus "xhigh" on grok-4.6 and on OpenAI/Anthropic
# models that advertise it. Set to "" or NA to omit the field entirely and take
# the model default.
#
# This is a QC pass over a comparatively small number of cells per sweep, and
# every wrong correction is worse than a slow one, so "high" is a reasonable
# default here even though the bulk generation script runs at "medium".
# Drop to "medium" if cost becomes the binding constraint.
REASONING_EFFORT <- "high"

# OpenAI Responses API settings
# gpt-5.6 is the current general recommendation; gpt-5.6-terra / gpt-5.6-luna are
# cheaper, and gpt-5.5 still works if you have prompts tuned to it.
OPENAI_MODEL <- "gpt-5.6"
OPENAI_API_KEY_ENV <- "OPENAI_API_KEY"
OPENAI_API_URL <- "https://api.openai.com/v1/responses"
OPENAI_SEARCH_CONTEXT_SIZE <- "medium" # low, medium, high if supported by chosen model
OPENAI_FORCE_WEB_SEARCH <- TRUE         # uses tool_choice="required"

# Anthropic Messages API settings
# claude-sonnet-4-6 has been superseded; claude-sonnet-5 is the current
# speed/intelligence pick, claude-opus-5 is stronger and pricier.
ANTHROPIC_MODEL <- "claude-sonnet-5"
ANTHROPIC_API_KEY_ENV <- "ANTHROPIC_API_KEY"
ANTHROPIC_API_URL <- "https://api.anthropic.com/v1/messages"
ANTHROPIC_VERSION <- "2023-06-01"

# Three web search tool versions now exist:
#   web_search_20250305  basic search (default here)
#   web_search_20260209  adds dynamic filtering
#   web_search_20260318  adds response_inclusion
# The newer two are appealing (dynamic filtering cuts token spend on
# search-heavy work), but from web_search_20260209 onward `allowed_callers`
# defaults to ["code_execution_20260120"] rather than ["direct"], i.e. search
# runs from inside code execution. If you switch versions without thinking
# about that you may get a tool-config error or an unexpected bill. Setting
# ANTHROPIC_WEB_SEARCH_ALLOWED_CALLERS to "direct" pins the old behaviour.
ANTHROPIC_WEB_SEARCH_TOOL_TYPE <- "web_search_20250305"
ANTHROPIC_WEB_SEARCH_ALLOWED_CALLERS <- NA_character_  # e.g. "direct" on 20260209+

# Raised from 4096. On Sonnet 5 / Opus 5 adaptive thinking is on by default and
# max_tokens is a hard cap on thinking PLUS visible output, so at high effort a
# 4096 budget can be spent thinking and leave the JSON truncated mid-object.
# That surfaces as a parse_error rather than an obvious API failure, which is
# why call_anthropic() now logs stop_reason == "max_tokens" explicitly.
ANTHROPIC_MAX_TOKENS <- 16000L
ANTHROPIC_MAX_SEARCH_USES <- 5L

# xAI/Grok Responses API settings
# grok-4.6 (Aug 2026) is the current frontier model: 500k context, reasoning
# effort low/medium/high/xhigh, default high.
#
# NOTE: grok-4.6 is NOT Batch-API-enabled, which is fine here because this
# script makes synchronous calls, but it is exactly why the FragranceDive
# batch script must stay on grok-4.3. Do not "helpfully" sync the two model
# strings; the batch job will be rejected.
XAI_MODEL <- "grok-4.6"
XAI_API_KEY_ENV <- "XAI_API_KEY"
XAI_BASE_URL <- "https://api.x.ai/v1"
XAI_API_URL <- paste0(XAI_BASE_URL, "/responses")
XAI_FORCE_WEB_SEARCH <- TRUE            # uses tool_choice="required" if accepted by endpoint

# Shared LLM request settings
#
# Temperature is off by default and you probably want to leave it that way.
# Claude 4.7 and later (Sonnet 5, Opus 5, Fable 5.1) return HTTP 400 when
# temperature/top_p/top_k are set to a non-default value, and the Grok 4.x
# reasoning models reject several sampling parameters outright. On a reasoning
# model temperature 0 is not doing the determinism work it used to anyway.
SEND_TEMPERATURE <- FALSE
LLM_TEMPERATURE <- 0.0              # only used when SEND_TEMPERATURE is TRUE
LLM_MAX_RETRIES <- 5L
LLM_TIMEOUT_SECONDS <- 300L         # raised: high effort + web search is slow
PRINT_RAW_LLM_ON_PARSE_ERROR <- TRUE
RAW_RESPONSE_DIR <- "fragrance_qc_raw_llm_responses"

# Source-quality guidance included in prompts
SOURCE_GUIDANCE <- paste(
  "Prefer, in roughly this order: the perfume house's own product page and press materials;",
  "specialist trade and editorial coverage (Perfumer & Flavorist, Now Smell This, Persolaise,",
  "Bois de Jasmin, Cafleurebon); Michael Edwards' Fragrances of the World; authorised retailer",
  "listings (Sephora, Bloomingdale's, LuckyScent, Jovoy, Aus Liebe zum Duft) for launch year,",
  "concentration and availability; and IFRA/regulatory documents for materials questions.",
  "Fragrantica, Parfumo and Basenotes are useful and often the only source for note pyramids,",
  "launch years and perfumer credits, but they are user-contributed: treat their community",
  "rankings, longevity/sillage votes and reformulation claims as weak evidence, and cap confidence",
  "accordingly. Note that houses frequently do not disclose full compositions, so an absence of",
  "evidence about a specific material is not evidence of absence."
)

FRAGRANCE_INTERPRETATION_POLICY <- paste(
  "Identify the exact release before answering. Fragrance names collide constantly, so:",
  "(1) If the item is given as 'Name (House)', the parenthetical is the house and is",
  "authoritative for disambiguation. If no house is given and several houses have a fragrance",
  "of that name, use 2=Unknown unless one release is overwhelmingly the referent.",
  "(2) Distinguish the base fragrance from its flankers and its concentration variants: Eau de",
  "Cologne, Eau de Toilette, Eau de Parfum, Parfum/Extrait, Elixir, Intense, Absolu, and",
  "gendered counterparts ('pour Homme', 'for Her') are different products with different",
  "compositions. If the item names no concentration, answer for the original/best-known release",
  "and say so in the rationale rather than averaging across the line.",
  "(3) Where a fragrance has been reformulated, answer for the current widely-sold formulation",
  "unless the question explicitly concerns a vintage. Where vintage and current differ on the",
  "point being asked, prefer 1 or 3 over 0 or 4.",
  "(4) Discontinued, exclusive-boutique, and region-limited releases are frequently thinly",
  "documented; use 2 rather than inferring from a sibling fragrance in the same line."
)

# --------------------------- PACKAGE SETUP ------------------------------------

required_packages <- c("data.table", "ranger", "httr", "jsonlite", "dotenv")
missing_packages <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_packages) > 0) {
  stop(
    "Missing required packages: ", paste(missing_packages, collapse = ", "),
    "\nInstall with: install.packages(c(",
    paste(sprintf('"%s"', missing_packages), collapse = ", "), "))",
    call. = FALSE
  )
}

library(data.table)
library(ranger)
library(httr)
library(jsonlite)

load_env_files <- function(env_files) {
  if (is.null(env_files) || all(is.na(env_files))) return(character())
  env_files <- env_files[!is.na(env_files) & nzchar(env_files)]
  loaded <- character()
  for (env_file in env_files) {
    if (file.exists(env_file)) {
      dotenv::load_dot_env(env_file)
      loaded <- c(loaded, normalizePath(env_file, winslash = "/", mustWork = FALSE))
    }
  }
  loaded
}

LOADED_ENV_FILES <- load_env_files(ENV_FILE)

assert_api_key <- function(env_var, provider_name) {
  api_key <- Sys.getenv(env_var)
  if (!nzchar(api_key)) {
    loaded_msg <- if (length(LOADED_ENV_FILES) > 0) {
      paste(LOADED_ENV_FILES, collapse = ", ")
    } else {
      "none"
    }
    stop(
      sprintf(
        paste0(
          "Missing %s for %s. Current working directory: %s. Env files loaded: %s. ",
          "Put a line like %s=your_key_here in grok_api_key.env, ",
          "or set it with Sys.setenv(%s = 'your_key_here') before running."
        ),
        env_var, provider_name, getwd(), loaded_msg, env_var, env_var
      ),
      call. = FALSE
    )
  }
  api_key
}

`%||%` <- function(x, y) if (is.null(x)) y else x

# ------------------------------ LOGGING ---------------------------------------

log_error <- function(msg) {
  timestamp <- format(Sys.time(), "%Y-%m-%d %H:%M:%S")
  line <- sprintf("[%s] %s\n", timestamp, msg)
  cat(line)
  cat(line, file = ERROR_LOG_FILE, append = TRUE)
}

log_info <- function(...) {
  cat(sprintf(...), "\n")
}

safe_dir_create <- function(path) {
  if (!dir.exists(path)) dir.create(path, recursive = TRUE, showWarnings = FALSE)
}

# ---------------------------- DATA HELPERS ------------------------------------

detect_missing_likert <- function(x, treat_2_as_missing = TREAT_2_AS_MISSING) {
  # TRUE means "this input cell was genuinely missing or invalid" and should be
  # prioritized for LLM filling. Literal-2 handling here is generic raw-input
  # compatibility; reported analysis inputs had original 2s converted to blanks.
  s <- trimws(as.character(x))
  explicit_missing <- is.na(x) | toupper(s) %in% toupper(MISSING_VALUE_STRINGS)

  numeric_value <- suppressWarnings(as.numeric(s))
  integerish <- !is.na(numeric_value) & abs(numeric_value - round(numeric_value)) < 1e-8
  valid_likert <- integerish & as.integer(round(numeric_value)) %in% 0:4

  missing <- explicit_missing | !valid_likert
  if (isTRUE(treat_2_as_missing)) {
    missing <- missing | (valid_likert & as.integer(round(numeric_value)) == 2L)
  }
  missing
}

sanitize_likert <- function(x, treat_2_as_missing = TREAT_2_AS_MISSING) {
  missing <- detect_missing_likert(x, treat_2_as_missing = treat_2_as_missing)
  y <- suppressWarnings(as.integer(round(as.numeric(trimws(as.character(x))))))
  y[missing] <- NA_integer_
  y
}

column_mode <- function(v) {
  v <- v[!is.na(v)]
  if (length(v) == 0L) return(2L)   # nothing to learn from; fall back to Unknown
  tab <- table(v)
  as.integer(names(tab)[which.max(tab)])
}

make_model_matrix <- function(mat) {
  out <- mat
  strategy <- tolower(trimws(as.character(MODEL_NA_STRATEGY %||% "constant")[1]))

  if (identical(strategy, "column_mode")) {
    for (j in seq_len(ncol(out))) {
      na_j <- is.na(out[, j])
      if (any(na_j)) out[na_j, j] <- column_mode(out[, j])
    }
  } else {
    fill <- suppressWarnings(as.integer(MISSING_FILL_VALUE_FOR_MODEL[1]))
    if (is.na(fill)) fill <- 2L
    out[is.na(out)] <- fill
  }

  # Belt and braces. ranger errors outright on NA predictors, and that error
  # gets swallowed upstream into "zero candidates", so never let one through.
  if (anyNA(out)) {
    out[is.na(out)] <- 2L
  }

  storage.mode(out) <- "integer"
  out
}

model_missing_handling_description <- function() {
  strategy <- tolower(trimws(as.character(MODEL_NA_STRATEGY %||% "constant")[1]))
  if (identical(strategy, "column_mode")) {
    "predictor NAs imputed with each column's modal value"
  } else {
    sprintf("predictor NAs imputed as %s", as.character(MISSING_FILL_VALUE_FOR_MODEL[1]))
  }
}

format_value_for_log <- function(x) {
  if (length(x) == 0L || is.na(x)) return("missing")
  as.character(as.integer(x))
}

normalize_name <- function(x) {
  # BUGFIX carried over from the PharmaDive version: the original applied
  # gsub("[^a-z0-9]+", "") BEFORE tolower(), so every uppercase letter was
  # deleted rather than lowercased ("Aventus" -> "ventus"). Both sides of a
  # comparison were mangled identically, so exact echoes still matched and the
  # bug stayed invisible, but any case difference between the name we sent and
  # the name the model echoed back ("CREED" vs "Creed") silently failed to
  # match and the row was logged as unmatched. tolower() now runs first.
  tolower_x <- tolower(trimws(as.character(x)))
  gsub("[^a-z0-9]+", "", tolower_x)
}

normalize_name_base <- function(x) {
  # House-stripped form: "Aventus (Creed)" -> "aventus".
  # Only ever used as a last-resort fallback, and only when unambiguous on both
  # sides, because the parenthetical is exactly what separates "Aventus (Creed)"
  # from "Aventus (Mancera)".
  normalize_name(gsub("\\(.*?\\)", "", as.character(x)))
}

make_cell_key <- function(row_index, target_local_idx) {
  # Use column position, not column name, so duplicate question headers are
  # treated as separate cells/targets throughout review skipping and logging.
  paste0(row_index, "||target_local_", as.integer(target_local_idx))
}

normalize_review_log_schema <- function(x) {
  # CSV round-tripping can change inferred classes. In particular, fread()
  # reads ISO-like timestamps as POSIXct by default, while newly buffered rows
  # use character timestamps from format(Sys.time()). Normalize every review-log
  # piece before rbindlist() so persisted and in-memory rows always bind safely.
  dt <- data.table::copy(as.data.table(x))
  if (ncol(dt) == 0L) return(dt)

  character_cols <- c(
    "timestamp", "provider", "model", "target_col_name", "question_text",
    "fragrance", "candidate_type", "llm_parse_status", "rationale", "sources",
    "raw_response_file"
  )
  integer_cols <- c(
    "sweep", "target_local_idx", "target_csv_col_idx", "row_index",
    "original_value_at_review", "recommended_value", "final_value",
    "search_uses", "api_status", "pred_class"
  )
  numeric_cols <- c(
    "confidence", "pred_expected", "true_prob", "outlier_score", "model_oob_error"
  )
  logical_cols <- c("change_made", "needs_change_llm")

  for (nm in intersect(character_cols, names(dt))) {
    data.table::set(dt, j = nm, value = as.character(dt[[nm]]))
  }
  for (nm in intersect(integer_cols, names(dt))) {
    data.table::set(dt, j = nm, value = suppressWarnings(as.integer(dt[[nm]])))
  }
  for (nm in intersect(numeric_cols, names(dt))) {
    data.table::set(dt, j = nm, value = suppressWarnings(as.numeric(dt[[nm]])))
  }
  for (nm in intersect(logical_cols, names(dt))) {
    data.table::set(dt, j = nm, value = as.logical(dt[[nm]]))
  }
  dt
}

read_review_log_file <- function(path) {
  if (!file.exists(path)) return(data.table())
  out <- data.table::fread(
    path,
    encoding = "UTF-8",
    showProgress = FALSE,
    # Explicitly keep the timestamp as text. This prevents fread() from
    # assigning POSIXct while new in-memory rows remain character.
    colClasses = list(character = "timestamp")
  )
  normalize_review_log_schema(out)
}

read_existing_review_log <- function() {
  tryCatch(
    read_review_log_file(REVIEW_LOG_FILE),
    error = function(e) {
      log_error(sprintf("Could not read existing review log: %s", e$message))
      data.table()
    }
  )
}

append_csv <- function(dt, path) {
  # Used for auxiliary append-only files such as OUTLIER_LOG_FILE.
  # The main REVIEW_LOG_FILE is buffered and flushed by flush_review_log_buffer()
  # at the configured persistence cadence.
  dt <- as.data.table(dt)
  if (nrow(dt) == 0L) return(invisible(NULL))
  data.table::fwrite(dt, path, append = file.exists(path), col.names = !file.exists(path))
}

atomic_write_file <- function(path, writer_fun) {
  # Write to a temporary file in the same directory, then replace the target.
  # This prevents partially written checkpoint/output/review-log files if the
  # process is interrupted during file writing.
  dir <- dirname(path)
  if (is.na(dir) || !nzchar(dir)) dir <- "."
  if (!dir.exists(dir)) dir.create(dir, recursive = TRUE, showWarnings = FALSE)

  tmp <- tempfile(pattern = paste0(".", basename(path), ".tmp."), tmpdir = dir)
  backup <- NA_character_
  committed <- FALSE

  on.exit({
    if (file.exists(tmp)) unlink(tmp)
    if (committed && !is.na(backup) && file.exists(backup)) unlink(backup)
  }, add = TRUE)

  writer_fun(tmp)
  if (!file.exists(tmp)) stop("Atomic writer did not create temporary file for ", path, call. = FALSE)

  if (file.exists(path)) {
    backup <- tempfile(pattern = paste0(".", basename(path), ".bak."), tmpdir = dir)
    if (!file.rename(path, backup)) {
      stop("Could not move existing file aside before replacing: ", path, call. = FALSE)
    }
    if (!file.rename(tmp, path)) {
      restored <- file.rename(backup, path)
      stop(
        "Could not replace ", path,
        if (restored) "; restored previous file." else "; previous file is at backup path: ",
        if (restored) "" else backup,
        call. = FALSE
      )
    }
    committed <- TRUE
  } else if (!file.rename(tmp, path)) {
    stop("Could not move temporary file into place: ", path, call. = FALSE)
  } else {
    committed <- TRUE
  }

  invisible(TRUE)
}

write_current_output <- function(df) {
  # Preserve duplicate column headers exactly and keep still-missing cells blank.
  atomic_write_file(OUTPUT_FILE, function(path) {
    write.csv(df, path, row.names = FALSE, na = "")
  })
}

save_checkpoint <- function(df, sweep, target_local_idx, extra = list()) {
  cp <- list(
    df = df,
    sweep = sweep,
    target_local_idx = target_local_idx,
    timestamp = Sys.time(),
    config = list(
      input_file = INPUT_FILE,
      output_file = OUTPUT_FILE,
      checkpoint_file = CHECKPOINT_FILE,
      review_log_file = REVIEW_LOG_FILE,
      save_progress_every = SAVE_PROGRESS_EVERY,
      id_col = ID_COL,
      first_question_col = FIRST_QUESTION_COL,
      last_question_col = LAST_QUESTION_COL,
      provider = PROVIDER,
      model = current_llm_model(),
      reasoning_effort = effort_value(),
      send_temperature = SEND_TEMPERATURE,
      model_type = MODEL_TYPE,
      n_sweeps = N_SWEEPS,
      outliers_per_col = OUTLIERS_PER_COL,
      llm_batch_size = LLM_BATCH_SIZE,
      apply_llm_corrections = APPLY_LLM_CORRECTIONS,
      dry_run = DRY_RUN,
      prioritize_missing_values = PRIORITIZE_MISSING_VALUES,
      min_confidence_to_apply_missing = MIN_CONFIDENCE_TO_APPLY_MISSING,
      treat_2_as_missing = TREAT_2_AS_MISSING,
      problem_log_file = PROBLEM_LOG_FILE,
      missing_fill_value_for_model = MISSING_FILL_VALUE_FOR_MODEL,
      model_na_strategy = MODEL_NA_STRATEGY
    ),
    extra = extra
  )
  atomic_write_file(CHECKPOINT_FILE, function(path) saveRDS(cp, path))
  write_current_output(df)
}

combine_review_logs <- function(existing_log, buffered_rows = list()) {
  pieces <- list()
  if (!is.null(existing_log) && nrow(as.data.table(existing_log)) > 0L) {
    pieces[[length(pieces) + 1L]] <- existing_log
  }
  if (length(buffered_rows) > 0L) {
    pieces <- c(pieces, buffered_rows)
  }
  pieces <- lapply(pieces, normalize_review_log_schema)
  pieces <- pieces[vapply(pieces, nrow, integer(1)) > 0L]
  if (length(pieces) == 0L) return(data.table())
  rbindlist(pieces, use.names = TRUE, fill = TRUE)
}

flush_review_log_buffer <- function(buffered_rows, path = REVIEW_LOG_FILE) {
  new_rows <- combine_review_logs(data.table(), buffered_rows)
  if (nrow(new_rows) == 0L) return(0L)

  if (file.exists(path)) {
    existing <- tryCatch(
      read_review_log_file(path),
      error = function(e) {
        stop(
          "Could not read existing review log before appending buffered rows; refusing to overwrite ",
          path, ": ", e$message,
          call. = FALSE
        )
      }
    )
    out <- combine_review_logs(existing, list(new_rows))
  } else {
    out <- new_rows
  }

  atomic_write_file(path, function(tmp) {
    data.table::fwrite(out, tmp)
  })
  nrow(new_rows)
}

PERSISTENCE_LEVELS <- c(batch = 1L, column = 2L, sweep = 3L, final = 4L)

validate_persistence_config <- function() {
  cadence <- tolower(trimws(SAVE_PROGRESS_EVERY))
  if (!cadence %in% names(PERSISTENCE_LEVELS)) {
    stop(
      "SAVE_PROGRESS_EVERY must be one of: ", paste(names(PERSISTENCE_LEVELS), collapse = ", "),
      call. = FALSE
    )
  }
  invisible(cadence)
}

validate_request_config <- function() {
  # Catch a typo'd effort level here rather than as an opaque 400 twenty
  # minutes into a run. Not an exhaustive whitelist per model: xhigh is
  # grok-4.6+ / GPT-5.x / Opus-class only, and max is narrower still, so an
  # unsupported-but-spelled-correctly value still relies on the API to reject
  # it (and then on the 400 fallback to drop it and continue).
  eff <- effort_value()
  known <- c("none", "minimal", "low", "medium", "high", "xhigh", "max")
  if (!is.na(eff) && !eff %in% known) {
    stop("REASONING_EFFORT must be one of: ", paste(known, collapse = ", "),
         ", or \"\" to use the model default. Got: ", eff, call. = FALSE)
  }
  strategy <- tolower(trimws(as.character(MODEL_NA_STRATEGY %||% "constant")[1]))
  if (!strategy %in% c("constant", "column_mode")) {
    stop("MODEL_NA_STRATEGY must be \"constant\" or \"column_mode\"", call. = FALSE)
  }
  invisible(TRUE)
}

should_persist_after <- function(event) {
  cadence <- validate_persistence_config()
  event <- tolower(trimws(event))
  if (!event %in% names(PERSISTENCE_LEVELS)) stop("Unknown persistence event: ", event, call. = FALSE)
  unname(PERSISTENCE_LEVELS[event]) >= unname(PERSISTENCE_LEVELS[cadence])
}

persist_progress <- function(df, sweep, target_local_idx, buffered_review_rows = list(), extra = list()) {
  save_checkpoint(df, sweep, target_local_idx, extra = extra)
  n_review_rows <- flush_review_log_buffer(buffered_review_rows, REVIEW_LOG_FILE)
  log_info(
    "Persisted progress after %s: %s, %s, and %d buffered review-log row(s).",
    extra$persistence_event %||% "manual checkpoint",
    OUTPUT_FILE,
    CHECKPOINT_FILE,
    n_review_rows
  )
  invisible(n_review_rows)
}

load_input_or_checkpoint <- function() {
  if (file.exists(CHECKPOINT_FILE)) {
    log_info("Loading checkpoint: %s", CHECKPOINT_FILE)
    cp <- readRDS(CHECKPOINT_FILE)
    df <- cp$df
    attr(df, "loaded_from_checkpoint") <- TRUE
    attr(df, "checkpoint_sweep") <- suppressWarnings(as.integer(cp$sweep %||% 0L))
    attr(df, "checkpoint_target_local_idx") <- suppressWarnings(as.integer(cp$target_local_idx %||% NA_integer_))
    attr(df, "checkpoint_persistence_event") <- as.character(cp$extra$persistence_event %||% NA_character_)
    return(df)
  }
  if (!file.exists(INPUT_FILE)) stop("INPUT_FILE not found: ", INPUT_FILE, call. = FALSE)
  log_info("Loading input: %s", INPUT_FILE)
  # Use base read.csv(check.names = FALSE) rather than fread so duplicate
  # question headers are preserved exactly. Duplicate headers are allowed and
  # treated as distinct columns by their numeric CSV/local column positions.
  df <- read.csv(
    INPUT_FILE,
    stringsAsFactors = FALSE,
    check.names = FALSE,
    fileEncoding = "UTF-8",
    na.strings = character()
  )
  attr(df, "loaded_from_checkpoint") <- FALSE
  attr(df, "checkpoint_sweep") <- 0L
  attr(df, "checkpoint_target_local_idx") <- NA_integer_
  attr(df, "checkpoint_persistence_event") <- NA_character_
  df
}

resolve_question_columns <- function(df) {
  last_col <- if (is.na(LAST_QUESTION_COL)) ncol(df) else LAST_QUESTION_COL
  if (last_col > ncol(df)) stop("LAST_QUESTION_COL exceeds ncol(df)", call. = FALSE)
  if (FIRST_QUESTION_COL < 1L || FIRST_QUESTION_COL > last_col) {
    stop("Invalid FIRST_QUESTION_COL/LAST_QUESTION_COL", call. = FALSE)
  }
  qcols <- seq.int(FIRST_QUESTION_COL, last_col)
  if (ID_COL %in% qcols) stop("ID_COL must not be inside question columns", call. = FALSE)
  qcols
}

resolve_target_local_indices <- function(df, qcols) {
  if (is.null(TARGET_COLUMNS_TO_RUN)) return(seq_along(qcols))

  if (is.character(TARGET_COLUMNS_TO_RUN)) {
    csv_idx <- integer(length(TARGET_COLUMNS_TO_RUN))
    for (i in seq_along(TARGET_COLUMNS_TO_RUN)) {
      hits <- which(names(df) == TARGET_COLUMNS_TO_RUN[i])
      if (length(hits) == 0L) stop("Some TARGET_COLUMNS_TO_RUN names were not found in the CSV", call. = FALSE)
      if (length(hits) > 1L) {
        stop(sprintf(
          "TARGET_COLUMNS_TO_RUN name '%s' matches %d duplicate columns; use numeric CSV column indices for duplicate-named targets.",
          TARGET_COLUMNS_TO_RUN[i], length(hits)
        ), call. = FALSE)
      }
      csv_idx[i] <- hits[1]
    }
  } else {
    csv_idx <- as.integer(TARGET_COLUMNS_TO_RUN)
  }

  local_idx <- match(csv_idx, qcols)
  if (anyNA(local_idx)) {
    stop("TARGET_COLUMNS_TO_RUN must be within FIRST_QUESTION_COL:LAST_QUESTION_COL", call. = FALSE)
  }
  local_idx
}

load_questions <- function(df, qcols) {
  if (!is.na(QUESTIONS_FILE) && nzchar(QUESTIONS_FILE)) {
    if (!file.exists(QUESTIONS_FILE)) stop("QUESTIONS_FILE not found: ", QUESTIONS_FILE, call. = FALSE)
    qs <- readLines(QUESTIONS_FILE, warn = FALSE, encoding = "UTF-8")
    qs <- qs[nzchar(trimws(qs))]
    if (length(qs) != length(qcols)) {
      stop(sprintf(
        "QUESTIONS_FILE has %d non-empty lines but there are %d question columns",
        length(qs), length(qcols)
      ), call. = FALSE)
    }
    return(qs)
  }
  names(df)[qcols]
}

coerce_question_columns <- function(df, qcols, treat_2_as_missing = TREAT_2_AS_MISSING) {
  for (j in qcols) df[[j]] <- sanitize_likert(df[[j]], treat_2_as_missing = treat_2_as_missing)
  df
}

# ---- Batch problem log -> "these 2s were never researched" -------------------

read_problem_log_rows <- function() {
  # Returns a 1-based vector of CSV data-row indices whose values the batch
  # script never actually established, or integer(0).
  #
  # The batch script's custom_ids are "req_<start_i>_<end_i>" where start_i and
  # end_i index the `fragrances` vector, and it writes the output CSV in that
  # same order (rownames(results_matrix) <- fragrances). So start_i:end_i maps
  # straight onto data rows of INPUT_FILE. This breaks if you reorder or filter
  # rows between the two scripts -- don't.
  if (is.null(PROBLEM_LOG_FILE) || length(PROBLEM_LOG_FILE) == 0L) return(integer(0))
  if (is.na(PROBLEM_LOG_FILE) || !nzchar(PROBLEM_LOG_FILE)) return(integer(0))
  if (!file.exists(PROBLEM_LOG_FILE)) {
    log_info("PROBLEM_LOG_FILE not found (%s); skipping unanswered-2 detection.", PROBLEM_LOG_FILE)
    return(integer(0))
  }

  pl <- tryCatch(
    utils::read.csv(PROBLEM_LOG_FILE, stringsAsFactors = FALSE, fileEncoding = "UTF-8"),
    error = function(e) {
      log_error(sprintf("Could not read PROBLEM_LOG_FILE: %s", e$message))
      NULL
    }
  )
  if (is.null(pl) || nrow(pl) == 0L) return(integer(0))
  if (!all(c("start_i", "end_i", "status") %in% names(pl))) {
    log_error("PROBLEM_LOG_FILE lacks start_i/end_i/status columns; skipping unanswered-2 detection.")
    return(integer(0))
  }

  keep <- trimws(as.character(pl$status)) %in% PROBLEM_LOG_UNANSWERED_STATUSES
  pl <- pl[keep, , drop = FALSE]
  if (nrow(pl) == 0L) return(integer(0))

  starts <- suppressWarnings(as.integer(pl$start_i))
  ends   <- suppressWarnings(as.integer(pl$end_i))
  idx <- unlist(
    Map(function(a, b) if (is.na(a) || is.na(b) || b < a) integer(0) else seq.int(a, b), starts, ends),
    use.names = FALSE
  )

  sort(unique(as.integer(idx)))
}

mark_unanswered_twos_as_missing <- function(df, qcols) {
  rows <- read_problem_log_rows()
  if (length(rows) == 0L) return(df)

  rows <- rows[rows >= 1L & rows <= nrow(df)]
  if (length(rows) == 0L) {
    log_error("Problem-log row ranges fall entirely outside the input CSV; row order may have changed. Skipping.")
    return(df)
  }

  n_marked <- 0L
  for (j in qcols) {
    col <- df[[j]]
    hit <- rows[!is.na(col[rows]) & col[rows] == 2L]
    if (length(hit) > 0L) {
      col[hit] <- NA_integer_
      df[[j]] <- col
      n_marked <- n_marked + length(hit)
    }
  }

  log_info(
    "Problem log: %d fragrance row(s) from failed batches; reclassified %d literal 2(s) as unanswered so they get filled first.",
    length(rows), n_marked
  )
  df
}

# ------------------------- OUTLIER MODELING -----------------------------------

fit_score_classification <- function(mat, target_local_idx, q_names) {
  n <- nrow(mat)
  q <- ncol(mat)
  y <- mat[, target_local_idx]
  classes_present <- sort(unique(y))

  if (length(classes_present) < MIN_TRAIN_CLASSES) {
    log_error(sprintf(
      "Skipping %s: only %d class(es) present",
      q_names[target_local_idx], length(classes_present)
    ))
    return(data.table())
  }

  feature_idx <- setdiff(seq_len(q), target_local_idx)
  x <- as.data.frame(mat[, feature_idx, drop = FALSE])
  names(x) <- paste0("V", feature_idx)

  train_df <- cbind(.target = factor(y, levels = as.character(classes_present)), x)
  mtry_val <- if (is.null(MTRY)) max(1L, floor(sqrt(ncol(x)))) else MTRY

  fit <- ranger::ranger(
    formula = .target ~ .,
    data = train_df,
    probability = TRUE,
    num.trees = N_TREES,
    mtry = mtry_val,
    min.node.size = MIN_NODE_SIZE_CLASSIFICATION,
    seed = SEED + target_local_idx,
    classification = TRUE,
    oob.error = TRUE,
    verbose = FALSE
  )

  probs_raw <- fit$predictions
  prob_full <- matrix(0, nrow = n, ncol = 5L)
  colnames(prob_full) <- as.character(0:4)
  if (!is.null(colnames(probs_raw))) {
    common <- intersect(colnames(prob_full), colnames(probs_raw))
    prob_full[, common] <- probs_raw[, common, drop = FALSE]
  } else if (ncol(probs_raw) == 5L) {
    prob_full <- probs_raw
    colnames(prob_full) <- as.character(0:4)
  } else {
    stop("Unexpected probability matrix returned by ranger")
  }

  prob_full[is.na(prob_full)] <- 0
  prob_full <- pmax(prob_full, EPS_PROB)
  prob_full <- prob_full / rowSums(prob_full)

  values <- 0:4
  expected <- as.numeric(prob_full %*% values)
  pred_value <- values[max.col(prob_full, ties.method = "first")]
  true_prob <- prob_full[cbind(seq_len(n), y + 1L)]
  ordinal_gap <- abs(y - expected)
  class_gap <- abs(y - pred_value)
  score <- -log(pmax(true_prob, EPS_PROB)) + ORDINAL_GAP_WEIGHT * ordinal_gap + CLASS_GAP_WEIGHT * class_gap

  data.table(
    row_index = seq_len(n),
    target_value = y,
    pred_expected = expected,
    pred_class = pred_value,
    true_prob = true_prob,
    outlier_score = score,
    model_oob_error = fit$prediction.error %||% NA_real_
  )
}

fit_score_regression <- function(mat, target_local_idx, q_names) {
  n <- nrow(mat)
  q <- ncol(mat)
  y <- as.numeric(mat[, target_local_idx])
  if (length(unique(y)) < MIN_TRAIN_CLASSES) {
    log_error(sprintf(
      "Skipping %s: only %d unique value(s) present",
      q_names[target_local_idx], length(unique(y))
    ))
    return(data.table())
  }

  feature_idx <- setdiff(seq_len(q), target_local_idx)
  x <- as.data.frame(mat[, feature_idx, drop = FALSE])
  names(x) <- paste0("V", feature_idx)
  train_df <- cbind(.target = y, x)
  mtry_val <- if (is.null(MTRY)) max(1L, floor(sqrt(ncol(x)))) else MTRY

  fit <- ranger::ranger(
    formula = .target ~ .,
    data = train_df,
    num.trees = N_TREES,
    mtry = mtry_val,
    min.node.size = MIN_NODE_SIZE_REGRESSION,
    seed = SEED + target_local_idx,
    oob.error = TRUE,
    verbose = FALSE
  )

  pred <- pmin(4, pmax(0, fit$predictions))
  resid <- y - pred
  scale <- stats::mad(resid, constant = 1.4826, na.rm = TRUE)
  if (is.na(scale) || scale < 0.25) scale <- stats::sd(resid, na.rm = TRUE)
  if (is.na(scale) || scale < 0.25) scale <- 0.25
  score <- abs(resid) / scale

  data.table(
    row_index = seq_len(n),
    target_value = as.integer(y),
    pred_expected = pred,
    pred_class = as.integer(round(pred)),
    true_prob = NA_real_,
    outlier_score = score,
    model_oob_error = fit$prediction.error %||% NA_real_
  )
}

score_target_column <- function(mat, target_local_idx, q_names) {
  if (MODEL_TYPE == "classification") {
    fit_score_classification(mat, target_local_idx, q_names)
  } else if (MODEL_TYPE == "regression") {
    fit_score_regression(mat, target_local_idx, q_names)
  } else {
    stop("MODEL_TYPE must be 'classification' or 'regression'", call. = FALSE)
  }
}

reviewed_keys_for_current_data <- function(review_log, current_mat, q_names) {
  if (!SKIP_PREVIOUSLY_REVIEWED || nrow(review_log) == 0L) return(character())
  needed <- c("row_index", "target_local_idx", "llm_parse_status", "final_value")
  if (!all(needed %in% names(review_log))) {
    # Old logs that identified cells only by duplicated column names are unsafe
    # for duplicate-header inputs, so ignore them rather than accidentally skip
    # the wrong duplicate column.
    return(character())
  }

  ok <- review_log[llm_parse_status == "parsed"]
  if (nrow(ok) == 0L) return(character())

  if (RECHECK_IF_VALUE_CHANGED_SINCE_REVIEW) {
    keep <- logical(nrow(ok))
    for (i in seq_len(nrow(ok))) {
      local_idx <- suppressWarnings(as.integer(ok$target_local_idx[i]))
      row_idx <- suppressWarnings(as.integer(ok$row_index[i]))
      if (!is.na(local_idx) && !is.na(row_idx) &&
          local_idx >= 1L && local_idx <= ncol(current_mat) &&
          row_idx >= 1L && row_idx <= nrow(current_mat)) {
        current_value <- suppressWarnings(as.integer(current_mat[row_idx, local_idx]))
        final_value <- suppressWarnings(as.integer(ok$final_value[i]))
        if (isTRUE(RETRY_UNFILLED_MISSING_VALUES) && is.na(current_value) && is.na(final_value)) {
          keep[i] <- FALSE
        } else {
          keep[i] <- identical(current_value, final_value)
        }
      }
    }
    ok <- ok[keep]
  }

  unique(make_cell_key(ok$row_index, ok$target_local_idx))
}

select_missing_for_column <- function(df, mat, qcols, q_names, target_local_idx, review_log) {
  target_col_name <- q_names[target_local_idx]
  missing_rows <- which(is.na(mat[, target_local_idx]))
  if (length(missing_rows) == 0L) return(data.table())

  out <- data.table(
    row_index = missing_rows,
    target_value = NA_integer_,
    pred_expected = NA_real_,
    pred_class = NA_integer_,
    true_prob = NA_real_,
    outlier_score = Inf,
    model_oob_error = NA_real_,
    fragrance = as.character(df[[ID_COL]])[missing_rows],
    target_local_idx = target_local_idx,
    target_csv_col_idx = qcols[target_local_idx],
    target_col_name = target_col_name,
    cell_key = make_cell_key(missing_rows, target_local_idx),
    candidate_type = "missing"
  )

  already <- reviewed_keys_for_current_data(review_log, mat, q_names)
  if (length(already) > 0L) out <- out[!cell_key %in% already]
  out
}

select_outliers_for_column <- function(df, mat, qcols, q_names, target_local_idx, review_log) {
  target_col_name <- q_names[target_local_idx]
  budget <- as.integer(OUTLIERS_PER_COL)
  if (is.na(budget) || budget <= 0L) return(data.table())

  missing_candidates <- data.table()
  if (isTRUE(PRIORITIZE_MISSING_VALUES)) {
    missing_candidates <- select_missing_for_column(df, mat, qcols, q_names, target_local_idx, review_log)
    if (nrow(missing_candidates) >= budget) {
      log_info(
        "Found %d unchecked missing/invalid value(s); selecting %d and skipping outlier model for this column.",
        nrow(missing_candidates), budget
      )
      return(head(missing_candidates, budget))
    }
  }

  remaining_budget <- budget - nrow(missing_candidates)
  outlier_candidates <- data.table()

  if (remaining_budget > 0L) {
    model_mat <- make_model_matrix(mat)
    non_missing_target_rows <- which(!is.na(mat[, target_local_idx]))

    # Train/score the outlier model only on rows whose target value is known.
    # Missing feature values in other columns are either preserved for ranger's
    # internal missing-value handling or constant-imputed, depending on
    # MISSING_FILL_VALUE_FOR_MODEL.
    if (length(non_missing_target_rows) > 0L) {
      scores <- score_target_column(model_mat[non_missing_target_rows, , drop = FALSE], target_local_idx, q_names)
      if (nrow(scores) > 0L) {
        scores[, row_index := non_missing_target_rows[row_index]]
        scores[, `:=`(
          fragrance = as.character(df[[ID_COL]])[row_index],
          target_local_idx = target_local_idx,
          target_csv_col_idx = qcols[target_local_idx],
          target_col_name = target_col_name,
          cell_key = make_cell_key(row_index, target_local_idx),
          candidate_type = "outlier"
        )]

        already <- reviewed_keys_for_current_data(review_log, mat, q_names)
        if (length(already) > 0L) scores <- scores[!cell_key %in% already]
        scores <- scores[outlier_score >= MIN_OUTLIER_SCORE]
        data.table::setorder(scores, -outlier_score, true_prob)
        outlier_candidates <- head(scores, remaining_budget)
      }
    }
  }

  out <- rbindlist(list(missing_candidates, outlier_candidates), fill = TRUE)

  if (nrow(out) > 0L) {
    append_csv(out[, .(
      timestamp = format(Sys.time(), "%Y-%m-%d %H:%M:%S"),
      sweep = NA_integer_,
      target_local_idx,
      target_csv_col_idx,
      target_col_name,
      row_index,
      fragrance,
      candidate_type,
      target_value,
      pred_expected,
      pred_class,
      true_prob,
      outlier_score,
      model_oob_error
    )], OUTLIER_LOG_FILE)
  }
  out
}

# ------------------------------ PROMPTS ---------------------------------------

SYSTEM_PROMPT <- paste(
  "You are a careful perfumery data-quality auditor with deep knowledge of fragrance houses,",
  "perfumers, olfactive families, note pyramids, concentrations, flankers and reformulations.",
  "You check one annotation question at a time for a small list of fragrances.",
  "Many of these releases are obscure, niche, discontinued or region-limited, so search the web",
  "before answering rather than relying on recall; a confident wrong answer is the worst outcome.",
  "Do not defer to the machine-learning prediction; it is only a triage signal.",
  "Return only valid JSON. Do not return markdown, prose outside JSON, apologies, or code fences.",
  "Use the exact 0-4 scale requested by the user. If reliable evidence is insufficient or the named",
  "fragrance cannot be pinned to a specific release, use 2=Unknown.",
  sep = "\n"
)

create_review_prompt <- function(question_text, target_col_name, candidates) {
  rows <- lapply(seq_len(nrow(candidates)), function(i) {
    list(
      fragrance = as.character(candidates$fragrance[i]),
      row_index = as.integer(candidates$row_index[i]),
      current_value = if (is.na(candidates$target_value[i])) NULL else as.integer(candidates$target_value[i]),
      candidate_type = as.character(candidates$candidate_type[i] %||% "outlier"),
      ml_predicted_expected_value = if (is.na(candidates$pred_expected[i])) NULL else round(as.numeric(candidates$pred_expected[i]), 3),
      ml_predicted_modal_value = if (is.na(candidates$pred_class[i])) NULL else as.integer(candidates$pred_class[i]),
      ml_outlier_score = if (!is.finite(candidates$outlier_score[i])) NULL else round(as.numeric(candidates$outlier_score[i]), 3)
    )
  })

  payload <- jsonlite::toJSON(rows, auto_unbox = TRUE, pretty = TRUE)

  paste0(
    "Task: audit the current 0-4 annotation values for ONE question and a small set of fragrances\n",
    "(perfumes, colognes, eaux de toilette, extraits, and similar).\n\n",
    "Question/column name:\n", target_col_name, "\n\n",
    "Question to answer:\n", question_text, "\n\n",
    "Scale:\n",
    "0 = Definitely No\n",
    "1 = Probably No\n",
    "2 = Unknown / insufficient reliable evidence / release cannot be pinned down\n",
    "3 = Probably Yes\n",
    "4 = Definitely Yes\n\n",
    "Source guidance:\n", SOURCE_GUIDANCE, "\n\n",
    "Fragrance interpretation policy:\n", FRAGRANCE_INTERPRETATION_POLICY, "\n\n",
    "Items to audit. The ML fields explain why these cells were selected but must not determine your answer:\n",
    payload, "\n\n",
    "Return ONLY this JSON shape, preserving the same fragrance order:\n",
    "{\"results\":[{\"fragrance\":\"...\",\"row_index\":1,\"current_value\":0,\"recommended_value\":4,\"confidence\":0.0,\"needs_change\":true,\"rationale\":\"one short sentence\",\"sources\":[\"source or URL\"]}]}\n\n",
    "Rules:\n",
    "- Echo the fragrance name exactly as given, including any parenthetical house.\n",
    "- current_value may be null when the input cell is genuinely missing.\n",
    "- recommended_value must be an integer 0,1,2,3,or 4; use 2 if the correct researched answer remains genuinely unknown.\n",
    "- confidence must be a number from 0 to 1 reflecting source strength and how certainly you identified the specific release.\n",
    "- needs_change should be true exactly when recommended_value differs from current_value; for null current_value, any valid recommended_value counts as a change.\n",
    "- If sources conflict, choose 1 or 3 rather than 0 or 4 unless the stronger source is decisive.\n",
    "- A current_value of 2 often means the earlier pass never researched this item rather than that it researched it and failed. Do not treat an existing 2 as evidence.\n",
    "- If you cannot identify the specific release or answer it after searching, recommended_value must be 2.\n",
    "- Include at least one source string when confidence is above 0.5.\n"
  )
}

# ----------------------------- API HELPERS ------------------------------------

extract_response_text <- function(parsed) {
  # OpenAI/xAI Responses often use output -> message -> content -> output_text.
  text_parts <- character()

  if (!is.null(parsed$output) && is.list(parsed$output)) {
    for (item in parsed$output) {
      if (is.list(item) && identical(item$type %||% NULL, "message") && !is.null(item$content)) {
        for (ct in item$content) {
          if (is.list(ct) && ((ct$type %||% "") %in% c("output_text", "text")) && is.character(ct$text %||% NULL)) {
            text_parts <- c(text_parts, ct$text)
          }
        }
      }
    }
  }

  # Anthropic Messages uses content blocks with type="text".
  if (length(text_parts) == 0L && !is.null(parsed$content) && is.list(parsed$content)) {
    for (block in parsed$content) {
      if (is.list(block) && identical(block$type %||% NULL, "text") && is.character(block$text %||% NULL)) {
        text_parts <- c(text_parts, block$text)
      }
    }
  }

  # Some SDK-compatible servers expose output_text directly.
  if (length(text_parts) == 0L && is.character(parsed$output_text %||% NULL)) {
    text_parts <- c(text_parts, parsed$output_text)
  }

  paste(text_parts, collapse = "\n")
}

count_search_uses <- function(parsed) {
  n <- 0L
  if (!is.null(parsed$usage$server_tool_use$web_search_requests)) {
    n <- n + as.integer(parsed$usage$server_tool_use$web_search_requests)
  }
  if (!is.null(parsed$output) && is.list(parsed$output)) {
    n <- n + sum(vapply(parsed$output, function(item) {
      is.list(item) && grepl("web_search", as.character(item$type %||% ""), fixed = TRUE)
    }, logical(1)))
  }
  n
}

raw_response_path <- function(prefix) {
  safe_dir_create(RAW_RESPONSE_DIR)
  file.path(
    RAW_RESPONSE_DIR,
    sprintf("%s_%s_%06d.json", prefix, format(Sys.time(), "%Y%m%d_%H%M%S"), sample.int(999999, 1))
  )
}

# ---- Request-shaping helpers -------------------------------------------------

effort_value <- function() {
  e <- REASONING_EFFORT
  if (is.null(e) || length(e) == 0L) return(NA_character_)
  e <- as.character(e)[1]
  if (is.na(e) || !nzchar(trimws(e))) return(NA_character_)
  tolower(trimws(e))
}

apply_temperature <- function(body) {
  # Deliberately conservative. Anthropic's current models 400 on any non-default
  # temperature, and the Grok/GPT reasoning models either reject sampling
  # parameters or ignore them, so nothing is gained by sending it blind.
  if (!isTRUE(SEND_TEMPERATURE)) return(body)
  if (is.null(LLM_TEMPERATURE) || is.na(LLM_TEMPERATURE)) return(body)
  if (identical(PROVIDER, "anthropic")) {
    log_error("SEND_TEMPERATURE is TRUE but temperature is not accepted by current Claude models; omitting it.")
    return(body)
  }
  body$temperature <- LLM_TEMPERATURE
  body
}

# Parameters that a provider may start rejecting outright. When a 400 names one
# of these, we strip it and try once more rather than failing the whole cell.
DROPPABLE_PARAMS <- c("temperature", "tool_choice", "reasoning", "effort", "output_config", "top_p", "top_k")

params_named_in_error <- function(raw_body) {
  if (is.null(raw_body) || !nzchar(raw_body)) return(character())
  lowered <- tolower(raw_body)
  DROPPABLE_PARAMS[vapply(DROPPABLE_PARAMS, function(p) grepl(p, lowered, fixed = TRUE), logical(1))]
}

drop_params <- function(body, params) {
  for (p in params) {
    if (!is.null(body[[p]])) body[[p]] <- NULL
    # Anthropic nests effort inside output_config.
    if (identical(p, "effort") && !is.null(body$output_config)) {
      body$output_config$effort <- NULL
      if (length(body$output_config) == 0L) body$output_config <- NULL
    }
  }
  body
}

post_json_with_retries <- function(url, headers, body, provider_name, allow_param_fallback = TRUE) {
  for (attempt in seq_len(LLM_MAX_RETRIES)) {
    response <- tryCatch({
      httr::POST(
        url = url,
        httr::add_headers(.headers = headers),
        body = jsonlite::toJSON(body, auto_unbox = TRUE, null = "null"),
        encode = "json",
        httr::timeout(LLM_TIMEOUT_SECONDS)
      )
    }, error = function(e) e)

    if (inherits(response, "error")) {
      log_error(sprintf("%s request error on attempt %d: %s", provider_name, attempt, response$message))
      Sys.sleep(min(60, 5 * attempt))
      next
    }

    status <- httr::status_code(response)
    raw_body <- httr::content(response, "text", encoding = "UTF-8")

    if (status >= 200L && status < 300L) {
      parsed <- tryCatch(jsonlite::fromJSON(raw_body, simplifyVector = FALSE), error = function(e) NULL)
      if (is.null(parsed)) {
        log_error(sprintf("%s returned non-JSON API body; saved raw response", provider_name))
        writeLines(raw_body, raw_response_path(paste0(provider_name, "_api_nonjson")))
        return(list(text = NULL, parsed = NULL, raw = raw_body, status = status, search_uses = NA_integer_))
      }
      return(list(
        text = extract_response_text(parsed),
        parsed = parsed,
        raw = raw_body,
        status = status,
        search_uses = count_search_uses(parsed)
      ))
    }

    retryable <- status %in% c(408L, 409L, 425L, 429L, 500L, 502L, 503L, 504L, 529L)
    log_error(sprintf(
      "%s API status %d on attempt %d: %s",
      provider_name, status, attempt, substr(raw_body, 1, 500)
    ))

    # A 400 that names one of our optional request parameters almost always
    # means that parameter was deprecated or restricted on the chosen model
    # (this is how Claude 4.7+ rejects temperature, for example). Strip the
    # named parameters and try once more, so a provider change degrades the
    # request instead of killing the run. Only done once per call, and only
    # when something was actually removed.
    if (status == 400L && isTRUE(allow_param_fallback)) {
      offending <- intersect(params_named_in_error(raw_body), names(body))
      if (length(offending) == 0L && grepl("effort", tolower(raw_body), fixed = TRUE) &&
          !is.null(body$output_config$effort)) {
        offending <- "effort"
      }
      if (length(offending) > 0L) {
        log_error(sprintf(
          "%s 400 names %s; retrying once without %s.",
          provider_name, paste(offending, collapse = ", "), paste(offending, collapse = ", ")
        ))
        return(post_json_with_retries(
          url = url,
          headers = headers,
          body = drop_params(body, offending),
          provider_name = paste0(provider_name, "_no_", paste(offending, collapse = "_")),
          allow_param_fallback = FALSE
        ))
      }
    }

    if (!retryable) {
      return(list(text = NULL, parsed = NULL, raw = raw_body, status = status, search_uses = NA_integer_))
    }
    Sys.sleep(min(90, 10 * attempt))
  }

  log_error(sprintf("%s max retries exceeded", provider_name))
  list(text = NULL, parsed = NULL, raw = NULL, status = NA_integer_, search_uses = NA_integer_)
}

call_openai <- function(prompt) {
  api_key <- assert_api_key(OPENAI_API_KEY_ENV, "openai")

  tool <- list(type = "web_search", search_context_size = OPENAI_SEARCH_CONTEXT_SIZE, external_web_access = TRUE)
  body <- list(
    model = OPENAI_MODEL,
    input = list(
      list(role = "system", content = SYSTEM_PROMPT),
      list(role = "user", content = prompt)
    ),
    tools = list(tool)
  )

  # OpenAI Responses API: reasoning.effort. Supported values are
  # model-dependent and can include none/minimal/low/medium/high/xhigh/max;
  # gpt-5.6 defaults to medium. Omitted entirely when REASONING_EFFORT is blank.
  eff <- effort_value()
  if (!is.na(eff)) body$reasoning <- list(effort = eff)

  body <- apply_temperature(body)
  if (OPENAI_FORCE_WEB_SEARCH) body$tool_choice <- "required"

  post_json_with_retries(
    url = OPENAI_API_URL,
    headers = c("Authorization" = paste("Bearer", api_key), "Content-Type" = "application/json"),
    body = body,
    provider_name = "openai"
  )
}

call_xai <- function(prompt) {
  api_key <- assert_api_key(XAI_API_KEY_ENV, "xai")

  body <- list(
    model = XAI_MODEL,
    input = list(
      list(role = "system", content = SYSTEM_PROMPT),
      list(role = "user", content = prompt)
    ),
    tools = list(list(type = "web_search"))
  )

  # xAI Responses API: reasoning.effort -> {"reasoning":{"effort":"high"}}.
  # This is the field that was missing from the PharmaDive script, and it is
  # the same shape the FragranceDive batch script uses in make_request().
  # Do NOT send the Chat Completions form (a top-level "reasoning_effort"
  # string) to /v1/responses; it is ignored rather than rejected, which is the
  # worst possible failure mode. grok-4.6 accepts low/medium/high/xhigh and
  # defaults to high; reasoning cannot be disabled on 4.5+.
  eff <- effort_value()
  if (!is.na(eff)) body$reasoning <- list(effort = eff)

  body <- apply_temperature(body)
  if (XAI_FORCE_WEB_SEARCH) body$tool_choice <- "required"

  post_json_with_retries(
    url = XAI_API_URL,
    headers = c("Authorization" = paste("Bearer", api_key), "Content-Type" = "application/json"),
    body = body,
    provider_name = "xai"
  )
}

anthropic_web_search_tool <- function() {
  tool <- list(
    type = ANTHROPIC_WEB_SEARCH_TOOL_TYPE,
    name = "web_search",
    max_uses = ANTHROPIC_MAX_SEARCH_USES
  )
  # On web_search_20260209 and later, allowed_callers defaults to
  # ["code_execution_20260120"], i.e. search runs from inside code execution
  # (dynamic filtering) rather than directly. Set
  # ANTHROPIC_WEB_SEARCH_ALLOWED_CALLERS <- "direct" to pin the older
  # behaviour when moving to a newer tool version.
  if (!is.na(ANTHROPIC_WEB_SEARCH_ALLOWED_CALLERS) &&
      nzchar(ANTHROPIC_WEB_SEARCH_ALLOWED_CALLERS)) {
    tool$allowed_callers <- list(ANTHROPIC_WEB_SEARCH_ALLOWED_CALLERS)
  }
  tool
}

anthropic_body <- function(messages) {
  body <- list(
    model = ANTHROPIC_MODEL,
    max_tokens = ANTHROPIC_MAX_TOKENS,
    system = SYSTEM_PROMPT,
    messages = messages,
    tools = list(anthropic_web_search_tool())
  )

  # Anthropic puts effort in a top-level output_config object. It must not go
  # inside `thinking` (that is a 400), and no beta header is needed on current
  # models. The API default is high, so passing "high" is equivalent to
  # omitting it. Note that on Sonnet 5 / Opus 5 adaptive thinking is on by
  # default and max_tokens caps thinking + visible output together, so higher
  # effort makes stop_reason "max_tokens" more likely -- see the comment on
  # ANTHROPIC_MAX_TOKENS in the config block.
  eff <- effort_value()
  if (!is.na(eff)) body$output_config <- list(effort = eff)

  apply_temperature(body)
}

call_anthropic <- function(prompt) {
  api_key <- assert_api_key(ANTHROPIC_API_KEY_ENV, "anthropic")

  headers <- c(
    "x-api-key" = api_key,
    "anthropic-version" = ANTHROPIC_VERSION,
    "Content-Type" = "application/json"
  )

  body <- anthropic_body(list(list(role = "user", content = prompt)))

  result <- post_json_with_retries(ANTHROPIC_API_URL, headers, body, "anthropic")

  # Handle Anthropic pause_turn, mirroring the continuation pattern in the user's example scripts.
  parsed <- result$parsed
  if (!is.null(parsed) && identical(parsed$stop_reason %||% "", "pause_turn")) {
    all_content <- parsed$content
    for (k in seq_len(3L)) {
      cont_body <- anthropic_body(list(
        list(role = "user", content = prompt),
        list(role = "assistant", content = all_content),
        list(role = "user", content = "Continue. Return only the final JSON object.")
      ))
      cont <- post_json_with_retries(ANTHROPIC_API_URL, headers, cont_body, "anthropic_continuation")
      if (is.null(cont$parsed)) break
      all_content <- c(all_content, cont$parsed$content)
      result$text <- paste(result$text, extract_response_text(cont$parsed), sep = "\n")
      result$raw <- paste(result$raw, cont$raw, sep = "\n")
      result$search_uses <- sum(result$search_uses %||% 0L, cont$search_uses %||% 0L, na.rm = TRUE)
      if (!identical(cont$parsed$stop_reason %||% "", "pause_turn")) break
    }
  }

  # Truncation is silent as far as JSON parsing is concerned: you just get an
  # unterminated object and a parse_error. Say so explicitly instead.
  if (!is.null(result$parsed) && identical(result$parsed$stop_reason %||% "", "max_tokens")) {
    log_error(sprintf(
      "Anthropic stopped at max_tokens (%d). Raise ANTHROPIC_MAX_TOKENS or lower REASONING_EFFORT; the JSON is probably truncated.",
      ANTHROPIC_MAX_TOKENS
    ))
  }

  result
}

call_llm <- function(prompt) {
  if (PROVIDER == "openai") return(call_openai(prompt))
  if (PROVIDER == "anthropic") return(call_anthropic(prompt))
  if (PROVIDER == "xai") return(call_xai(prompt))
  stop("Unsupported PROVIDER: ", PROVIDER, call. = FALSE)
}

# ---------------------------- JSON PARSING ------------------------------------

strip_code_fences <- function(text) {
  text <- gsub("```json", "", text, fixed = TRUE)
  text <- gsub("```", "", text, fixed = TRUE)
  trimws(text)
}

extract_json_object_text <- function(text) {
  text <- strip_code_fences(text)
  if (!nzchar(text)) return(text)
  first <- regexpr("\\{", text)[1]
  last_matches <- gregexpr("\\}", text)[[1]]
  if (first < 0L || length(last_matches) == 0L || last_matches[1] < 0L) return(text)
  last <- tail(last_matches, 1)
  substr(text, first, last)
}

parse_llm_review <- function(response_text, candidates) {
  expected <- as.character(candidates$fragrance)
  expected_norm <- normalize_name(expected)

  if (is.null(response_text) || !nzchar(trimws(response_text))) {
    return(list(status = "empty_response", dt = data.table(), cleaned = ""))
  }

  cleaned <- extract_json_object_text(response_text)
  parsed <- tryCatch(jsonlite::fromJSON(cleaned, simplifyVector = FALSE), error = function(e) e)
  if (inherits(parsed, "error")) {
    return(list(status = paste0("parse_error: ", parsed$message), dt = data.table(), cleaned = cleaned))
  }

  results <- parsed$results %||% NULL
  if (is.null(results) && is.list(parsed) && all(vapply(parsed, is.list, logical(1)))) {
    results <- parsed
  }
  if (is.null(results) || !is.list(results)) {
    return(list(status = "parse_error: JSON did not contain a results array", dt = data.table(), cleaned = cleaned))
  }

  item_to_dt <- function(item) {
    sources <- item$sources %||% character()
    if (is.list(sources)) sources <- vapply(sources, function(z) paste(unlist(z), collapse = " "), character(1))
    data.table(
      fragrance = as.character(item$fragrance %||% NA_character_),
      row_index = suppressWarnings(as.integer(item$row_index %||% NA_integer_)),
      current_value_llm = suppressWarnings(as.integer(item$current_value %||% NA_integer_)),
      recommended_value = suppressWarnings(as.integer(item$recommended_value %||% NA_integer_)),
      confidence = suppressWarnings(as.numeric(item$confidence %||% NA_real_)),
      needs_change_llm = as.logical(item$needs_change %||% NA),
      rationale = as.character(item$rationale %||% NA_character_),
      sources = paste(as.character(sources), collapse = " | ")
    )
  }

  dt <- rbindlist(lapply(results, item_to_dt), fill = TRUE)
  if (nrow(dt) == 0L) return(list(status = "parse_error: empty results", dt = data.table(), cleaned = cleaned))

  dt[, fragrance_norm := normalize_name(fragrance)]
  dt[, fragrance_base := normalize_name_base(fragrance)]
  expected_base <- normalize_name_base(expected)

  # Align to candidates by row_index first, then fragrance name. This prevents one malformed name
  # from shifting all later rows.
  aligned <- vector("list", nrow(candidates))
  for (i in seq_len(nrow(candidates))) {
    ri <- as.integer(candidates$row_index[i])
    hit <- dt[row_index == ri]
    if (nrow(hit) == 0L) hit <- dt[fragrance_norm == expected_norm[i]]
    if (nrow(hit) == 0L) {
      # Last resort: the model dropped or reworded the house parenthetical
      # ("Aventus" for "Aventus (Creed)"). Accept only when the stripped name
      # is unique among both the fragrances we asked about and the rows we got
      # back, so a batch containing two same-named fragrances from different
      # houses can never be cross-assigned.
      if (sum(expected_base == expected_base[i]) == 1L) {
        base_hits <- dt[fragrance_base == expected_base[i]]
        if (nrow(base_hits) == 1L) hit <- base_hits
      }
    }
    if (nrow(hit) == 0L) {
      aligned[[i]] <- data.table(
        fragrance = expected[i],
        row_index = ri,
        current_value_llm = NA_integer_,
        recommended_value = NA_integer_,
        confidence = NA_real_,
        needs_change_llm = NA,
        rationale = "No matching LLM result for this fragrance/row_index.",
        sources = NA_character_
      )
    } else {
      aligned[[i]] <- hit[1, .(
        fragrance,
        row_index,
        current_value_llm,
        recommended_value,
        confidence,
        needs_change_llm,
        rationale,
        sources
      )]
    }
  }

  out <- rbindlist(aligned, fill = TRUE)
  bad_value <- is.na(out$recommended_value) | !(out$recommended_value %in% 0:4)
  if (any(bad_value)) {
    out[bad_value, recommended_value := NA_integer_]
  }
  out[is.na(confidence) | confidence < 0 | confidence > 1, confidence := NA_real_]
  list(status = "parsed", dt = out, cleaned = cleaned)
}

should_apply_change <- function(current_value, recommended_value, confidence, candidate_type = "outlier") {
  if (!APPLY_LLM_CORRECTIONS || DRY_RUN) return(FALSE)
  if (is.na(recommended_value) || !(recommended_value %in% 0:4)) return(FALSE)

  current_is_missing <- is.na(current_value) || identical(candidate_type, "missing")
  if (!current_is_missing && identical(as.integer(current_value), as.integer(recommended_value))) return(FALSE)

  if (current_is_missing && isTRUE(APPLY_ANY_VALID_RECOMMENDATION_TO_MISSING)) return(TRUE)

  threshold <- MIN_CONFIDENCE_TO_APPLY
  if (current_is_missing && !is.na(MIN_CONFIDENCE_TO_APPLY_MISSING)) {
    threshold <- MIN_CONFIDENCE_TO_APPLY_MISSING
  }
  if (is.na(confidence) || confidence < threshold) return(FALSE)
  TRUE
}

manual_approve <- function(fragrance, question, current_value, recommended_value, confidence, rationale, sources) {
  if (!MANUAL_APPROVAL) return(TRUE)
  cat("\n---------------- Manual approval ----------------\n")
  cat("Fragrance: ", fragrance, "\n", sep = "")
  cat("Question: ", question, "\n", sep = "")
  cat(sprintf("Current: %s | Recommended: %s | Confidence: %.2f\n", current_value, recommended_value, confidence))
  cat("Rationale: ", rationale, "\n", sep = "")
  cat("Sources: ", sources, "\n", sep = "")
  ans <- readline("Apply change? [y/N]: ")
  tolower(trimws(ans)) %in% c("y", "yes")
}

current_llm_model <- function() {
  switch(PROVIDER, openai = OPENAI_MODEL, anthropic = ANTHROPIC_MODEL, xai = XAI_MODEL, PROVIDER)
}

make_batch_failure_log <- function(batch, sweep, target_local_idx, csv_col_idx, target_col_name,
                                   question_text, parse_status, api_result = list(), raw_path = NA_character_) {
  batch <- as.data.table(batch)
  t_local <- as.integer(target_local_idx)
  csv_col <- as.integer(csv_col_idx)
  status_text <- as.character(parse_status)
  raw_path_text <- raw_path %||% NA_character_

  batch[, .(
    timestamp = format(Sys.time(), "%Y-%m-%d %H:%M:%S"),
    provider = PROVIDER,
    model = current_llm_model(),
    sweep = as.integer(sweep),
    target_local_idx = t_local,
    target_csv_col_idx = csv_col,
    target_col_name = target_col_name,
    question_text = question_text,
    row_index = row_index,
    fragrance = fragrance,
    candidate_type = candidate_type,
    original_value_at_review = target_value,
    recommended_value = NA_integer_,
    final_value = target_value,
    confidence = NA_real_,
    change_made = FALSE,
    llm_parse_status = status_text,
    needs_change_llm = NA,
    rationale = NA_character_,
    sources = NA_character_,
    search_uses = api_result$search_uses %||% NA_integer_,
    api_status = api_result$status %||% NA_integer_,
    raw_response_file = raw_path_text,
    pred_expected = pred_expected,
    pred_class = pred_class,
    true_prob = true_prob,
    outlier_score = outlier_score,
    model_oob_error = model_oob_error
  )]
}

process_llm_batch <- function(df, batch, question_text, target_col_name, target_local_idx,
                              qcols, sweep, batch_number, batch_count) {
  # Treat each LLM request as a small transaction: build log rows and edits on a
  # working copy, then return them together. The caller keeps the pre-batch df if
  # this function throws, so one bad query cannot leave partial in-memory edits.
  df_work <- df
  csv_col_idx <- qcols[target_local_idx]
  batch_review_rows <- list()

  prompt <- create_review_prompt(question_text, target_col_name, batch)
  api_result <- tryCatch(call_llm(prompt), error = function(e) {
    log_error(sprintf("LLM call failed for %s: %s", target_col_name, e$message))
    list(text = NULL, raw = NULL, status = NA_integer_, search_uses = NA_integer_)
  })

  parsed <- parse_llm_review(api_result$text, batch)
  if (!identical(parsed$status, "parsed")) {
    log_error(sprintf("LLM parse failed for %s: %s", target_col_name, parsed$status))
    if (PRINT_RAW_LLM_ON_PARSE_ERROR && !is.null(api_result$raw)) {
      p <- raw_response_path("parse_error")
      writeLines(api_result$raw, p)
      log_error(sprintf("Raw LLM response saved to %s", p))
    }
  }

  # Save raw API body for auditability. This remains per-query because the raw
  # responses are separate audit artifacts, not the main checkpoint/output/log files.
  raw_path <- NA_character_
  if (!is.null(api_result$raw) && nzchar(api_result$raw)) {
    raw_path <- raw_response_path("llm")
    writeLines(api_result$raw, raw_path)
  }

  if (identical(parsed$status, "parsed")) {
    for (i in seq_len(nrow(batch))) {
      row_idx <- as.integer(batch$row_index[i])
      current_value <- as.integer(df_work[[csv_col_idx]][row_idx])

      rec <- parsed$dt[i]
      recommended_value <- as.integer(rec$recommended_value)
      confidence <- as.numeric(rec$confidence)
      candidate_type <- as.character(batch$candidate_type[i] %||% "outlier")
      apply_now <- should_apply_change(current_value, recommended_value, confidence, candidate_type)

      if (apply_now) {
        apply_now <- manual_approve(
          fragrance = as.character(df_work[[ID_COL]][row_idx]),
          question = question_text,
          current_value = current_value,
          recommended_value = recommended_value,
          confidence = confidence,
          rationale = rec$rationale,
          sources = rec$sources
        )
      }

      final_value <- current_value
      change_made <- FALSE
      if (apply_now) {
        final_value <- recommended_value
        df_work[[csv_col_idx]][row_idx] <- final_value
        change_made <- TRUE
      }

      log_row <- data.table(
        timestamp = format(Sys.time(), "%Y-%m-%d %H:%M:%S"),
        provider = PROVIDER,
        model = current_llm_model(),
        sweep = sweep,
        target_local_idx = target_local_idx,
        target_csv_col_idx = csv_col_idx,
        target_col_name = target_col_name,
        question_text = question_text,
        row_index = row_idx,
        fragrance = as.character(df_work[[ID_COL]][row_idx]),
        candidate_type = candidate_type,
        original_value_at_review = current_value,
        recommended_value = recommended_value,
        final_value = final_value,
        confidence = confidence,
        change_made = change_made,
        llm_parse_status = parsed$status,
        needs_change_llm = rec$needs_change_llm,
        rationale = rec$rationale,
        sources = rec$sources,
        search_uses = api_result$search_uses %||% NA_integer_,
        api_status = api_result$status %||% NA_integer_,
        raw_response_file = raw_path,
        pred_expected = as.numeric(batch$pred_expected[i]),
        pred_class = as.integer(batch$pred_class[i]),
        true_prob = as.numeric(batch$true_prob[i]),
        outlier_score = as.numeric(batch$outlier_score[i]),
        model_oob_error = as.numeric(batch$model_oob_error[i])
      )
      batch_review_rows[[length(batch_review_rows) + 1L]] <- log_row

      if (change_made) {
        log_info(
          "Changed row %d (%s), %s: %s -> %d [confidence %.2f; %s]",
          row_idx, as.character(df_work[[ID_COL]][row_idx]), target_col_name,
          format_value_for_log(current_value), final_value, confidence, candidate_type
        )
      } else {
        log_info(
          "No applied change row %d (%s), %s: current %s, recommended %s [confidence %s; %s]",
          row_idx, as.character(df_work[[ID_COL]][row_idx]), target_col_name,
          format_value_for_log(current_value),
          ifelse(is.na(recommended_value), "NA", as.character(recommended_value)),
          ifelse(is.na(confidence), "NA", sprintf("%.2f", confidence)),
          candidate_type
        )
      }
    }
  } else {
    # Log parse/API failures without marking them as successfully reviewed, so they can be retried.
    batch_review_rows[[1L]] <- make_batch_failure_log(
      batch = batch,
      sweep = sweep,
      target_local_idx = target_local_idx,
      csv_col_idx = csv_col_idx,
      target_col_name = target_col_name,
      question_text = question_text,
      parse_status = parsed$status,
      api_result = api_result,
      raw_path = raw_path
    )
  }

  list(
    df = df_work,
    review_log = combine_review_logs(data.table(), batch_review_rows)
  )
}

# ------------------------------- MAIN -----------------------------------------

main <- function() {
  set.seed(SEED)
  safe_dir_create(RAW_RESPONSE_DIR)
  validate_persistence_config()
  validate_request_config()

  df <- load_input_or_checkpoint()
  loaded_from_checkpoint <- isTRUE(attr(df, "loaded_from_checkpoint"))
  checkpoint_sweep <- suppressWarnings(as.integer(attr(df, "checkpoint_sweep") %||% 0L))
  checkpoint_event <- as.character(attr(df, "checkpoint_persistence_event") %||% NA_character_)
  if (nrow(df) == 0L) stop("Input has zero rows", call. = FALSE)
  if (ID_COL < 1L || ID_COL > ncol(df)) stop("Invalid ID_COL", call. = FALSE)

  qcols <- resolve_question_columns(df)
  target_local_indices <- resolve_target_local_indices(df, qcols)
  q_names <- names(df)[qcols]
  dup_q_names <- q_names[duplicated(q_names)]
  if (length(dup_q_names) > 0L) {
    log_info(
      "Duplicate question headers detected: %d repeated header instance(s) across %d distinct question text(s). They will be preserved in output and treated as separate columns by position.",
      sum(duplicated(q_names)), length(unique(dup_q_names))
    )
  }
  questions <- load_questions(df, qcols)

  # Generic raw-run behavior: literal 2s can be reclassified only on the original
  # input load. Reported statistical-analysis inputs convert all Unknown outcomes
  # to blanks, so no original code-2 observation is analyzed as a category.
  treat_2_this_load <- isTRUE(TREAT_2_AS_MISSING) && !loaded_from_checkpoint
  df <- coerce_question_columns(df, qcols, treat_2_as_missing = treat_2_this_load)

  # Same reasoning applies to the problem-log pass: it reclassifies 2s that the
  # batch script never researched, and that judgement only holds for the raw
  # batch output. Once we are resuming from a checkpoint, a 2 in those rows may
  # be an answer this script already verified, so leave it alone.
  if (!loaded_from_checkpoint) {
    df <- mark_unanswered_twos_as_missing(df, qcols)
  }

  n_missing_initial <- sum(is.na(as.matrix(df[, qcols, drop = FALSE])))
  log_info("Rows: %d | Question columns: %d | Target columns this run: %d | Provider: %s | Model type: %s", nrow(df), length(qcols), length(target_local_indices), PROVIDER, MODEL_TYPE)
  log_info("Missing-priority mode: %s | Initial genuinely missing/invalid cells: %d", PRIORITIZE_MISSING_VALUES, n_missing_initial)
  log_info("Treat literal 2s as missing on this load: %s | RF missing handling: %s", treat_2_this_load, model_missing_handling_description())
  log_info("Model: %s | Reasoning effort: %s | temperature sent: %s",
           current_llm_model(),
           ifelse(is.na(effort_value()), "(model default)", effort_value()),
           isTRUE(SEND_TEMPERATURE))
  log_info("Output: %s | Checkpoint: %s | Review log: %s | Persistence cadence: %s", OUTPUT_FILE, CHECKPOINT_FILE, REVIEW_LOG_FILE, SAVE_PROGRESS_EVERY)

  pending_review_rows <- list()
  has_persisted <- loaded_from_checkpoint
  last_target_local_idx <- NA_integer_

  start_sweep <- 1L
  if (loaded_from_checkpoint && !is.na(checkpoint_sweep) && checkpoint_sweep >= 1L) {
    # With the recommended sweep-level cadence, a saved sweep is complete, so
    # resume at the following sweep. For a non-sweep checkpoint, repeat the same
    # sweep rather than silently skipping unfinished work.
    if (identical(checkpoint_event, "sweep") || identical(checkpoint_event, "final")) {
      start_sweep <- checkpoint_sweep + 1L
    } else {
      start_sweep <- checkpoint_sweep
    }
    log_info(
      "Checkpoint records sweep %d after event '%s'; resuming at sweep %d.",
      checkpoint_sweep,
      ifelse(is.na(checkpoint_event), "unknown", checkpoint_event),
      start_sweep
    )
  }

  sweeps_to_run <- if (start_sweep <= N_SWEEPS) seq.int(start_sweep, N_SWEEPS) else integer()
  for (sweep in sweeps_to_run) {
    log_info("\n===== Sweep %d of %d =====", sweep, N_SWEEPS)
    n_cols_done <- 0L

    for (target_pos in seq_along(target_local_indices)) {
      if (n_cols_done >= MAX_COLUMNS_PER_SWEEP) break
      target_local_idx <- target_local_indices[target_pos]
      last_target_local_idx <- target_local_idx
      target_col_name <- q_names[target_local_idx]
      question_text <- questions[target_local_idx]
      log_info("\n--- Sweep %d | Target %d/%d (question column %d/%d): %s ---", sweep, target_pos, length(target_local_indices), target_local_idx, length(qcols), target_col_name)

      # Always refresh from df because previous columns may have been corrected.
      mat <- as.matrix(df[, qcols, drop = FALSE])
      storage.mode(mat) <- "integer"

      # Include both already-flushed review rows and rows buffered during the
      # current not-yet-persisted interval, so skip logic remains correct even
      # when REVIEW_LOG_FILE is only written once per sweep.
      review_log <- combine_review_logs(read_existing_review_log(), pending_review_rows)
      candidates <- tryCatch(
        select_outliers_for_column(df, mat, qcols, q_names, target_local_idx, review_log),
        error = function(e) {
          log_error(sprintf("Outlier scoring failed for %s: %s", target_col_name, e$message))
          data.table()
        }
      )

      if (nrow(candidates) == 0L) {
        log_info("No unchecked candidates for this column.")
        n_cols_done <- n_cols_done + 1L
        if (should_persist_after("column")) {
          persist_progress(
            df, sweep, target_local_idx, pending_review_rows,
            extra = list(persistence_event = "column", columns_completed = n_cols_done)
          )
          pending_review_rows <- list()
          has_persisted <- TRUE
        }
        if (SLEEP_BETWEEN_COLUMNS > 0) Sys.sleep(SLEEP_BETWEEN_COLUMNS)
        next
      }

      candidates[, sweep := sweep]
      n_missing_selected <- sum(candidates$candidate_type == "missing", na.rm = TRUE)
      n_outlier_selected <- sum(candidates$candidate_type == "outlier", na.rm = TRUE)
      finite_scores <- candidates$outlier_score[is.finite(candidates$outlier_score)]
      if (length(finite_scores) > 0L) {
        log_info(
          "Selected %d candidate(s): %d missing/invalid, %d outlier(s). Top outlier score: %.3f",
          nrow(candidates), n_missing_selected, n_outlier_selected, max(finite_scores, na.rm = TRUE)
        )
      } else {
        log_info(
          "Selected %d candidate(s): %d missing/invalid, %d outlier(s). Outlier model skipped.",
          nrow(candidates), n_missing_selected, n_outlier_selected
        )
      }

      batch_ids <- ceiling(seq_len(nrow(candidates)) / LLM_BATCH_SIZE)
      for (b in unique(batch_ids)) {
        batch <- candidates[batch_ids == b]
        log_info("LLM batch %d/%d for %s: %s", b, max(batch_ids), target_col_name, paste(batch$fragrance, collapse = "; "))

        batch_result <- tryCatch(
          process_llm_batch(
            df = df,
            batch = batch,
            question_text = question_text,
            target_col_name = target_col_name,
            target_local_idx = target_local_idx,
            qcols = qcols,
            sweep = sweep,
            batch_number = b,
            batch_count = max(batch_ids)
          ),
          error = function(e) {
            # Roll back this LLM request by keeping df unchanged and recording
            # a failure row for every candidate in the failed batch.
            log_error(sprintf(
              "LLM batch %d/%d failed for %s; retaining pre-batch data and continuing sweep: %s",
              b, max(batch_ids), target_col_name, e$message
            ))
            list(
              df = df,
              review_log = make_batch_failure_log(
                batch = batch,
                sweep = sweep,
                target_local_idx = target_local_idx,
                csv_col_idx = qcols[target_local_idx],
                target_col_name = target_col_name,
                question_text = question_text,
                parse_status = paste0("batch_error: ", e$message),
                api_result = list(status = NA_integer_, search_uses = NA_integer_),
                raw_path = NA_character_
              )
            )
          }
        )

        df <- batch_result$df
        if (!is.null(batch_result$review_log) && nrow(batch_result$review_log) > 0L) {
          pending_review_rows[[length(pending_review_rows) + 1L]] <- batch_result$review_log
        }

        if (should_persist_after("batch")) {
          persist_progress(
            df, sweep, target_local_idx, pending_review_rows,
            extra = list(persistence_event = "batch", batch = b, batches_in_column = max(batch_ids))
          )
          pending_review_rows <- list()
          has_persisted <- TRUE
        }

        if (SLEEP_BETWEEN_API_CALLS > 0) Sys.sleep(SLEEP_BETWEEN_API_CALLS)
      }

      n_cols_done <- n_cols_done + 1L
      if (should_persist_after("column")) {
        persist_progress(
          df, sweep, target_local_idx, pending_review_rows,
          extra = list(persistence_event = "column", columns_completed = n_cols_done)
        )
        pending_review_rows <- list()
        has_persisted <- TRUE
      }
      if (SLEEP_BETWEEN_COLUMNS > 0) Sys.sleep(SLEEP_BETWEEN_COLUMNS)
    }

    if (should_persist_after("sweep")) {
      persist_progress(
        df, sweep, last_target_local_idx, pending_review_rows,
        extra = list(persistence_event = "sweep", columns_completed = n_cols_done)
      )
      pending_review_rows <- list()
      has_persisted <- TRUE
    }
  }

  if (length(pending_review_rows) > 0L || !has_persisted) {
    persist_progress(
      df, N_SWEEPS, last_target_local_idx, pending_review_rows,
      extra = list(persistence_event = "final")
    )
    pending_review_rows <- list()
    has_persisted <- TRUE
  }

  log_info("\nDone. Corrected CSV saved to: %s", OUTPUT_FILE)
  log_info("Review log saved to: %s", REVIEW_LOG_FILE)
  log_info("Checkpoint saved to: %s", CHECKPOINT_FILE)

  # Final summary
  mat_final <- as.matrix(df[, qcols, drop = FALSE])
  storage.mode(mat_final) <- "integer"
  counts <- table(factor(as.vector(mat_final), levels = 0:4))
  n_missing_final <- sum(is.na(mat_final))
  log_info("Final response distribution: %s | missing=%d", paste(sprintf("%s=%s", names(counts), as.integer(counts)), collapse = ", "), n_missing_final)
}

main()

#########################################################################################################
#########################################################################################################
# This concludes the correction script.
# There is an additional phase not shown here because it is essentially identical to above: the factor 
# score plotting of the 50 well-known fragrances is repeated using factor scores calculated from the 
# corrected data.
# The next phase shown here generates confusion matrices for repeated items in the full data set, the 
# 200 random fragrances, and the 200 well-known fragrances.
#########################################################################################################
#########################################################################################################

# ============================================================================
# REPEATED FRAGRANCE ITEMS: ORIGINAL VERSUS CORRECTED DATA
# Creates a six-panel main-text figure, complete Supplement figures, and tables.
#
# RUN: Put this script and both CSV files in the same directory, then:
#      source("fragrance_repeat_agreement.R")
# Or:  Rscript fragrance_repeat_agreement.R original.csv corrected.csv output_dir
#
# Only two add-on packages are needed: polycor and irr. Graphics use base R/grid.
# The input CSV files are never modified.
# ============================================================================

# ----------------------------- USER SETTINGS --------------------------------
script_directory <- function() {
  # source() in R/RStudio and Rscript are both supported.
  for (fr in rev(sys.frames())) {
    if (!is.null(fr$ofile)) {
      return(dirname(normalizePath(fr$ofile, mustWork = TRUE)))
    }
  }
  a <- grep("^--file=", commandArgs(trailingOnly = FALSE), value = TRUE)
  if (length(a)) return(dirname(normalizePath(sub("^--file=", "", a[1]),
                                             mustWork = TRUE)))
  getwd()
}
INPUT_DIR <- script_directory()
ORIGINAL_FILE <- file.path(INPUT_DIR, "fragrance_data_raw_200.csv")
CORRECTED_FILE <- file.path(INPUT_DIR, "fragrance_data_raw_200_corrected.csv")
OUTPUT_DIR <- file.path(INPUT_DIR, "fragrance_repeat_agreement_output")

AUTO_INSTALL_PACKAGES <- TRUE
ID_COLUMN <- 1L                    # Fragrance name; excluded from duplicate search.
RESPONSE_LEVELS <- c(0, 1, 3, 4)   # The gap between 1 and 3 is retained for ICC.
PAIRING <- "all_pairs"            # Or "first_vs_repeats".
POLYCHORIC_METHOD <- "two_step"   # Or "ML"; two-step is stable for sparse tables.
POLYCHORIC_MAX_ABS <- 0.9999
CONFIDENCE_LEVEL <- 0.95           # ICC confidence limits in CSV, not on figures.

# Main comparisons hold the fragrance sample fixed WITHIN each item pair.
# The requested Supplement uses every available pair within each file.
MAIN_SAMPLE <- "available"          # "matched" or "available".
WRITE_MATCHED_SUPPLEMENT <- TRUE   # Additional sensitivity figure, same sample as main.
MAIN_PAIR_IDS <- c("P14", "P11", "P20")              # E.g., c("P26", "P23", "P01"); see pair_map.csv.
MIN_MAIN_N <- 100L                 # Automatic selection prefers adequate sample sizes.

MAIN_WIDTH_IN <- 8.5
MAIN_HEIGHT_IN <- 11
SUPP_WIDTH_IN <- 11.7              # Approximately A4 landscape.
SUPP_HEIGHT_IN <- 8.3
SUPP_PAIR_COLUMNS <- 3L            # Three side-by-side pairs = SIX matrix columns.
SUPP_ROWS <- 3L                    # 18 matrices per landscape page.
PNG_DPI <- 400
TIFF_DPI <- 600
WRITE_TIFF <- TRUE
WRITE_SUPPLEMENT_PNGS <- TRUE
FONT_FAMILY <- "sans"

# Optional positional command-line arguments override the file paths.
args <- commandArgs(trailingOnly = TRUE)
if (length(args) >= 2L) {
  ORIGINAL_FILE <- args[1]
  CORRECTED_FILE <- args[2]
}
if (length(args) >= 3L) OUTPUT_DIR <- args[3]
if (!file.exists(CORRECTED_FILE) && length(args) < 2L) {
  alt <- file.path(INPUT_DIR, "fragrance_data_raw_200_corrected.csv")
  if (file.exists(alt)) CORRECTED_FILE <- alt
}

# ----------------------------- SETUP / INPUT --------------------------------
stopifnot(PAIRING %in% c("all_pairs", "first_vs_repeats"),
          POLYCHORIC_METHOD %in% c("two_step", "ML"),
          MAIN_SAMPLE %in% c("matched", "available"),
          identical(as.numeric(RESPONSE_LEVELS), c(0, 1, 3, 4)),
          SUPP_PAIR_COLUMNS >= 1L, SUPP_ROWS >= 1L,
          MAIN_WIDTH_IN >= 7.5, MAIN_HEIGHT_IN >= 10)
if (!file.exists(ORIGINAL_FILE) || !file.exists(CORRECTED_FILE)) {
  stop("Cannot locate both input files. Edit ORIGINAL_FILE and CORRECTED_FILE.\n",
       "Original: ", ORIGINAL_FILE, "\nCorrected: ", CORRECTED_FILE)
}
needed <- c("polycor", "irr")
missing_packages <- needed[!vapply(needed, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_packages)) {
  if (!AUTO_INSTALL_PACKAGES) {
    stop("Please install these packages: ", paste(missing_packages, collapse = ", "))
  }
  install.packages(missing_packages, repos = "https://cloud.r-project.org")
}
if (!all(vapply(needed, requireNamespace, logical(1), quietly = TRUE))) {
  stop("Package installation failed. Install polycor and irr, then rerun.")
}
for (d in c("figures", "tables", "diagnostics")) {
  dir.create(file.path(OUTPUT_DIR, d), recursive = TRUE, showWarnings = FALSE)
}
write_csv <- function(x, name, subdir = "tables") {
  utils::write.csv(x, file.path(OUTPUT_DIR, subdir, name), row.names = FALSE,
                   na = "NA", fileEncoding = "UTF-8")
}

read_ratings <- function(path) {
  d <- utils::read.csv(path, header = TRUE, check.names = FALSE,
                       colClasses = "character", stringsAsFactors = FALSE,
                       na.strings = c("", "NA", "NaN", "NULL","2"),
                       strip.white = TRUE, fileEncoding = "UTF-8-BOM")
  if (nrow(d) < 3L || ncol(d) < 3L) stop("Input is too small: ", path)
  if (ID_COLUMN < 1L || ID_COLUMN > ncol(d)) stop("Invalid ID_COLUMN.")
  id <- trimws(d[[ID_COLUMN]])
  if (anyNA(id) || any(!nzchar(id)) || anyDuplicated(id)) {
    stop("Fragrance identifiers must be nonmissing and unique: ", path)
  }
  d[[ID_COLUMN]] <- id
  for (j in setdiff(seq_len(ncol(d)), ID_COLUMN)) {
    raw <- trimws(d[[j]])
    raw[!is.na(raw) & raw == ""] <- NA_character_
    z <- suppressWarnings(as.numeric(raw))
    bad <- !is.na(raw) & (is.na(z) | !is.finite(z) | !(z %in% RESPONSE_LEVELS))
    if (any(bad)) {
      stop("Unexpected response in ", basename(path), ", column ", j, " (",
           names(d)[j], "): ", paste(unique(raw[bad]), collapse = ", "),
           ". Only 0, 1, 3, 4, and missing values are permitted.")
    }
    d[[j]] <- z
  }
  d
}
original <- read_ratings(ORIGINAL_FILE)
corrected <- read_ratings(CORRECTED_FILE)
if (!identical(names(original), names(corrected))) {
  stop("The CSV headers differ. Use matching columns, in the same order, in both files.")
}
if (!setequal(original[[ID_COLUMN]], corrected[[ID_COLUMN]])) {
  stop("The two files do not contain exactly the same fragrance IDs.\n",
       "Resolve this explicitly before comparing the correction.")
}
was_reordered <- !identical(original[[ID_COLUMN]], corrected[[ID_COLUMN]])
corrected <- corrected[match(original[[ID_COLUMN]], corrected[[ID_COLUMN]]), , drop = FALSE]
stopifnot(identical(original[[ID_COLUMN]], corrected[[ID_COLUMN]]))

# Detect both literal repeated headers and .1/.2/... suffixes. Work by column
# POSITION, so literally identical headers are not lost or silently renamed.
make_pair_map <- function(headers) {
  idx <- setdiff(seq_along(headers), ID_COLUMN)
  base <- sub("\\.[0-9]+$", "", trimws(headers[idx]))
  names(base) <- idx
  out <- list()
  for (item in unique(base)) {
    ii <- idx[base == item]
    if (length(ii) < 2L) next
    # File order defines occurrence 1, occurrence 2, etc.
    combinations <- if (PAIRING == "all_pairs") {
      t(utils::combn(seq_along(ii), 2L))
    } else {
      cbind(1L, seq.int(2L, length(ii)))
    }
    for (k in seq_len(nrow(combinations))) {
      a <- combinations[k, 1]; b <- combinations[k, 2]
      out[[length(out) + 1L]] <- data.frame(
        item = item, column1_index = ii[a], column2_index = ii[b],
        column1 = headers[ii[a]], column2 = headers[ii[b]],
        occurrence1 = a, occurrence2 = b, occurrences_in_item = length(ii),
        stringsAsFactors = FALSE, check.names = FALSE)
    }
  }
  if (!length(out)) stop("No repeated items detected. Inspect the column headers.")
  ans <- do.call(rbind, out)
  ans$pair_id <- sprintf("P%02d", seq_len(nrow(ans)))
  ans[, c("pair_id", setdiff(names(ans), "pair_id")), drop = FALSE]
}
pairs <- make_pair_map(names(original))
write_csv(pairs, "pair_map.csv")
message("Found ", length(unique(pairs$item)), " repeated items and ", nrow(pairs),
        " pairs (", 2L * nrow(pairs), " matrices per full Supplement).")

numeric_cols <- setdiff(seq_len(ncol(original)), ID_COLUMN)
a <- as.matrix(original[, numeric_cols, drop = FALSE])
b <- as.matrix(corrected[, numeric_cols, drop = FALSE])
quality <- data.frame(
  metric = c("fragrances", "rating_columns", "repeated_items", "item_pairs",
             "original_missing_cells", "corrected_missing_cells",
             "missing_to_observed_cells", "observed_to_missing_cells",
             "observed_values_changed", "corrected_rows_reordered_to_match"),
  value = c(nrow(original), length(numeric_cols), length(unique(pairs$item)), nrow(pairs),
            sum(is.na(a)), sum(is.na(b)), sum(is.na(a) & !is.na(b)),
            sum(!is.na(a) & is.na(b)), sum(a != b, na.rm = TRUE), was_reordered))
write_csv(quality, "input_quality.csv", "diagnostics")

# ------------------------------ STATISTICS ----------------------------------
# Capture warnings in the diagnostics table rather than concealing them.
capture_conditions <- function(expr) {
  messages <- character()
  error <- NULL
  value <- tryCatch(withCallingHandlers(expr, warning = function(w) {
    messages <<- c(messages, conditionMessage(w)); invokeRestart("muffleWarning")
  }), error = function(e) { error <<- conditionMessage(e); NULL })
  list(value = value, warnings = unique(messages), error = error)
}

polychoric_stat <- function(tab) {
  used <- tab[rowSums(tab) > 0, colSums(tab) > 0, drop = FALSE]
  sparse <- sum(tab == 0L)
  if (sum(tab) < 3L || nrow(used) < 2L || ncol(used) < 2L) {
    return(list(r = NA_real_, boundary = FALSE, method = POLYCHORIC_METHOD,
                note = "Not estimable: too few pairs or a constant occurrence.",
                occupied_rows = nrow(used), occupied_columns = ncol(used),
                zero_cells = sparse))
  }
  # Removing empty margins is for estimation only. The plotted table always
  # retains all four categories, and no pseudocount/continuity correction is used.
  fit <- capture_conditions(polycor::polychor(
    used, ML = identical(POLYCHORIC_METHOD, "ML"), std.err = FALSE,
    maxcor = POLYCHORIC_MAX_ABS, control = list(maxit = 2000L)))
  rho <- if (is.null(fit$value)) NA_real_ else as.numeric(fit$value)
  if (length(rho) != 1L || !is.finite(rho) || abs(rho) > 1) rho <- NA_real_
  boundary <- is.finite(rho) && abs(rho) >= 0.999
  sparse_note <- if (min(c(rowSums(used), colSums(used))) < 5L) {
    "Sparse margins: an occupied category has fewer than 5 observations; r may be unstable."
  } else NULL
  notes <- c(fit$error, fit$warnings, sparse_note,
             if (boundary) "Near correlation boundary; interpret cautiously." else NULL)
  if (!is.finite(rho) && !length(notes)) notes <- "Non-finite estimate."
  list(r = rho, boundary = boundary, method = POLYCHORIC_METHOD,
       note = paste(notes, collapse = " | "), occupied_rows = nrow(used),
       occupied_columns = ncol(used), zero_cells = sparse)
}

icc_stat <- function(x, y) {
  # ICC(A,1): two-way, absolute agreement, SINGLE measurement. NOT ICC(C,1),
  # NOT ICC(3,1) consistency, and NOT an average-measures ICC.
  # The numeric scores remain 0,1,3,4. Do not replace them with factor integers.
  n <- length(x); z <- cbind(x, y)
  if (n < 3L) return(list(icc = NA_real_, lower = NA_real_, upper = NA_real_,
                         note = "Too few complete pairs."))
  grand <- mean(z)
  ms_subject <- 2 * sum((rowMeans(z) - grand)^2) / (n - 1)
  ms_occurrence <- n * sum((colMeans(z) - grand)^2)
  residual <- sweep(sweep(z, 1L, rowMeans(z), "-"), 2L, colMeans(z), "-") + grand
  ms_error <- sum(residual^2) / (n - 1)
  denominator <- ms_subject + ms_error + 2 * (ms_occurrence - ms_error) / n
  if (!is.finite(denominator) || abs(denominator) <= .Machine$double.eps * 100) {
    return(list(icc = NA_real_, lower = NA_real_, upper = NA_real_,
                note = "Not estimable: zero variance/zero ICC denominator."))
  }
  formula_value <- (ms_subject - ms_error) / denominator
  fit <- capture_conditions(irr::icc(z, model = "twoway", type = "agreement",
                                     unit = "single", conf.level = CONFIDENCE_LEVEL))
  notes <- c(fit$error, fit$warnings)
  lower <- upper <- NA_real_
  if (!is.null(fit$value)) {
    if (is.finite(fit$value$value) && abs(fit$value$value - formula_value) > 1e-8) {
      stop("Internal check failed: ICC formula and irr::icc disagree.")
    }
    lower <- fit$value$lbound; upper <- fit$value$ubound
    if (length(lower) != 1L || !is.finite(lower)) lower <- NA_real_
    if (length(upper) != 1L || !is.finite(upper)) upper <- NA_real_
  }
  if (length(unique(x)) < 2L || length(unique(y)) < 2L) {
    notes <- c(notes, "At least one occurrence is constant; inspect the matrix.")
  }
  # Retain negative estimates; clipping to zero would hide poor agreement.
  list(icc = formula_value, lower = lower, upper = upper,
       note = paste(unique(notes), collapse = " | "))
}

analyze_pair <- function(x, y) {
  stopifnot(length(x) == length(y), !anyNA(x), !anyNA(y))
  tab <- table(factor(x, levels = RESPONSE_LEVELS), factor(y, levels = RESPONSE_LEVELS),
               dnn = c("earlier_occurrence", "later_occurrence"))
  n <- length(x)
  rho <- polychoric_stat(tab); ic <- icc_stat(x, y)
  exact <- if (n) sum(diag(tab)) / n else NA_real_
  stopifnot(sum(tab) == n)
  if (n) stopifnot(abs(exact - mean(x == y)) < 1e-12)
  list(table = tab, n = n, exact = exact, rho = rho, icc = ic,
       mean1 = if (n) mean(x) else NA_real_, mean2 = if (n) mean(y) else NA_real_,
       bias = if (n) mean(y - x) else NA_real_,
       mae = if (n) mean(abs(y - x)) else NA_real_,
       polarity_agreement = if (n) mean((x >= 3) == (y >= 3)) else NA_real_,
       dominant_cell_share = if (n) max(tab) / n else NA_real_)
}

records <- list(); summaries <- list(); cells <- list(); memberships <- list()
for (i in seq_len(nrow(pairs))) {
  p <- pairs[i, ]; j <- p$column1_index; k <- p$column2_index
  good_original <- stats::complete.cases(original[, c(j, k), drop = FALSE])
  good_corrected <- stats::complete.cases(corrected[, c(j, k), drop = FALSE])
  common <- good_original & good_corrected
  memberships[[i]] <- data.frame(pair_id = p$pair_id,
    fragrance = original[[ID_COLUMN]], available_original = good_original,
    available_corrected = good_corrected, matched = common, stringsAsFactors = FALSE)
  for (sample in c("matched", "available")) {
    for (dataset in c("Original", "Corrected")) {
      d <- if (dataset == "Original") original else corrected
      keep <- if (sample == "matched") common else {
        if (dataset == "Original") good_original else good_corrected
      }
      result <- analyze_pair(d[[j]][keep], d[[k]][keep])
      key <- paste(sample, p$pair_id, dataset, sep = "::")
      records[[key]] <- result
      summaries[[length(summaries) + 1L]] <- data.frame(
        sample = sample, pair_id = p$pair_id, item = p$item,
        occurrence1 = p$occurrence1, occurrence2 = p$occurrence2, dataset = dataset,
        n = result$n, excluded = nrow(d) - result$n,
        exact_agreement = result$exact, polychoric_r = result$rho$r,
        polychoric_method = result$rho$method, polychoric_boundary = result$rho$boundary,
        icc_a1 = result$icc$icc, icc_ci_lower = result$icc$lower,
        icc_ci_upper = result$icc$upper, mean_occurrence1 = result$mean1,
        mean_occurrence2 = result$mean2, mean_difference_2_minus_1 = result$bias,
        mean_absolute_difference = result$mae,
        polarity_agreement = result$polarity_agreement,
        dominant_cell_share = result$dominant_cell_share,
        occupied_rows = result$rho$occupied_rows,
        occupied_columns = result$rho$occupied_columns,
        zero_cells_out_of_16 = result$rho$zero_cells,
        polychoric_note = result$rho$note, icc_note = result$icc$note,
        stringsAsFactors = FALSE)
      ct <- as.data.frame(result$table, stringsAsFactors = FALSE)
      names(ct) <- c("earlier_response", "later_response", "count")
      ct$proportion_of_all_pairs <- if (result$n) ct$count / result$n else NA_real_
      ct$sample <- sample; ct$pair_id <- p$pair_id; ct$dataset <- dataset
      cells[[length(cells) + 1L]] <- ct[, c("sample", "pair_id", "dataset",
        "earlier_response", "later_response", "count", "proportion_of_all_pairs")]
    }
  }
  message("Analyzed ", p$pair_id, " (", i, "/", nrow(pairs), "): ", p$item)
}
statistics <- do.call(rbind, summaries)
write_csv(statistics, "agreement_statistics.csv")
write_csv(do.call(rbind, cells), "confusion_matrix_cells.csv")
write_csv(do.call(rbind, memberships), "sample_membership.csv", "diagnostics")
get_record <- function(sample, pair_id, dataset) {
  records[[paste(sample, pair_id, dataset, sep = "::")]]
}

comparison <- list()
for (sample in c("matched", "available")) {
  for (i in seq_len(nrow(pairs))) {
    id <- pairs$pair_id[i]
    x <- get_record(sample, id, "Original"); y <- get_record(sample, id, "Corrected")
    comparison[[length(comparison) + 1L]] <- data.frame(
      sample = sample, pair_id = id, item = pairs$item[i],
      n_original = x$n, n_corrected = y$n,
      polychoric_original = x$rho$r, polychoric_corrected = y$rho$r,
      delta_polychoric = y$rho$r - x$rho$r,
      icc_original = x$icc$icc, icc_corrected = y$icc$icc,
      delta_icc = y$icc$icc - x$icc$icc,
      exact_original = x$exact, exact_corrected = y$exact,
      delta_exact_percentage_points = 100 * (y$exact - x$exact),
      mae_original = x$mae, mae_corrected = y$mae,
      delta_mae = y$mae - x$mae, stringsAsFactors = FALSE)
  }
}
comparison <- do.call(rbind, comparison)
write_csv(comparison, "original_vs_corrected.csv")
write_csv(statistics[nzchar(statistics$polychoric_note) | nzchar(statistics$icc_note), ],
          "estimation_notes.csv", "diagnostics")

# ------------------ TRANSPARENT MAIN-TEXT ITEM SELECTION ---------------------
# These are ILLUSTRATIVE, data-selected examples, not an unbiased sample of items.
# The all-pair Supplement and comparison table prevent selective presentation.
choose_main <- function(comp) {
  d <- comp[comp$sample == MAIN_SAMPLE, , drop = FALSE]
  if (!is.null(MAIN_PAIR_IDS)) {
    if (length(MAIN_PAIR_IDS) != 3L || anyDuplicated(MAIN_PAIR_IDS) ||
        !all(MAIN_PAIR_IDS %in% d$pair_id)) {
      stop("MAIN_PAIR_IDS must contain exactly three distinct IDs from pair_map.csv.")
    }
    out <- d[match(MAIN_PAIR_IDS, d$pair_id), , drop = FALSE]
    out$selection_role <- c("User-selected example 1", "User-selected example 2",
                            "User-selected example 3")
    out$selection_rule <- "Manual override in MAIN_PAIR_IDS."
    return(out)
  }
  ok <- is.finite(d$icc_original) & is.finite(d$icc_corrected) &
        is.finite(d$polychoric_original) & is.finite(d$polychoric_corrected)
  pool <- d[ok & pmin(d$n_original, d$n_corrected) >= MIN_MAIN_N, , drop = FALSE]
  if (length(unique(pool$item)) < 3L) pool <- d[ok, , drop = FALSE]
  if (length(unique(pool$item)) < 3L) {
    stop("Fewer than three distinct items have estimable ICCs and polychoric correlations. ",
         "Specify MAIN_PAIR_IDS manually to include non-estimable examples.")
  }
  # A: largest ICC gain without a loss of exact agreement or polychoric r.
  candidates <- pool[pool$delta_icc > 0 & pool$delta_exact_percentage_points >= -1e-8 &
                       pool$delta_polychoric >= -1e-8, , drop = FALSE]
  rule_a <- "Largest ICC increase among items with no decrease in exact agreement or polychoric r."
  if (!nrow(candidates)) {
    candidates <- pool
    rule_a <- "Fallback: largest ICC change; inspect whether it actually improved."
  }
  one <- candidates[order(-candidates$delta_icc, -candidates$delta_exact_percentage_points,
                           candidates$pair_id)[1], , drop = FALSE]
  one$selection_role <- if (one$delta_icc > 0) "Improved agreement" else "Largest ICC change"
  one$selection_rule <- rule_a
  pool <- pool[pool$item != one$item, , drop = FALSE]
  # B: essentially unchanged; prefer a non-ceiling example with visible disagreement.
  candidates <- pool[abs(pool$delta_icc) <= 0.01 &
                       abs(pool$delta_exact_percentage_points) <= 1 &
                       abs(pool$delta_polychoric) <= 0.01, , drop = FALSE]
  rule_b <- paste("Absolute ICC and polychoric changes <= .01 and exact-agreement change",
                  "<= 1 percentage point; among these, lowest corrected exact agreement.")
  if (nrow(candidates)) {
    two <- candidates[order(candidates$exact_corrected, candidates$pair_id)[1], , drop = FALSE]
  } else {
    candidates <- pool
    two <- candidates[order(abs(candidates$delta_icc),
                             abs(candidates$delta_exact_percentage_points),
                             candidates$pair_id)[1], , drop = FALSE]
    rule_b <- "Fallback: smallest absolute ICC change."
  }
  two$selection_role <- "Little or no improvement"; two$selection_rule <- rule_b
  pool <- pool[pool$item != two$item, , drop = FALSE]
  # C: counterexample where exact matches increase but absolute-agreement ICC falls.
  candidates <- pool[pool$delta_icc < -0.02 &
                       pool$delta_exact_percentage_points > 1, , drop = FALSE]
  rule_c <- "Exact agreement increases > 1 percentage point but ICC decreases > .02; largest ICC decrease."
  if (nrow(candidates)) {
    three <- candidates[order(candidates$delta_icc, candidates$pair_id)[1], , drop = FALSE]
    role_c <- "More exact matches, lower ICC"
  } else {
    candidates <- pool[pool$delta_icc > 0.02 &
                         pool$delta_exact_percentage_points < -1, , drop = FALSE]
    if (nrow(candidates)) {
      three <- candidates[order(-candidates$delta_icc, candidates$pair_id)[1], , drop = FALSE]
      rule_c <- "ICC increases > .02 but exact agreement decreases > 1 percentage point; largest ICC increase."
      role_c <- "Higher ICC, fewer exact matches"
    } else {
      three <- pool[order(pool$delta_icc, pool$pair_id)[1], , drop = FALSE]
      rule_c <- "Fallback: lowest ICC change among remaining items."
      role_c <- "Contrasting example"
    }
  }
  three$selection_role <- role_c; three$selection_rule <- rule_c
  rbind(one, two, three)
}
selected <- choose_main(comparison)
selected$figure_row <- c("A", "B", "C")
write_csv(selected, "main_text_selection.csv")
message("Main figure: ", paste(selected$pair_id, collapse = ", "))

# ---------------------- PRINT-ORIENTED GRID GRAPHICS -------------------------
# Vector PDFs use physical dimensions. Raster copies preserve those dimensions.
# Every matrix uses the same grayscale mapping from 0% to 100% of its complete
# pairs. Percentages are NOT row-normalized. Diagonal borders mark exact agreement.
# Rows run 0,1,3,4 from TOP to BOTTOM; columns run 0,1,3,4 from LEFT to RIGHT.
gtext <- function(label, x, y, size = 10, face = "plain", just = "centre", rot = 0,
                  colour = "black") {
  grid::grid.text(label, x = grid::unit(x, "inches"), y = grid::unit(y, "inches"),
    just = just, rot = rot, gp = grid::gpar(fontsize = size, fontface = face,
    fontfamily = FONT_FAMILY, col = colour, lineheight = 1.05))
}
grect <- function(x, y, width, height, fill = NA, colour = "black", lwd = 0.5) {
  grid::grid.rect(x = grid::unit(x, "inches"), y = grid::unit(y, "inches"),
    width = grid::unit(width, "inches"), height = grid::unit(height, "inches"),
    gp = grid::gpar(fill = fill, col = colour, lwd = lwd))
}
wrap_text <- function(x, width) paste(strwrap(x, width = width), collapse = "\n")
fmt <- function(x) if (length(x) != 1L || !is.finite(x)) "NA" else sprintf("%.3f", x)
fmt_r <- function(x) {
  if (!is.finite(x$rho$r)) return("NA")
  if (x$rho$r > 0.9995) return(">0.999*")
  if (x$rho$r < -0.9995) return("<-0.999*")
  paste0(fmt(x$rho$r), if (x$rho$boundary) "*" else "")
}
fmt_exact <- function(x) if (is.finite(x$exact)) sprintf("%.1f%%", 100 * x$exact) else "NA"

matrix_grid <- function(x, centre_x, bottom, side, show_percent = FALSE,
                        count_size = 12, tick_size = 10, axis_labels = NULL) {
  tab <- x$table; cell <- side / length(RESPONSE_LEVELS); left <- centre_x - side / 2
  for (r in seq_along(RESPONSE_LEVELS)) {
    for (c in seq_along(RESPONSE_LEVELS)) {
      ncell <- as.numeric(tab[r, c]); p <- if (x$n) ncell / x$n else 0
      cx <- left + (c - 0.5) * cell
      cy <- bottom + (length(RESPONSE_LEVELS) - r + 0.5) * cell
      # A fixed, linear grayscale scale shared across panels and output files.
      grect(cx, cy, cell, cell, fill = grDevices::gray(1 - 0.75 * p),
            colour = grDevices::gray(0.55), lwd = 0.45)
      ink <- if (p > 0.60) "white" else "black"
      if (show_percent) {
        gtext(as.character(ncell), cx, cy + cell * 0.16, size = count_size,
              face = "bold", colour = ink)
        label <- if (x$n) sprintf("%.1f%%", 100 * p) else "--"
        gtext(label, cx, cy - cell * 0.22, size = count_size * 0.70, colour = ink)
      } else {
        gtext(as.character(ncell), cx, cy, size = count_size,
              face = if (r == c) "bold" else "plain", colour = ink)
      }
    }
  }
  # Draw diagonal borders last, so neighbouring cells cannot cover the outline.
  for (r in seq_along(RESPONSE_LEVELS)) {
    grect(left + (r - 0.5) * cell,
          bottom + (length(RESPONSE_LEVELS) - r + 0.5) * cell,
          cell, cell, lwd = 1.1)
  }
  offset <- if (show_percent) 0.13 else 0.09
  for (k in seq_along(RESPONSE_LEVELS)) {
    gtext(RESPONSE_LEVELS[k], left + (k - 0.5) * cell, bottom - offset, size = tick_size)
    gtext(RESPONSE_LEVELS[k], left - offset,
          bottom + (length(RESPONSE_LEVELS) - k + 0.5) * cell, size = tick_size)
  }
  if (!is.null(axis_labels)) {
    gtext(axis_labels[2], centre_x, bottom - 0.35, size = tick_size)
    gtext(axis_labels[1], left - 0.41, bottom + side / 2, size = tick_size, rot = 90)
  }
}

shade_key <- function(x, y, width = 1.4, label_size = 7.8) {
  n <- 80L
  for (j in seq_len(n)) {
    p <- (j - 0.5) / n
    grect(x + width * (j - 0.5) / n, y, width / n, 0.075,
          fill = grDevices::gray(1 - 0.75 * p), colour = NA)
  }
  gtext("0%", x, y - 0.10, label_size, just = "left")
  gtext("100%", x + width, y - 0.10, label_size, just = "right")
}

draw_main <- function() {
  grid::grid.newpage()
  w <- MAIN_WIDTH_IN; h <- MAIN_HEIGHT_IN
  gtext("Agreement between repeated fragrance ratings", w / 2, h - 0.37, 16, "bold")
  subtitle <- if (MAIN_SAMPLE == "matched") {
    "Illustrative items; the same complete-case fragrances are compared within each row"
  } else "Illustrative items; all available complete pairs within each dataset"
  gtext(subtitle, w / 2, h - 0.65, 10.3)
  centres <- c(w * 0.285, w * 0.735)
  row_height <- (h - 1.70) / 3
  side <- min(1.83, row_height - 1.24, w * 0.30)
  first_top <- h - 0.99
  for (i in 1:3) {
    top <- first_top - (i - 1) * row_height
    p <- pairs[match(selected$pair_id[i], pairs$pair_id), ]
    gtext(paste0(LETTERS[i], ".  ", p$item), w / 2, top, 12.0, "bold")
    for (j in 1:2) {
      dataset <- c("Original", "Corrected")[j]
      x <- get_record(MAIN_SAMPLE, p$pair_id, dataset)
      gtext(paste0(dataset, "  |  n = ", x$n, "  |  Exact = ", fmt_exact(x)),
            centres[j], top - 0.29, 11.2, "bold")
      gtext(paste0("Polychoric r = ", fmt_r(x), "    ICC(A,1) = ", fmt(x$icc$icc)),
            centres[j], top - 0.52, 10.9)
      matrix_grid(x, centres[j], bottom = top - 0.72 - side, side = side,
                  show_percent = TRUE, count_size = 13.5, tick_size = 10.5,
                  axis_labels = paste("Occurrence", c(p$occurrence1, p$occurrence2)))
    }
  }
  gtext("Cells show count and percentage of complete pairs; outlined diagonal cells indicate exact agreement.",
        w / 2, 0.47, 8.6)
  gtext("0 = definitely no; 1 = probably no; 3 = probably yes; 4 = definitely yes. ICC = absolute agreement, single measurement.",
        w / 2, 0.29, 8.2)
}

draw_supplement_page <- function(page, sample) {
  grid::grid.newpage()
  w <- SUPP_WIDTH_IN; h <- SUPP_HEIGHT_IN
  per_page <- SUPP_PAIR_COLUMNS * SUPP_ROWS
  n_pages <- ceiling(nrow(pairs) / per_page)
  gtext(paste0("Repeated-item agreement: all ", nrow(pairs), " pairs"),
        0.32, h - 0.25, 13.5, "bold", just = "left")
  subtitle <- if (sample == "matched") {
    "Matched complete cases: the same fragrances on both sides of each pair"
  } else "All available complete pairs within each dataset; sample sizes can differ"
  gtext(subtitle, 0.32, h - 0.49, 9.2, just = "left")
  gtext(paste0("Page ", page, " / ", n_pages), w - 0.32, h - 0.27,
        9.5, just = "right")
  left <- 0.24; top <- h - 0.78; bottom_margin <- 0.55
  group_w <- (w - 2 * left) / SUPP_PAIR_COLUMNS
  group_h <- (top - bottom_margin) / SUPP_ROWS
  side <- min((group_w / 2) - 0.36, group_h - 1.24, 1.35)
  if (side < 0.72) stop("Supplement layout is too dense; reduce columns/rows or enlarge page.")
  first <- (page - 1L) * per_page + 1L
  last <- min(page * per_page, nrow(pairs))
  for (ii in seq.int(first, last)) {
    local <- ii - first
    rr <- local %/% SUPP_PAIR_COLUMNS; cc <- local %% SUPP_PAIR_COLUMNS
    gx <- left + cc * group_w; gy <- top - rr * group_h
    p <- pairs[ii, ]
    question <- wrap_text(p$item, max(24L, floor(group_w * 14)))
    gtext(question, gx + group_w / 2, gy - 0.08, 8.9, "bold")
    gtext(paste0(p$pair_id, "  |  occurrences ", p$occurrence1, " vs ", p$occurrence2,
                 " (rows vs columns)"), gx + group_w / 2, gy - 0.36, 7.4)
    for (j in 1:2) {
      dataset <- c("Original", "Corrected")[j]
      cx <- gx + group_w * c(0.265, 0.765)[j]
      x <- get_record(sample, p$pair_id, dataset)
      gtext(dataset, cx, gy - 0.54, 8.6, "bold")
      gtext(paste0("r = ", fmt_r(x), "   ICC = ", fmt(x$icc$icc)), cx, gy - 0.70, 8.0)
      gtext(paste0("n = ", x$n, "   Exact = ", fmt_exact(x)), cx, gy - 0.85, 7.7)
      matrix_grid(x, cx, bottom = gy - 1.00 - side, side = side,
                  count_size = 9.3, tick_size = 7.7)
    }
  }
  gtext("Cells = counts; shading = share of complete pairs. All panels use the same scale. Diagonal borders = exact agreement.",
        0.32, 0.37, 7.7, just = "left")
  gtext("r = polychoric r; ICC = ICC(A,1). NA = not estimable; * = near-boundary r. Codes: 0 definitely no; 1 probably no; 3 probably yes; 4 definitely yes.",
        0.32, 0.19, 7.4, just = "left")
}

# ------------------------------- EXPORT -------------------------------------
open_device <- function(path, kind, width, height, dpi = PNG_DPI) {
  if (kind == "pdf") {
    if (capabilities("cairo")) {
      grDevices::cairo_pdf(filename = path, width = width, height = height,
                           family = FONT_FAMILY, onefile = TRUE, bg = "white")
    } else {
      grDevices::pdf(file = path, width = width, height = height, family = FONT_FAMILY,
                      onefile = TRUE, paper = "special", useDingbats = FALSE)
    }
  } else if (kind == "png") {
    if (capabilities("cairo")) {
      grDevices::png(path, width = width, height = height, units = "in", res = dpi,
                      type = "cairo", bg = "white")
    } else {
      grDevices::png(path, width = width, height = height, units = "in", res = dpi,
                      bg = "white")
    }
  } else if (kind == "tiff") {
    if (capabilities("cairo")) {
      grDevices::tiff(path, width = width, height = height, units = "in", res = dpi,
                       type = "cairo", compression = "lzw", bg = "white")
    } else {
      grDevices::tiff(path, width = width, height = height, units = "in", res = dpi,
                       compression = "lzw", bg = "white")
    }
  }
}
render_file <- function(filename, kind, width, height, draw, dpi = PNG_DPI) {
  path <- file.path(OUTPUT_DIR, "figures", filename)
  open_device(path, kind, width, height, dpi)
  on.exit(grDevices::dev.off(), add = TRUE)
  draw()
  invisible(path)
}
render_file("Main_text_confusion_matrices.pdf", "pdf", MAIN_WIDTH_IN, MAIN_HEIGHT_IN, draw_main)
render_file("Main_text_confusion_matrices.png", "png", MAIN_WIDTH_IN, MAIN_HEIGHT_IN, draw_main)
if (WRITE_TIFF && capabilities("tiff")) {
  render_file("Main_text_confusion_matrices.tiff", "tiff", MAIN_WIDTH_IN, MAIN_HEIGHT_IN,
              draw_main, dpi = TIFF_DPI)
} else if (WRITE_TIFF) {
  warning("This R installation does not support TIFF; PDF and PNG were still saved.")
}
for (sample in c("available", if (WRITE_MATCHED_SUPPLEMENT) "matched" else NULL)) {
  stem <- if (sample == "available") "Supplement_all_available" else "Supplement_matched_complete"
  npages <- ceiling(nrow(pairs) / (SUPP_PAIR_COLUMNS * SUPP_ROWS))
  render_file(paste0(stem, ".pdf"), "pdf", SUPP_WIDTH_IN, SUPP_HEIGHT_IN,
              function() for (page in seq_len(npages)) draw_supplement_page(page, sample))
  if (WRITE_SUPPLEMENT_PNGS) {
    for (page in seq_len(npages)) {
      render_file(sprintf("%s_page_%02d.png", stem, page), "png",
                  SUPP_WIDTH_IN, SUPP_HEIGHT_IN,
                  function() draw_supplement_page(page, sample), dpi = PNG_DPI)
    }
  }
}

# --------------------------- DOCUMENTATION / AUDIT ---------------------------
main_rows <- paste0(selected$figure_row, ": ", selected$item, " [", selected$pair_id,
                    "]. ", selected$selection_role, ".")
main_caption <- c(
  "MAIN-TEXT FIGURE CAPTION",
  "Agreement between repeated ratings of three illustrative fragrance items, before (left) and after (right) correction.",
  if (MAIN_SAMPLE == "matched") {
    "For each row, both matrices use exactly the fragrances with both occurrences observed in BOTH files. Sample size may differ across rows."
  } else "Each matrix uses all complete pairs in its dataset; sample sizes may differ between datasets.",
  "Rows and columns represent the earlier and later occurrence, respectively. Cells report counts and percentages of the matrix's complete pairs, not row percentages.",
  "The four categories are 0 (definitely no), 1 (probably no), 3 (probably yes), and 4 (definitely yes). Outlined diagonal cells indicate exact agreement.",
  paste0("Polychoric correlations use the ", POLYCHORIC_METHOD,
         " estimator with no continuity correction. ICC(A,1) is the two-way, absolute-agreement, single-measurement intraclass correlation on the original 0,1,3,4 scores."),
  "Items were selected after examining the data to illustrate different patterns, not as a random or unbiased sample of the repeated items.",
  main_rows,
  "If present, NA indicates that a statistic is not estimable, and * flags a polychoric estimate near its numerical boundary."
)
supp_caption <- c(
  "SUPPLEMENT FIGURE CAPTION",
  paste0("Confusion matrices for all ", nrow(pairs), " item pairs (",
         length(unique(pairs$item)), " distinct repeated questions)."),
  "Original and corrected matrices are adjacent. In Supplement_all_available, each matrix uses every complete pair available in that file.",
  "Supplement_matched_complete instead holds the fragrance sample fixed within each original/corrected comparison, matching the default main-text analysis.",
  "For questions administered three times, all three pairwise comparisons are shown under the default all_pairs setting. These comparisons are not statistically independent.",
  "Cells contain counts; shading is based on the fraction of all complete pairs and uses one fixed linear scale across every matrix.",
  "The occurrence numbers above each pair identify rows versus columns. Category order is 0,1,3,4, left to right and top to bottom.",
  "r = polychoric correlation; ICC = ICC(A,1), two-way absolute agreement for a single measurement. Exact = exact agreement. NA = not estimable; * = near-boundary polychoric estimate."
)
readme <- c(
  "FRAGRANCE REPEAT-AGREEMENT ANALYSIS", "",
  paste("Original file:", normalizePath(ORIGINAL_FILE)),
  paste("Corrected file:", normalizePath(CORRECTED_FILE)),
  paste("Generated:", format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z")), "",
  "KEY OUTPUTS",
  "figures/Main_text_confusion_matrices.pdf: vector, two columns x three rows, US-letter portrait.",
  "figures/Main_text_confusion_matrices.png: same figure at the configured PNG resolution.",
  "figures/Main_text_confusion_matrices.tiff: LZW-compressed 600-dpi figure when supported.",
  "figures/Supplement_all_available.pdf: every matrix, all available pairs in each file; six matrix columns by default.",
  "figures/Supplement_matched_complete.pdf: all pairs using matched complete cases, when enabled.",
  "tables/pair_map.csv: exact question names, column indices, occurrence numbers, and pair IDs.",
  "tables/agreement_statistics.csv: all coefficients, ICC confidence limits, sample sizes, exact agreement, and diagnostics.",
  "tables/original_vs_corrected.csv: side-by-side statistics and corrected-minus-original differences for both sample definitions.",
  "tables/main_text_selection.csv: selected rows, statistics, and explicit selection rules.",
  "tables/confusion_matrix_cells.csv: all cell counts and whole-matrix proportions.",
  "diagnostics/input_quality.csv: missingness and changes in response availability.",
  "diagnostics/sample_membership.csv: which fragrances enter each available/matched comparison.",
  "diagnostics/estimation_notes.csv: warnings, constant variables, and near-boundary estimates.", "",
  "MISSINGNESS",
  "Blanks are missing, not zero, not an additional category, and not an implicit response of 2.",
  "Available means pairwise complete within one file, not complete across all items.",
  "Matched means both occurrences are observed in both files for that particular pair.",
  "The corrected data change response availability, so all-available differences combine rating changes and sample-composition changes.",
  "Matched comparisons isolate changes on a fixed subset but do not establish unbiased generalization beyond complete cases.", "",
  "STATISTICAL DEFINITIONS",
  paste0("Polychoric r: polycor::polychor, method = ", POLYCHORIC_METHOD,
         ", maxcor = ", POLYCHORIC_MAX_ABS, ", no pseudocount/continuity correction."),
  "Polychoric correlation treats 0<1<3<4 as ordered categories and assumes an underlying bivariate normal model; sparse margins can make the estimate unstable.",
  "Unused marginal categories are dropped only for polychoric estimation. All four categories remain visible in every matrix.",
  "ICC(A,1): irr::icc(model='twoway', type='agreement', unit='single'), checked against its ANOVA formula.",
  "ICC = (MS_subject-MS_error)/(MS_subject+(k-1)*MS_error+k*(MS_occurrence-MS_error)/n), with k=2.",
  "ICC retains the stated 0,1,3,4 metric, so 1-to-3 disagreement is larger than 0-to-1 disagreement. Its interpretation assumes those numeric distances are meaningful.",
  "ICC confidence limits are model-based limits from irr, not tests of corrected-minus-original changes; ordinal, sparse data warrant caution.",
  "Negative ICC estimates are retained. Constant identical ratings can give 100% exact agreement while both correlations are undefined; undefined coefficients are not set to 1.",
  "Exact agreement = diagonal count / number of complete pairs. Polarity agreement compares no (0/1) versus yes (3/4).", "",
  "INTERPRETATION",
  "These analyses quantify agreement between repeated responses, not accuracy against an external truth criterion.",
  "Higher polychoric correlation does not necessarily imply greater absolute or exact agreement.",
  "If a correction procedure used the duplicate responses themselves, its improvement is not independent validation.",
  "The three main-text examples are deliberately illustrative. Use the all-pair tables/Supplement to describe the full distribution of changes.",
  "Do not treat the three pairings from a triplicated question as independent replications.", "",
  "REFERENCES / SOFTWARE DOCUMENTATION",
  "https://cran.r-project.org/web/packages/polycor/refman/polycor.html",
  "https://search.r-project.org/CRAN/refmans/irr/html/icc.html",
  "Olsson, U. (1979). Maximum likelihood estimation of the polychoric correlation coefficient. Psychometrika, 44, 443-460.",
  "McGraw, K. O., & Wong, S. P. (1996). Forming inferences about some intraclass correlation coefficients. Psychological Methods, 1, 30-46.",
  "Shrout, P. E., & Fleiss, J. L. (1979). Intraclass correlations: Uses in assessing rater reliability. Psychological Bulletin, 86, 420-428.", "",
  main_caption, "", supp_caption)
writeLines(readme, file.path(OUTPUT_DIR, "README_and_figure_captions.txt"), useBytes = TRUE)
writeLines(capture.output(sessionInfo()), file.path(OUTPUT_DIR, "diagnostics", "sessionInfo.txt"))
write_csv(data.frame(input = c("original", "corrected"),
                     md5 = unname(tools::md5sum(c(ORIGINAL_FILE, CORRECTED_FILE)))),
          "input_checksums.csv", "diagnostics")
saveRDS(list(settings = list(pairing = PAIRING, polychoric_method = POLYCHORIC_METHOD,
                             main_sample = MAIN_SAMPLE, response_levels = RESPONSE_LEVELS),
             pairs = pairs, statistics = statistics, comparison = comparison,
             selected = selected, records = records),
        file.path(OUTPUT_DIR, "diagnostics", "analysis_objects.rds"))
message("\nFinished. All outputs are in:\n", normalizePath(OUTPUT_DIR))
print(selected[, c("figure_row", "pair_id", "item", "n_original", "n_corrected",
                   "polychoric_original", "polychoric_corrected", "icc_original",
                   "icc_corrected", "exact_original", "exact_corrected")], row.names = FALSE)

#########################################################################################################
#########################################################################################################
# This concludes the confusion matrix analyses.
# Finally, I asked ChatGPT to check the accuracy of the 200 random fragrances before and after the 
# correction script. This accuracy check by ChatGPT was then reviewed by another LLM, Claude. Details 
# of this process are contained in a .zip file, but I paste the GUI portion of the conversations below.
# Thus, all below is normal text rather than R script.
#########################################################################################################
#########################################################################################################

# Historical conversation, commented for valid R parsing. FragranceDive now contains 35,449 records; older population descriptions in the original audit are superseded.
# MY INITIAL MESSAGE/PROMPT: 
#
# I have a large table of fragrance characteristic ratings, and I want to get an idea of how accurate they are. Each characteristic is in the form of a question, and the possible response options are: 0 = Definitely No 1 = Probably No 3 = Probably Yes 4 = Definitely Yes I've selected a random sample of 200 of these fragrances (out of 35,449) and am attaching their data. I would like you to do the absolute best you can, given your existing knowledge and all tools available to you, to check what percentage of these responses are correct. Use whatever approach you think will yield the most useful approximation of overall accuracy across the 35,449. There are some duplicate columns (questions) included for inter-rater reliability, but that reliability is irrelevant for your task. Treat duplicate columns like any other column when assessing accuracy, weighted no more or less than non-duplicate columns. Note also that many of these fragrances are quite obscure and you might be tempted to focus on well-known fragrances with the most easily-accessible data. Please don't do this; instead, try to make your accuracy checks as "random" as possible while still being systematic enough to make summary judgments. In presenting your results, please do your best to summarize them in one document (e.g., a .docx or markdown). If you provide additional (e.g., larger or more detailed) tables, please put them in a .zip as supplementary information separate from the single document I'm requesting.
#
# CHAT-GPT RESPONSE:
#
# [**Download the audit report (.docx)**](sandbox:/mnt/data/Fragrance_Rating_Accuracy_Audit.docx)
#
# [**Download the supplementary audit tables and reproducibility files (.zip)**](sandbox:/mnt/data/Fragrance_Accuracy_Audit_Supplement.zip)
#
# **My tentative working estimate is approximately 88% accuracy at the yes/no level—but the evidence does not support treating that as a firmly established accuracy rate for the entire table.**
#
# I randomly selected **1,200 ratings: six nonblank cells from every one of the 200 fragrances**, retaining duplicate question columns with the same treatment as other columns. Difficult-to-research fragrances were not replaced.
#
# The key results are:
#
# | Result                                                                         |                     Finding |
# | ------------------------------------------------------------------------------ | --------------------------: |
# | Ratings receiving a source-based assessment, including provisional judgments   |                433 of 1,200 |
# | Weighted agreement among those assessable ratings                              |                   **90.5%** |
# | Agreement in the stronger, direct-evidence subset                              | **87.4%** across 70 ratings |
# | Working estimate after adjusting for the file’s mix of 0, 1, 3 and 4 responses |                   **87.6%** |
# | Weighted proportion of the sample assessable                                   |                   **36.8%** |
#
# **The coverage limitation matters substantially.** Many questions remained unresolved, especially subtle odor impressions, performance, pricing, and ambiguous editions. Moreover, many negative reference judgments were provisional: a note not appearing in published descriptions is not proof that it is absent.
#
# The report explains exactly how those limitations affect the estimate, including sensitivity analyses showing how accuracy could fall if provisional negative judgments are wrong or unresolved ratings are less accurate. It also distinguishes **directional correctness**—No versus Yes—from the unvalidated distinction between “probably” and “definitely.”
#
# The ZIP includes the complete 1,200-cell audit trail, sources and explanations, coverage for all 200 fragrances, **44 flagged disagreements separated by evidence strength**, and a tested script that reproduces the reported calculations.
#
# MY RESPONSE MESSAGE/PROMPT:
#
# I have another version of that same file, with the same fragrances but many of the values updated (including some missing values filled in). Please run the same accuracy assessments on the attached file, being as consistent with your previous approach as possible (e.g., using the same cells). Provide the same audit document for this file, and include in that document some details about how and where the accuracy of this new document differs from the previous document.
#
# CHAT-GPT RESPONSE:
#
# [**Download the updated audit and comparison (.docx)**](sandbox:/mnt/data/Fragrance_Rating_Accuracy_Audit_Corrected_Comparison.docx)
#
# [**Download supplementary tables and reproducibility files (.zip)**](sandbox:/mnt/data/Fragrance_Corrected_Accuracy_Audit_Supplement.zip)
#
# **The corrected file performs substantially better against the original audit benchmark.** I reused the exact same 1,200 sampled cells, frozen reference judgments, question classifications and original sampling weights.
#
# ### Main comparison
#
# | Measure                                                                         | Original | Corrected |                     Change |
# | ------------------------------------------------------------------------------- | -------: | --------: | -------------------------: |
# | Weighted directional agreement on the **432 cells assessable in both versions** |    90.7% | **96.9%** | **+6.2 percentage points** |
# | Direct-evidence agreement, 70 cells                                             |    87.4% | **97.2%** |                +9.8 points |
# | Sensitivity to affirmative reference judgments                                  |    70.5% | **90.2%** |               +19.7 points |
# | Same-method, response-standardized working estimate*                            |    87.6% | **95.6%** |               +8.0 points* |
#
# *The standardized estimates use each version’s response distribution. They are assumption-dependent extrapolations, not a strictly paired comparison or independently established full-table accuracy.*
#
# The earlier **90.5%** result used 433 assessable cells. One previously disagreeing response is now blank; excluding that cell from **both** versions gives the comparable original baseline of **90.7%**. Deleting an answer is not counted as correcting it. 
#
# ### What changed
#
# Of the **44 previously flagged disagreements**, **30 now agree** with the unchanged reference, **13 still disagree**, and **one is now blank**. One previously agreeing cell now disagrees with a provisional reference. That leaves **14 current flags: two direct-tier and 12 provisional**.
#
# Most gains concern **note presence and prominence**. Seven of the nine original direct-tier disagreements resolve. Performance/evolution agreement remains **60.9%**, although that result rests on only eight assessable cells.
#
# ### Important limits
#
# The file fills **2,194 previously blank cells**, but those cells were outside the original sampling frame and **were not accuracy-tested in this comparison**. Evidence coverage remains approximately **36.8% of the weighted, still-answered original audit targets**.
#
# Also, because the previous audit was available before these revisions, this is a **repeated-benchmark comparison, not an independent holdout validation**. The approximately **96% working estimate** therefore requires substantial assumptions about unresolved judgments, provisional reference labels and newly filled cells.
#
# The document includes the paired comparisons, concrete corrections and remaining errors, response- and question-specific results, uncertainty intervals, and updated sensitivity analyses. The ZIP contains the complete cell-level comparison and a tested script reproducing the calculations.
#
# MY SUBSEQUENT MESSAGE TO A DIFFERENT LLM (CLAUDE):
#
# I asked ChatGPT to assess the accuracy of two fragrance data files (attached), and I'm attaching the conversation and results. Please review these results for accuracy. Check that the analysis was performed appropriately, and spot-check some proportion of judgments using your own investigation. Provide an assessment of the accuracy of ChatGPT's judgments and appropriateness of its approach and methods.
#
# CLAUDE RESPONSE:
#
# The short version:
#
# **The arithmetic is impeccable.** I recomputed every headline number from the raw CSVs — the sample reproduces exactly from the stated seed, all 1,200 cells match the source data, the weights sum correctly, and `Analyze.py` runs clean. 90.5%, 87.4%, 36.8%, 87.6%, the confusion matrices, the paired 432-cell comparison — all confirmed. No rounding errors or inconsistencies anywhere.
#
# **The second report's central finding is contaminated, and ChatGPT could have proven it but didnt.** Directional flip rates in the corrected file: 68.2% among cells the first audit flagged as wrong, 0.3% among cells it blessed as correct, 6.7% everywhere else. One flip among 389 "correct" cells where the base rate predicts 26 (p ≈ 5×10⁻¹¹). The corrections were made from the flag list. Simulating independent revision at the base rate gives 85.2%, not 96.9% — so essentially the whole 6.2-point "gain" is compliance with the benchmark, not improved annotation. ChatGPT warns about this in prose but never tests it, then headlines the number anyway.
#
# **The reference standard is narrower than presented.** 75% of assessed evidence is Parfumo alone, 87% Parfumo plus Fragrantica, with 182 distinct URLs for 200 fragrances — one source per fragrance, no triangulation. If your original ratings drew on the same public corpus, agreement is partly circular. And 88% of assessed cells are note questions, so the 88% headline extrapolates from note lists to price, longevity, projection, and "would consumers call this strange."
#
# **Spot-checks: 7 of 8 direct-tier judgments held.** I confirmed real errors in your data (Armani Code Sport vetiver, Loewe 7 Plata citrus, Dark Crude grapefruit, Metal Chypre bergamot, Blossom of Vetiver vanilla). One failed: Dancing With Strangers "fruity" is called direct-tier evidence, but the brand's note list has no fruit and the cited Parfumo page classifies it sweet-powdery — your "0" is probably right.
#
# The first report's core claim survives. The ~96% figures shouldn't be used.
#
# MY RESPONSE:
#
# Wait, this is interesting. Are you saying you strongly suspect the sequence was?:
#
# 1. Collect data
# 2. Audit Data
# 3. Run Correction Script that learned from the audit
# 4. Audit corrected data
# I can 100% assure you this is not what happened. The sequence of events was:
#
# 1. Collect data
# 2. Run correction script (which has outlier-detection)
# 3. Audit original data
# 4. Audit corrected data using the same approach and cells
# At least, I assume #4 happened that way - is there a reason to think it didn't use the same cells to check accuracy the second time?
#
# I'm attaching the correction script so you can see how the second ordering is feasible. I assure you both audits happened after the data was completely collected and corrected. Does this change your report and conclusion?
#
# CLAUDE RESPONSE:
#
# Yes, it changes the report. I was wrong, and I stated it too confidently.
#
# **Where I erred.** I found a pattern consistent with teaching-to-the-test and stopped there, without testing the obvious alternative: that same pattern is produced by *any* accurate independent error detector. Your script is one. Modeled properly, the evidence now favors your account:
#
# - Corrector sensitivity on audit-confirmed errors: 69.8%; false-positive rate on audit-confirmed-correct cells: 0.26%
# - File-wide flip rate 6.80% → **implied original error rate 9.4%**
# - Audit's own independent estimate: **9.5%**
#
# Two completely different routes to the original file's error rate, agreeing to a tenth of a point. And the simplest tell: patching the flag list would have changed ~30 cells. Your file changed 1,023 directionally and 4,293 total.
#
# The 0.3% figure is explained mechanically by the script. An RF predicting each column from the other 89 only flags cells inconsistent with their row — audit-confirmed-correct cells are consistent cells, so they never entered the review budget.
#
# **On your direct question:** yes, audit 2 used the same cells. Identical 1,200 row/column pairs, identical reference codes and tiers, and `Reference_frozen.json` is byte-identical across both ZIPs (SHA-256 matches the hash printed in report 1).
#
# **Two new findings replace the one I withdrew.**
#
# The first is the one I'd act on today. Your duplicate columns are the RF's own strongest predictors, so twin-disagreeing cells got flagged and resolved into agreement. Directional κ across duplicate pairs went **0.650 → 0.850**; cells whose twin disagreed were changed 42.7% of the time versus 10.9% for cells whose twin agreed. The corrected file's reliability statistic is an artifact of the correction. Archive the original for IRR, and hold duplicates out of the predictor set next run.
#
# The second: outlier detection can only find errors inconsistent with their row, so internally consistent systematic errors — a uniform misreading of "leading note," say — survive untouched and now make up a larger share of what remains. The audit can't see them either, since 88% of assessed cells are note questions scored against the same pyramids.
#
# The caveat that survives all this is method-sharing rather than leakage: your corrector is a web-search LLM reading Fragrantica and Parfumo, and ChatGPT's reference is 87% Fragrantica and Parfumo. The 96.9% is closer to an upper bound on what any database-mediated audit could find than to an accuracy estimate.