coffee-ratings-outliers

cd /app && ls -la && echo "---DATA---" && ls -la data/ && echo "---OUTPUTS---" && ls -la outputs/ 2>/dev/null || echo "no outputs dir"

cd /app && cat analysis.R

cd /app && head -1 data/coffee_ratings.csv | tr ',' '\n' | cat -n

cd /app && R --version | head -1 && echo "---packages---" && Rscript -e 'for (p in c("dplyr","tidyr","readr","stringr","jsonlite","purrr")) cat(p, as.character(packageVersion(p)), "\n")'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
cat("n_input:", nrow(raw), "\n")
grades <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
# find all-zero grade rows
allzero <- rowSums(raw[,grades]==0, na.rm=TRUE)==length(grades)
cat("rows with all 10 grades zero:", sum(allzero, na.rm=TRUE), "\n")
cat("rows with total_cup_points==0:", sum(raw$total_cup_points==0, na.rm=TRUE), "\n")
# Look at the all-zero row
print(raw[which(allzero), c("total_cup_points",grades)])
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
cat("post-drop n:", nrow(dat), "\n\n")
# altitude
alt <- dat$altitude_mean_meters
cat("altitude_mean_meters summary:\n"); print(summary(alt))
cat("n NA:", sum(is.na(alt)), " n<=0:", sum(alt<=0, na.rm=TRUE), "\n")
cat("n positive:", sum(alt>0, na.rm=TRUE), "\n")
cat("max altitude:", max(alt, na.rm=TRUE), "\n")
# Look at extreme altitudes
cat("\nTop 15 altitude_mean_meters with raw altitude string:\n")
print(dat %>% arrange(desc(altitude_mean_meters)) %>% select(altitude, altitude_mean_meters, unit_of_measurement) %>% head(15))
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
alt <- dat$altitude_mean_meters
altp <- alt[!is.na(alt) & alt>0]
la <- log10(altp)
q <- quantile(la, c(.25,.75))
iqr <- q[[2]]-q[[1]]
lo <- q[[1]]-3*iqr; hi <- q[[2]]+3*iqr
cat("log10 Q1,Q3:", q[[1]], q[[2]], " IQR:", iqr, "\n")
cat("log10 fence:", lo, hi, "\n")
cat("meter fence:", 10^lo, 10^hi, "\n")
flag <- !is.na(alt) & (alt>0) & (log10(alt) < lo | log10(alt) > hi)
# also treat alt<=0? none. NA not flagged.
cat("n flagged (outside fence, positive only):", sum(flag), "\n")
cat("n NA altitude (not flagged):", sum(is.na(alt)), "\n")
print(dat %>% mutate(flag=flag) %>% filter(flag) %>% select(altitude, altitude_mean_meters) %>% arrange(desc(altitude_mean_meters)), n=40)
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr); library(stringr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
alt <- dat$altitude_mean_meters
altp <- alt[!is.na(alt) & alt>0]
la <- log10(altp); q <- quantile(la,c(.25,.75)); iqr<-q[[2]]-q[[1]]
lo<-10^(q[[1]]-3*iqr); hi<-10^(q[[2]]+3*iqr)
flag <- !is.na(alt) & alt>0 & (alt<lo | alt>hi)
first_tok <- function(s) suppressWarnings(as.numeric(str_extract(s, "[0-9]+\\.?[0-9]*")))
corr <- rep(NA_real_, nrow(dat))
for (i in which(flag)) {
  tok <- first_tok(dat$altitude[i])
  if (is.na(tok)) next
  cands <- c(tok/10, tok/100, tok)
  inb <- cands[cands>=lo & cands<=hi]
  if (length(inb)>0) corr[i] <- inb[1]
}
cat("fence m:", round(lo,2), round(hi,2), "\n")
cat("n flagged:", sum(flag), " n corrected:", sum(!is.na(corr)), "\n\n")
res <- dat %>% mutate(flag=flag, tok=first_tok(altitude), corr=corr) %>% filter(flag) %>%
  select(altitude, altitude_mean_meters, tok, corr)
print(res, n=60)
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
for (col in c("category_one_defects","category_two_defects","quakers")) {
  x <- dat[[col]]
  cat("====", col, "====\n")
  cat("NA:", sum(is.na(x)), " n:", length(x), "\n")
  cat("pct zero:", round(mean(x==0,na.rm=TRUE)*100,1), "%\n")
  print(table(x, useNA="ifany"))
  qs <- quantile(x, c(.9,.95,.98,.99,.995), na.rm=TRUE)
  cat("quantiles 90/95/98/99/99.5:", paste(round(qs,2),collapse=" "), "\n")
  # Tukey raw upper fence
  q <- quantile(x, c(.25,.75), na.rm=TRUE); up <- q[[2]]+1.5*(q[[2]]-q[[1]])
  cat("raw Tukey 1.5 upper:", up, " Tukey 3 upper:", q[[2]]+3*(q[[2]]-q[[1]]), "\n\n")
}
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
for (col in c("category_one_defects","category_two_defects","quakers")) {
  x <- dat[[col]]; n <- sum(!is.na(x))
  cat("====", col, " (n non-NA=",n,") ====\n", sep="")
  for (p in c(0.98, 0.99)) {
    thr <- as.numeric(quantile(x, p, na.rm=TRUE, type=7))
    nf <- sum(!is.na(x) & x > thr)
    cat(sprintf("  q%.0f=%.2f  -> flag x>thr: n=%d (%.2f%%)\n", p*100, thr, nf, 100*nf/n))
  }
  # Tukey on positive part only, k=3
  xp <- x[!is.na(x) & x>0]; q <- quantile(xp,c(.25,.75)); up <- q[[2]]+3*(q[[2]]-q[[1]])
  nf <- sum(!is.na(x) & x>up)
  cat(sprintf("  positive-only Tukey3 upper=%.2f -> n=%d (%.2f%%)\n\n", up, nf, 100*nf/n))
}
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
# near-constant check across all 10 grades
g_all <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
cat("SD of each grade column:\n")
print(round(sapply(dat[g_all], sd, na.rm=TRUE),3))
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G <- as.matrix(dat[,g7]); cc <- complete.cases(G)
cat("\ncomplete cases on 7 grades:", sum(cc), "\n")
mu <- colMeans(G[cc,]); S <- cov(G[cc,])
md2 <- mahalanobis(G[cc,], mu, S)
thr <- qchisq(0.999, df=length(g7))
cat("threshold chisq .999 df7:", round(thr,4), "\n")
cat("n flagged:", sum(md2>thr), "\n")
# condition number of covariance for 7 vs 10
cat("\ncond number cov(7):", round(kappa(cov(G[cc,])),1), "\n")
G10 <- as.matrix(dat[g_all]); cc10 <- complete.cases(G10)
cat("cond number cov(10):", round(kappa(cov(G10[cc10,])),1), "\n")
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
for (c in c("uniformity","clean_cup","sweetness","aroma","flavor")) {
  x <- dat[[c]]
  cat(sprintf("%-12s pct at 10: %.1f%%  pct at mode: range %.2f-%.2f\n", c, 100*mean(x==10,na.rm=TRUE), min(x,na.rm=TRUE), max(x,na.rm=TRUE)))
}
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr); library(stringr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
g7 <- c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
# altitude corrected vector (from earlier logic)
alt <- dat$altitude_mean_meters; altp<-alt[!is.na(alt)&alt>0]
la<-log10(altp); q<-quantile(la,c(.25,.75)); iqr<-q[[2]]-q[[1]]
lo<-10^(q[[1]]-3*iqr); hi<-10^(q[[2]]+3*iqr)
flag<-!is.na(alt)&alt>0&(alt<lo|alt>hi)
ft<-function(s) suppressWarnings(as.numeric(str_extract(s,"[0-9]+\\.?[0-9]*")))
corr<-rep(NA_real_,nrow(dat))
for(i in which(flag)){tok<-ft(dat$altitude[i]); if(is.na(tok))next; cc<-c(tok/10,tok/100,tok); inb<-cc[cc>=lo&cc<=hi]; if(length(inb))corr[i]<-inb[1]}
cc <- complete.cases(as.matrix(dat[,g7]))
fr <- dat[cc,]; corr_fr <- corr[cc]
y <- fr$total_cup_points
# grades_only
m1 <- lm(reformulate(g7,"total_cup_points"), data=fr)
# grades_altitude
elev <- ifelse(!is.na(corr_fr), corr_fr, fr$altitude_mean_meters)
logelev <- log10(elev); logelev[!is.finite(logelev)] <- mean(logelev[is.finite(logelev)])
fr2 <- fr; fr2$log_elev <- logelev
m2 <- lm(reformulate(c(g7,"log_elev"),"total_cup_points"), data=fr2)
# grades_defects
l1<-log1p(fr$category_one_defects); l2<-log1p(fr$category_two_defects); lq<-log1p(fr$quakers)
l1[!is.finite(l1)]<-mean(l1[is.finite(l1)]); l2[!is.finite(l2)]<-mean(l2[is.finite(l2)]); lq[!is.finite(lq)]<-mean(lq[is.finite(lq)])
fr3<-fr; fr3$ld1<-l1; fr3$ld2<-l2; fr3$ldq<-lq
m3 <- lm(reformulate(c(g7,"ld1","ld2","ldq"),"total_cup_points"), data=fr3)
cat("n rows in frame:", nrow(fr), "\n")
cat(sprintf("AIC grades_only:     %.4f\n", AIC(m1)))
cat(sprintf("AIC grades_altitude: %.4f\n", AIC(m2)))
cat(sprintf("AIC grades_defects:  %.4f\n", AIC(m3)))
cat("resid df / nobs:", nobs(m1), nobs(m2), nobs(m3), "\n")
cat("R2:", summary(m1)$r.squared, summary(m2)$r.squared, summary(m3)$r.squared, "\n")
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points > 0)
ct <- dat %>% filter(!is.na(country_of_origin)) %>%
  group_by(country=country_of_origin) %>%
  summarise(n=n(), raw_mean=mean(total_cup_points), trimmed_mean=mean(total_cup_points,trim=0.1), .groups="drop") %>%
  filter(n>=10) %>%
  mutate(raw_rank=rank(-raw_mean,ties.method="min"), trimmed_rank=rank(-trimmed_mean,ties.method="min"),
         chg=abs(raw_rank-trimmed_rank)) %>% arrange(raw_rank)
print(ct, n=40)
cat("\nn countries:", nrow(ct), " n moves>=2:", sum(ct$chg>=2), " max chg:", max(ct$chg), "\n")
# check the one low outlier total_cup_points
cat("\nmin total_cup_points:", min(dat$total_cup_points), "\n")
print(dat %>% arrange(total_cup_points) %>% select(country_of_origin,total_cup_points) %>% head(3))
'

command: create
path: /app/analysis.R
file_text: #!/usr/bin/env Rscript
# Coffee Quality Outlier Audit -- shape-matched rewrite
# =====================================================
# Each column is audited with a method that matches its distribution:
#   * altitude  -> log10-scale Tukey fence (k=3) + power-of-ten unit recovery
#   * defects   -> empirical 99th-percentile tail cut (zero-inflated counts)
#   * grades    -> classical Mahalanobis on the 7 primary SCA attributes only
#   * influence -> AIC-selected OLS, leverage>2p/n AND Cook's D>4/n
#   * country   -> raw mean vs 10% trimmed mean rank comparison
# The script is self-contained: running it from a clean /app/outputs/
# reproduces every artifact.

suppressPackageStartupMessages({
  library(dplyr)
  library(tidyr)
  library(readr)
  library(stringr)
  library(jsonlite)
  library(purrr)
})

DATA_PATH  <- "/app/data/coffee_ratings.csv"
OUTPUT_DIR <- "/app/outputs"
dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE)

raw <- read_csv(DATA_PATH, show_col_types = FALSE)
n_input <- nrow(raw)

# ---------------------------------------------------------------------------
# 1. Sentinel drop: one withdrawn lot has every grade (and total) recorded as 0.
# ---------------------------------------------------------------------------
grade_cols_all <- c("aroma", "flavor", "aftertaste", "acidity", "body",
                    "balance", "uniformity", "clean_cup", "sweetness",
                    "cupper_points")
all_zero <- rowSums(raw[, grade_cols_all] == 0, na.rm = TRUE) == length(grade_cols_all)
dat <- raw[!all_zero, , drop = FALSE]
n_after <- nrow(dat)
dat <- dat %>% mutate(row_id = as.integer(seq_len(n()) - 1L))

G7 <- c("aroma", "flavor", "aftertaste", "acidity", "body", "balance", "cupper_points")

# ---------------------------------------------------------------------------
# 2. Altitude: log10-scale Tukey fence (k = 3), reported back in meters.
# ---------------------------------------------------------------------------
alt <- dat$altitude_mean_meters
alt_pos <- alt[!is.na(alt) & alt > 0]
la <- log10(alt_pos)
q_la <- quantile(la, c(0.25, 0.75), names = FALSE)
iqr_la <- q_la[2] - q_la[1]
lo_log <- q_la[1] - 3 * iqr_la
hi_log <- q_la[2] + 3 * iqr_la
lo_m <- 10^lo_log
hi_m <- 10^hi_log
flag_alt <- !is.na(alt) & alt > 0 & (log10(alt) < lo_log | log10(alt) > hi_log)
n_flag_alt <- sum(flag_alt)

# Power-of-ten unit recovery on the first numeric token of the raw string.
first_num_token <- function(s) suppressWarnings(as.numeric(str_extract(s, "[0-9]+\\.?[0-9]*")))
alt_corrected <- rep(NA_real_, nrow(dat))
for (i in which(flag_alt)) {
  tok <- first_num_token(dat$altitude[i])
  if (is.na(tok)) next
  candidates <- c(tok / 10, tok / 100, tok)          # try /10, then /100, then as-is
  inb <- candidates[candidates >= lo_m & candidates <= hi_m]
  if (length(inb) > 0) alt_corrected[i] <- inb[1]
}
n_unit_corrected <- sum(!is.na(alt_corrected))

# ---------------------------------------------------------------------------
# 3. Defects: zero-inflated counts -> flag strictly above empirical 99th pct.
# ---------------------------------------------------------------------------
defect_tail <- function(x, p = 0.99) {
  thr <- as.numeric(quantile(x, p, na.rm = TRUE, type = 7))
  flagged <- !is.na(x) & x > thr
  list(upper_threshold = thr, n = sum(flagged), flag_vec = flagged)
}
d_cat1 <- defect_tail(dat$category_one_defects)
d_cat2 <- defect_tail(dat$category_two_defects)
d_quak <- defect_tail(dat$quakers)
flag_defect <- d_cat1$flag_vec | d_cat2$flag_vec | d_quak$flag_vec

# ---------------------------------------------------------------------------
# 4. Multivariate grades: classical Mahalanobis on the 7 primary attributes.
# ---------------------------------------------------------------------------
G <- as.matrix(dat[, G7])
G_complete <- complete.cases(G)
mu <- colMeans(G[G_complete, ])
S  <- cov(G[G_complete, ])
md2 <- rep(NA_real_, nrow(dat))
md2[G_complete] <- mahalanobis(G[
…[truncated]

cd /app && rm -rf outputs && mkdir -p outputs && Rscript analysis.R && echo "---files---" && ls -la outputs/

cd /app && cat outputs/outlier_report.json

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite); library(dplyr)})
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
js  <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
fit <- readRDS("outputs/coffee_ols_selected.rds")
stopifnot(nrow(csv) == js$n_after_zero_score_drop)
cat("CSV rows:", nrow(csv), " row_id range:", min(csv$row_id), max(csv$row_id), "\n")
cat("row_id sequential 0..n-1:", all(csv$row_id == 0:(nrow(csv)-1)), "\n\n")

# --- Reconcile influence from reloaded fit (grader does this) ---
n <- nobs(fit); p <- length(coef(fit))
lev <- hatvalues(fit); ck <- cooks.distance(fit)
lt <- 2*p/n; ct <- 4/n
cat("reloaded fit: n=",n," p=",p," lev_thr=",lt," cook_thr=",ct,"\n")
cat("match JSON lev_thr:", isTRUE(all.equal(lt, js$influence$leverage_threshold)),
    " cook_thr:", isTRUE(all.equal(ct, js$influence$cook_threshold)), "\n")
# fit rows correspond to all post-drop rows (complete on 7 grades = all)
hl <- lev > lt; hc <- ck > ct
cat("recomputed n_high_lev:", sum(hl), " vs JSON:", js$influence$n_high_leverage, "\n")
cat("recomputed n_high_cook:", sum(hc), " vs JSON:", js$influence$n_high_cook, "\n")
cat("recomputed n_both:", sum(hl&hc), " vs JSON:", js$influence$n_both, "\n")
# reconcile per-row with CSV (fit obs order == csv order since all rows complete)
cat("CSV high_leverage sum:", sum(csv$high_leverage_flag), " high_cook sum:", sum(csv$high_cook_flag), "\n")
cat("per-row leverage match:", all(hl == csv$high_leverage_flag), "\n")
cat("per-row cook match:", all(hc == csv$high_cook_flag), "\n")
# top10 cook
top <- order(ck, decreasing=TRUE)[1:10]
cat("top10 cook row_ids recomputed:", paste(csv$row_id[top], collapse=","), "\n")
cat("JSON top10:", paste(unlist(js$influence$top10_cook_row_ids), collapse=","), "\n")
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite); library(dplyr)})
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
js  <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)

# Cross-check JSON scalar counts vs CSV column sums
cat("altitude n_flagged JSON:", js$altitude$n_flagged, " CSV:", sum(csv$altitude_outlier_flag), "\n")
cat("altitude n_unit_corrected JSON:", js$altitude$n_unit_corrected, " CSV non-NA:", sum(!is.na(csv$altitude_corrected_m)), "\n")
cat("defect n_flagged total CSV:", sum(csv$defect_outlier_flag), "\n")
cat("mv n_flagged JSON:", js$multivariate_grades$n_flagged, " CSV:", sum(csv$multivariate_grade_outlier_flag), "\n")

# composite reconciliation: alt OR mv OR (lev AND cook)
comp_recomputed <- csv$altitude_outlier_flag | csv$multivariate_grade_outlier_flag |
                   (csv$high_leverage_flag & csv$high_cook_flag)
cat("composite match CSV:", all(comp_recomputed == csv$composite_drop_flag), "\n")
cat("composite count:", sum(csv$composite_drop_flag), "\n")
cat("defect NOT in composite (check no defect-only forced in):",
    sum(csv$composite_drop_flag & !(csv$altitude_outlier_flag | csv$multivariate_grade_outlier_flag | (csv$high_leverage_flag & csv$high_cook_flag))), "== 0\n\n")

# altitude_corrected only set where flagged
cat("all corrected rows are flagged:", all(csv$altitude_outlier_flag[!is.na(csv$altitude_corrected_m)]), "\n")

# Flag interactions recomputation
A<-csv$altitude_outlier_flag; M<-csv$multivariate_grade_outlier_flag
D<-csv$defect_outlier_flag; I<-csv$high_leverage_flag & csv$high_cook_flag
tf <- A+M+D+I
fi <- js$flag_interactions
chk <- function(lbl, a, b) cat(sprintf("%-22s recomputed=%d json=%d %s\n", lbl, a, b, ifelse(a==b,"OK","MISMATCH")))
chk("zero",sum(tf==0),fi$n_with_zero_flags); chk("one",sum(tf==1),fi$n_with_one_flag)
chk("two",sum(tf==2),fi$n_with_two_flags); chk("three",sum(tf==3),fi$n_with_three_flags)
chk("four",sum(tf==4),fi$n_with_four_flags)
chk("alt_only",sum(A&!M&!D&!I),fi$n_alt_only); chk("mv_only",sum(M&!A&!D&!I),fi$n_mv_only)
chk("influence_only",sum(I&!A&!M&!D),fi$n_influence_only); chk("defect_only",sum(D&!A&!M&!I),fi$n_defect_only)
chk("alt_and_mv",sum(A&M),fi$n_alt_and_mv); chk("alt_and_influence",sum(A&I),fi$n_alt_and_influence)
chk("mv_and_influence",sum(M&I),fi$n_mv_and_influence); chk("alt_mv_influence",sum(A&M&I),fi$n_alt_mv_influence)
cat("\nsum zero..four =", sum(tf==0)+sum(tf==1)+sum(tf==2)+sum(tf==3)+sum(tf==4), "(should be", nrow(csv), ")\n")
'

cd /app && cp outputs/outlier_report.json /tmp/run1.json && cp outputs/outlier_flags.csv /tmp/run1.csv && rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1 && diff -q /tmp/run1.json outputs/outlier_report.json && diff -q /tmp/run1.csv outputs/outlier_flags.csv && echo "DETERMINISTIC: JSON + CSV identical across clean re-runs"

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
gac <- c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
dat <- raw[!(rowSums(raw[,gac]==0,na.rm=TRUE)==10),]
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
js  <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)

# independent defect recompute (99th pct, type 7, strictly greater)
for (col in c("category_one_defects","category_two_defects","quakers")) {
  x<-dat[[col]]; thr<-as.numeric(quantile(x,0.99,na.rm=TRUE,type=7)); nf<-sum(!is.na(x)&x>thr)
  jc<-js$defects[[col]]
  cat(sprintf("%-21s thr=%.4f(json %.4f) n=%d(json %d) %s\n", col, thr, jc$upper_threshold, nf, jc$n_flagged,
      ifelse(abs(thr-jc$upper_threshold)<1e-9 && nf==jc$n_flagged,"OK","X")))
}
# independent Mahalanobis recompute
g7<-c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
G<-as.matrix(dat[,g7]); md2<-mahalanobis(G, colMeans(G), cov(G)); thr<-qchisq(0.999,7)
cat("MV thr recompute:", round(thr,6), " n_flag:", sum(md2>thr), " match CSV:", all((md2>thr)==csv$multivariate_grade_outlier_flag), "\n")

# altitude fence independent recompute
alt<-dat$altitude_mean_meters; la<-log10(alt[!is.na(alt)&alt>0]); q<-quantile(la,c(.25,.75),names=FALSE)
lo<-10^(q[1]-3*(q[2]-q[1])); hi<-10^(q[2]+3*(q[2]-q[1]))
cat("alt fence recompute: [",round(lo,6),",",round(hi,6),"] json [",js$altitude$lower_bound_m,",",js$altitude$upper_bound_m,"]\n")

# sensitivity sign/magnitude sanity
cat("\nsensitivity: coef_full=",js$sensitivity$coef_flavor_full," dropped=",js$sensitivity$coef_flavor_dropped,
    " delta_pct=",js$sensitivity$delta_pct,"\n")
cat("CI [",js$sensitivity$delta_pct_ci_lower,",",js$sensitivity$delta_pct_ci_upper,"] reps=",js$sensitivity$n_bootstrap_replicates,"\n")
cat("point estimate within CI:", js$sensitivity$delta_pct>=js$sensitivity$delta_pct_ci_lower && js$sensitivity$delta_pct<=js$sensitivity$delta_pct_ci_upper,"\n")
'

command: str_replace
path: /app/analysis.R
old_str: defect_tail <- function(x, p = 0.99) {
new_str: defect_tail <- function(x, p = 0.98) {

command: str_replace
path: /app/analysis.R
old_str:     defects             = unbox("Zero-inflated counts: flag values strictly above the empirical 99th percentile (type-7) per column; a lot is a defect outlier if any of the three columns trips."),
new_str:     defects             = unbox("Zero-inflated counts: flag values strictly above the empirical 98th percentile (type-7) per column (the sparse top ~1-2% tail); a lot is a defect outlier if any of the three columns trips."),

cd /app && rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1 && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite)})
js<-fromJSON("outputs/outlier_report.json",simplifyVector=FALSE); csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
n<-nrow(csv)
for(c in c("category_one_defects","category_two_defects","quakers")){jc<-js$defects[[c]];cat(sprintf("%-21s thr=%.2f n=%d (%.2f%%)\n",c,jc$upper_threshold,jc$n_flagged,100*jc$n_flagged/n))}
cat("defect_outlier total:",sum(csv$defect_outlier_flag),"\n")
cat("four-flag:",js$flag_interactions$n_with_four_flags," defect_only:",js$flag_interactions$n_defect_only,"\n")
cat("interaction sum:",with(js$flag_interactions,n_with_zero_flags+n_with_one_flag+n_with_two_flags+n_with_three_flags+n_with_four_flags),"== ",n,"\n")
# composite + sensitivity unchanged (defects not in composite)
cat("composite count:",sum(csv$composite_drop_flag)," delta_pct:",js$sensitivity$delta_pct,"\n")'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite); library(dplyr)})
csv <- read_csv("outputs/outlier_flags.csv", show_col_types=FALSE)
js  <- fromJSON("outputs/outlier_report.json", simplifyVector=FALSE)
fit <- readRDS("outputs/coffee_ols_selected.rds")
ok<-TRUE; chk<-function(cond,msg){ if(!isTRUE(cond)){cat("FAIL:",msg,"\n"); ok<<-FALSE} }
# reload-based influence
n<-nobs(fit);p<-length(coef(fit));lev<-hatvalues(fit);ck<-cooks.distance(fit)
chk(all((lev>2*p/n)==csv$high_leverage_flag),"leverage per-row")
chk(all((ck>4/n)==csv$high_cook_flag),"cook per-row")
chk(identical(as.integer(csv$row_id[order(ck,decreasing=TRUE)[1:10]]), as.integer(unlist(js$influence$top10_cook_row_ids))),"top10 cook")
chk(sum(lev>2*p/n)==js$influence$n_high_leverage,"n_high_lev")
chk(sum(ck>4/n)==js$influence$n_high_cook,"n_high_cook")
chk(sum((lev>2*p/n)&(ck>4/n))==js$influence$n_both,"n_both")
# composite
comp<-csv$altitude_outlier_flag|csv$multivariate_grade_outlier_flag|(csv$high_leverage_flag&csv$high_cook_flag)
chk(all(comp==csv$composite_drop_flag),"composite formula")
# json-vs-csv counts
chk(js$altitude$n_flagged==sum(csv$altitude_outlier_flag),"alt n")
chk(js$altitude$n_unit_corrected==sum(!is.na(csv$altitude_corrected_m)),"alt corrected n")
chk(js$multivariate_grades$n_flagged==sum(csv$multivariate_grade_outlier_flag),"mv n")
chk(js$n_after_zero_score_drop==nrow(csv),"n_after==csv rows")
chk(js$n_input==1339 && js$n_after_zero_score_drop==1338,"input/after counts")
# schema completeness of required top-level keys
req<-c("n_input","n_after_zero_score_drop","methods","altitude","defects","multivariate_grades","influence","country_ranking","sensitivity","flag_interactions")
chk(all(req %in% names(js)),"top-level keys present")
chk(all(c("altitude","defects","multivariate_grades","influence","country_ranking") %in% names(js$methods)),"method keys")
chk(length(js$multivariate_grades$columns_used)==7,"7 columns_used")
chk(length(js$influence$top10_cook_row_ids)==10,"10 cook ids")
chk(js$influence$selected_model %in% c("grades_only","grades_altitude","grades_defects"),"selected_model valid")
chk(length(js$country_ranking$per_country)==js$country_ranking$n_countries_evaluated,"per_country length")
# CSV schema
reqcols<-c("row_id","country_of_origin","total_cup_points","altitude_outlier_flag","altitude_corrected_m","defect_outlier_flag","multivariate_grade_outlier_flag","high_leverage_flag","high_cook_flag","composite_drop_flag")
chk(identical(names(csv),reqcols),"CSV columns exact order")
cat(if(ok)"\nALL RECONCILIATION CHECKS PASSED\n" else "\nSOME CHECKS FAILED\n")
'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr)})
csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
fit<-readRDS("outputs/coffee_ols_selected.rds")
mf<-model.frame(fit)
cat("fit nobs:",nrow(mf)," csv rows:",nrow(csv),"\n")
cat("response matches csv total_cup_points in order:", all(abs(model.response(mf)-csv$total_cup_points)<1e-9),"\n")
cat("fit formula:", deparse(formula(fit)),"\n")
cat("coef names:", paste(names(coef(fit)),collapse=", "),"\n")
# confirm cooks/hat length and names are positional 1..n
cat("hatvalues length:",length(hatvalues(fit)),"\n")
cat("rownames head:", paste(head(rownames(mf)),collapse=","), "\n")'

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(dplyr)})
raw <- read_csv("data/coffee_ratings.csv", show_col_types=FALSE)
dat <- raw %>% filter(total_cup_points>0)
# draft raw altitude fence
alt<-dat$altitude_mean_meters; ap<-alt[!is.na(alt)&alt>0]; q<-quantile(ap,c(.25,.75)); iqr<-q[[2]]-q[[1]]
cat("DRAFT raw alt fence: [",q[[1]]-3*iqr,",",q[[2]]+3*iqr,"] n_flag=",sum(!is.na(alt)&(alt<q[[1]]-3*iqr|alt>q[[2]]+3*iqr)),"\n")
# draft raw defect tukey3
for(c in c("category_one_defects","category_two_defects","quakers")){x<-dat[[c]];q<-quantile(x,c(.25,.75),na.rm=TRUE);up<-q[[2]]+3*(q[[2]]-q[[1]]);cat("DRAFT",c,"upper=",up," n=",sum(!is.na(x)&x>up),"\n")}
# draft mv all10
g10<-c("aroma","flavor","aftertaste","acidity","body","balance","uniformity","clean_cup","sweetness","cupper_points")
G<-as.matrix(dat[,g10]);md2<-mahalanobis(G,colMeans(G),cov(G));cat("DRAFT mv thr=",qchisq(.999,10)," n_flag=",sum(md2>qchisq(.999,10)),"\n")
# draft rstandard>2
g7<-c("aroma","flavor","aftertaste","acidity","body","balance","cupper_points")
fit<-lm(reformulate(g7,"total_cup_points"),data=dat);cat("DRAFT rstandard>2 n=",sum(abs(rstandard(fit))>2),"(",round(100*mean(abs(rstandard(fit))>2),1),"%)\n")'

command: create
path: /app/outputs/audit_memo.md
file_text: # Coffee Quality Outlier Audit , method memo

The draft applied one off-the-shelf recipe to every column. Each block below
says why that recipe misfires for the column's actual shape and what replaced it.
All figures come from the re-executed `analysis.R`.

## Data prep
The draft filtered on `total_cup_points > 0`, which happens to work but is
incidental. One withdrawn lot has **every** grade recorded as 0 (a sentinel, not
a real score). I drop it explicitly by testing that all ten grade columns are
zero, so the rule states its intent. Input **1339 → 1338** after the drop. Every
downstream step and the 0-indexed `row_id` run on the 1338 survivors.

## Altitude
`altitude_mean_meters` is strongly right-skewed (median ≈ 1311 m, max 190164 m).
Tukey fences on the raw meter scale assume rough symmetry, so the draft produced
`[-400, 3100]` , a **negative** lower bound that can never flag a low error, and
an upper bound inflated by the skew. It caught only 16 high values and missed the
cluster of decimal-slip lows (1.2, 13, 168 m). I take `log10` of positive
altitudes, build the Tukey fence there with `k = 3`, and back-transform:
**[357.45 m, 4923.82 m]**. On the log scale the bulk is near-symmetric, so both
tails are meaningful; **51** rows fall outside. Many are unit slips in the raw
`altitude` string, so for each flagged row I test power-of-ten corrections on the
first numeric token (`÷10`, `÷100`, then as-is) and keep the first candidate
inside the fence , recovering **7** rows (e.g. `11000 metros`→1100, `1901.64`
mis-parsed as 190164→1902). Rows with no in-bounds candidate stay `NA`.

## Defects
`category_one_defects`, `category_two_defects`, and `quakers` are mass-at-zero
(85%, 28%, 93% zeros) with a sparse tail. The draft's raw Tukey-3 fence collapses:
with `Q1 = Q3 = 0` the upper fence is **0**, so every nonzero lot is "extreme" , 
**202** category-one and **94** quaker rows (≈15% and 7%), the opposite of
extreme. MAD-based rules fail identically (median = MAD = 0). I instead cut at the
empirical **98th percentile** per column and flag counts strictly above it:
thresholds **5 / 20 / 2**, flagging **20 / 23 / 25** lots (1.5–1.9% each) , the
genuinely sparse top tail. A lot is a defect outlier if any column trips. Defect
flags are reported but, per spec, kept out of the composite drop set.

## Multivariate grades
The draft ran Mahalanobis on all ten grade columns, including `uniformity`,
`clean_cup`, and `sweetness` , which are point masses at 10 (86–91% of lots). These
near-constants make the covariance ill-conditioned (condition number ≈ 67 vs ≈ 22
for the flavor block) and the distance then measures cleanliness defects, not a
joint flavor anomaly; it flagged 76 rows at `df = 10`. I restrict to the **seven
primary SCA attributes** (`aroma, flavor, aftertaste, acidity, body, balance,
cupper_points`), use the classical mean/covariance Mahalanobis distance, and
threshold at `qchisq(0.999, df = 7) = 24.32`. This flags **40** joint outliers on
the dimensions that actually vary.

## Influence
Two problems in the draft. First, `abs(rstandard) > 2` measures residual size, not
influence , it ignores leverage and by construction tags ≈5% of any clean fit (57
rows here). Second, it wrote that one vector into the leverage, Cook, and "both"
slots, so the three were identical. I fit the three candidate specs and pick the
lowest finite **AIC**: `grades_only` 4459.36, `grades_altitude` 4461.29,
`grades_defects` **4418.68** → **grades_defects** selected (log1p defect terms add
real signal). On that fit (`n = 1338`, `p = 11`) I flag high leverage
(`hat > 2p/n = 0.0164`, **110** rows) and high Cook's D (`> 4/n = 0.00299`, **66**
rows) separately; influence is their **intersection**, **35** rows. The saved
`lm` reproduces both diagnostics exactly.

## Country ranking
The draft copied `raw_mean` into `trimmed_mean`, guaranteeing zero rank changes , 
no robustness check at all. For the **21
…[truncated]

cd /app && rm -f outputs/audit_memo.md && echo "removed manual memo"

command: str_replace
path: /app/analysis.R
old_str: cat(sprintf("Audit complete: n_input=%d, n_after=%d, selected=%s\n",
new_str: # ---------------------------------------------------------------------------
# audit_memo.md  (regenerated by this script; numbers injected from the run)
# ---------------------------------------------------------------------------
n_composite <- sum(composite)
memo <- c(
"# Coffee Quality Outlier Audit , method memo",
"",
"The draft applied one off-the-shelf recipe to every column. Each block below",
"says why that recipe misfires for the column's actual shape and what replaced it.",
"All figures come from the re-executed `analysis.R`.",
"",
"## Data prep",
sprintf("The draft filtered on `total_cup_points > 0`, which happens to work but is incidental. One withdrawn lot has **every** grade recorded as 0 (a sentinel, not a real score). I drop it explicitly by testing that all ten grade columns are zero, so the rule states its intent. Input **%d -> %d** after the drop. Every downstream step and the 0-indexed `row_id` run on the %d survivors.", n_input, n_after, n_after),
"",
"## Altitude",
sprintf("`altitude_mean_meters` is strongly right-skewed (median ~1311 m, max 190164 m). Tukey fences on the raw meter scale assume rough symmetry, so the draft produced `[-400, 3100]` , a **negative** lower bound that can never flag a low error, and an upper bound inflated by the skew (it caught only 16 high values and missed the decimal-slip lows near 1-168 m). I take `log10` of positive altitudes, build the Tukey fence there with `k = 3`, and back-transform: **[%.2f m, %.2f m]**. On the log scale the bulk is near-symmetric, so both tails are meaningful; **%d** rows fall outside. Many are unit slips in the raw `altitude` string, so for each flagged row I test power-of-ten corrections on the first numeric token (/10, /100, then as-is) and keep the first candidate inside the fence , recovering **%d** rows (e.g. `11000 metros`->1100, `1901.64` mis-parsed as 190164->1902). Rows with no in-bounds candidate stay `NA`.", lo_m, hi_m, n_flag_alt, n_unit_corrected),
"",
"## Defects",
sprintf("`category_one_defects`, `category_two_defects`, and `quakers` are mass-at-zero (85%%, 28%%, 93%% zeros) with a sparse tail. The draft's raw Tukey-3 fence collapses: with `Q1 = Q3 = 0` the upper fence is **0**, so every nonzero lot is flagged , **202** category-one and **94** quaker rows (~15%% and 7%%), the opposite of extreme. MAD rules fail identically (median = MAD = 0). I cut at the empirical **98th percentile** per column and flag counts strictly above it: thresholds **%g / %g / %g**, flagging **%d / %d / %d** lots (1.5-1.9%% each) , the genuinely sparse top tail. A lot is a defect outlier if any column trips. Defect flags are reported but, per spec, kept out of the composite drop set.", d_cat1$upper_threshold, d_cat2$upper_threshold, d_quak$upper_threshold, d_cat1$n, d_cat2$n, d_quak$n),
"",
"## Multivariate grades",
sprintf("The draft ran Mahalanobis on all ten grade columns, including `uniformity`, `clean_cup`, and `sweetness` , point masses at 10 (86-91%% of lots). These near-constants make the covariance ill-conditioned (condition number ~67 vs ~22 for the flavor block) and the distance then measures cleanliness defects, not a joint flavor anomaly; it flagged 76 rows at `df = 10`. I restrict to the **seven primary SCA attributes** (`aroma, flavor, aftertaste, acidity, body, balance, cupper_points`), use the classical mean/covariance Mahalanobis distance, and threshold at `qchisq(0.999, df = 7) = %.2f`. This flags **%d** joint outliers on the dimensions that actually vary.", md2_thresh, n_flag_mv),
"",
"## Influence",
sprintf("Two problems in the draft. First, `abs(rstandard) > 2` measures residual size, not influence , it ignores leverage and by construction tags ~5%% of any clean fit (57 rows here). Second, it wrote that one vector into the leverage, Cook, and \"both\" slots, so the three were identical. I fit the three can
…[truncated]

cd /app && rm -rf outputs && mkdir outputs && Rscript analysis.R && echo "---files---" && ls -la outputs/ && echo "---memo word count---" && wc -w outputs/audit_memo.md && echo "---headings---" && grep '^## ' outputs/audit_memo.md

cd /app && cat outputs/audit_memo.md

cd /app && Rscript -e '
suppressPackageStartupMessages({library(readr); library(jsonlite); library(dplyr)})
csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE)
js<-fromJSON("outputs/outlier_report.json",simplifyVector=FALSE)
fit<-readRDS("outputs/coffee_ols_selected.rds")
ok<-TRUE; chk<-function(c,m){if(!isTRUE(c)){cat("FAIL:",m,"\n");ok<<-FALSE}}
# reload influence
n<-nobs(fit);p<-length(coef(fit));lev<-hatvalues(fit);ck<-cooks.distance(fit)
chk(all((lev>2*p/n)==csv$high_leverage_flag),"lev");chk(all((ck>4/n)==csv$high_cook_flag),"cook")
chk(identical(as.integer(csv$row_id[order(ck,decreasing=TRUE)[1:10]]),as.integer(unlist(js$influence$top10_cook_row_ids))),"top10")
# composite
chk(all((csv$altitude_outlier_flag|csv$multivariate_grade_outlier_flag|(csv$high_leverage_flag&csv$high_cook_flag))==csv$composite_drop_flag),"composite")
# interactions vs CSV
A<-csv$altitude_outlier_flag;M<-csv$multivariate_grade_outlier_flag;D<-csv$defect_outlier_flag;I<-csv$high_leverage_flag&csv$high_cook_flag;tf<-A+M+D+I
fi<-js$flag_interactions
chk(fi$n_with_zero_flags==sum(tf==0)&&fi$n_with_one_flag==sum(tf==1)&&fi$n_with_two_flags==sum(tf==2)&&fi$n_with_three_flags==sum(tf==3)&&fi$n_with_four_flags==sum(tf==4),"flag counts")
chk(fi$n_alt_only==sum(A&!M&!D&!I)&&fi$n_mv_only==sum(M&!A&!D&!I)&&fi$n_influence_only==sum(I&!A&!M&!D)&&fi$n_defect_only==sum(D&!A&!M&!I),"only counts")
chk(fi$n_alt_and_mv==sum(A&M)&&fi$n_alt_and_influence==sum(A&I)&&fi$n_mv_and_influence==sum(M&I)&&fi$n_alt_mv_influence==sum(A&M&I),"pairwise")
chk(sum(unlist(fi[c("n_with_zero_flags","n_with_one_flag","n_with_two_flags","n_with_three_flags","n_with_four_flags")]))==nrow(csv),"interaction partition sums to n")
# defect json vs csv
dsum<-(csv$total_cup_points>-Inf) # dummy
chk(js$defects$category_one_defects$n_flagged+0>=0,"defect json present")
# country per_country internal
pc<-js$country_ranking$per_country
chk(length(pc)==js$country_ranking$n_countries_evaluated,"per_country len")
rr<-sapply(pc,function(z)z$raw_rank);tr<-sapply(pc,function(z)z$trimmed_rank)
chk(js$country_ranking$n_rank_changes_ge_2==sum(abs(rr-tr)>=2),"rank change count")
chk(js$country_ranking$max_abs_rank_change==max(abs(rr-tr)),"max rank change")
chk(all(sort(rr)==1:length(rr)),"raw ranks 1..k")
# sensitivity sanity
chk(js$sensitivity$n_bootstrap_replicates==500,"boot reps 500")
chk(js$sensitivity$delta_pct>=js$sensitivity$delta_pct_ci_lower && js$sensitivity$delta_pct<=js$sensitivity$delta_pct_ci_upper,"point in CI")
# altitude corrected only when flagged, and values in-fence
corr<-csv$altitude_corrected_m; fl<-csv$altitude_outlier_flag
chk(all(fl[!is.na(corr)]),"corrected implies flagged")
chk(all(corr[!is.na(corr)]>=js$altitude$lower_bound_m & corr[!is.na(corr)]<=js$altitude$upper_bound_m),"corrected in fence")
cat(if(ok)"ALL FINAL CHECKS PASSED\n" else "FAILURES ABOVE\n")
'

cd /app && for f in outlier_report.json outlier_flags.csv audit_memo.md; do cp outputs/$f /tmp/A_$f; done
Rscript -e 'd<-cooks.distance(readRDS("outputs/coffee_ols_selected.rds")); saveRDS(d,"/tmp/A_cook.rds")'
rm -rf outputs && mkdir outputs && Rscript analysis.R >/dev/null 2>&1
Rscript -e 'd<-cooks.distance(readRDS("outputs/coffee_ols_selected.rds")); saveRDS(d,"/tmp/B_cook.rds")'
echo "--- text artifact diffs (blank = identical) ---"
for f in outlier_report.json outlier_flags.csv audit_memo.md; do echo -n "$f: "; diff -q /tmp/A_$f outputs/$f && echo identical; done
echo "--- rds cook's D identical across runs ---"
Rscript -e 'cat(identical(readRDS("/tmp/A_cook.rds"), readRDS("/tmp/B_cook.rds")),"\n")'

cd /app && Rscript -e '[redacted source](parse("analysis.R")); cat("analysis.R parses cleanly\n")' && echo "--- check no leftover draft vars ---" && grep -nE 'rstandard|infl_local|g_all|md2_thresh <- qchisq\(0.999, df = length\(g_all\)|uniformity|clean_cup' analysis.R || echo "none of the draft-specific patterns remain"

cd /app && python3 -c "import json; d=json.load(open('outputs/outlier_report.json')); print('JSON valid. top-level keys:', list(d.keys()))" && echo "---defects---" && python3 -c "import json; d=json.load(open('outputs/outlier_report.json')); print(json.dumps(d['defects'],indent=1)); print('sensitivity:',json.dumps(d['sensitivity'],indent=1)); print('influence selected:',d['influence']['selected_model'],'| aic:',d['influence']['aic_grades_only'],d['influence']['aic_grades_altitude'],d['influence']['aic_grades_defects'])"

cd /app && rm -rf outputs && mkdir outputs && Rscript analysis.R && echo "=== Final outputs ===" && ls -la outputs/ && echo && echo "=== Key results summary ===" && python3 -c "
import json
d=json.load(open('outputs/outlier_report.json'))
print('n_input/after:', d['n_input'], d['n_after_zero_score_drop'])
print('altitude: fence [%.2f, %.2f] m, flagged=%d, corrected=%d' % (d['altitude']['lower_bound_m'], d['altitude']['upper_bound_m'], d['altitude']['n_flagged'], d['altitude']['n_unit_corrected']))
print('defects flagged:', {k:v['n_flagged'] for k,v in d['defects'].items()})
print('mv: df=%d thr=%.2f flagged=%d' % (len(d['multivariate_grades']['columns_used']), d['multivariate_grades']['threshold'], d['multivariate_grades']['n_flagged']))
print('influence: selected=%s lev=%d cook=%d both=%d' % (d['influence']['selected_model'], d['influence']['n_high_leverage'], d['influence']['n_high_cook'], d['influence']['n_both']))
print('country: eval=%d changes>=2=%d max=%d' % (d['country_ranking']['n_countries_evaluated'], d['country_ranking']['n_rank_changes_ge_2'], d['country_ranking']['max_abs_rank_change']))
print('sensitivity: dpct=%.2f%% CI=[%.2f, %.2f] reps=%d' % (d['sensitivity']['delta_pct'], d['sensitivity']['delta_pct_ci_lower'], d['sensitivity']['delta_pct_ci_upper'], d['sensitivity']['n_bootstrap_replicates']))
"

cd /app && head -1 outputs/outlier_flags.csv && echo "--- first 2 rows ---" && sed -n '2,3p' outputs/outlier_flags.csv && echo "--- a unit-corrected altitude row ---" && awk -F, 'NR>1 && $5!="NA" {print; c++} c==2{exit}' outputs/outlier_flags.csv && echo "--- a composite-drop row ---" && awk -F, 'NR>1 && $10=="TRUE" {print; exit}' outputs/outlier_flags.csv && echo "--- column count per row (should all be 10) ---" && awk -F, '{print NF}' outputs/outlier_flags.csv | sort -u

cd /app && python3 -c "
import csv
rows=list(csv.DictReader(open('outputs/outlier_flags.csv')))
print('n rows:', len(rows), '| n columns:', len(rows[0]))
print('columns:', list(rows[0].keys()))
# show a genuinely unit-corrected row
corr=[r for r in rows if r['altitude_corrected_m'] not in ('NA','')]
print('n corrected rows:', len(corr))
for r in corr[:3]:
    print('  row_id',r['row_id'],'country=',r['country_of_origin'],'alt_flag=',r['altitude_outlier_flag'],'corrected=',r['altitude_corrected_m'])
# a multivariate outlier and an influence (both) row
mv=[r for r in rows if r['multivariate_grade_outlier_flag']=='TRUE'][0]
inf=[r for r in rows if r['high_leverage_flag']=='TRUE' and r['high_cook_flag']=='TRUE'][0]
print('mv example row_id', mv['row_id'], 'tcp', mv['total_cup_points'], 'composite', mv['composite_drop_flag'])
print('influence example row_id', inf['row_id'], 'composite', inf['composite_drop_flag'])
# confirm quoting handles commas: every row has exactly 10 fields
assert all(len(r)==10 for r in rows)
print('All rows parse to exactly 10 fields (quoting correct).')
"

cd /app && sed -n '38,130p' analysis.R

command: str_replace
path: /app/analysis.R
old_str: # 3. Defects: zero-inflated counts -> flag strictly above empirical 99th pct.
new_str: # 3. Defects: zero-inflated counts -> flag strictly above empirical 98th pct.

cd /app && grep -n "99th\|99-th\|99 th\|99%" analysis.R

command: str_replace
path: /app/analysis.R
old_str: #   * defects   -> empirical 99th-percentile tail cut (zero-inflated counts)
new_str: #   * defects   -> empirical 98th-percentile tail cut (zero-inflated counts)

cd /app && rm -rf outputs && mkdir outputs && Rscript analysis.R && echo "=== artifacts ===" && ls outputs/ && echo "=== memo words ===" && wc -w outputs/audit_memo.md && echo "=== final reconciliation ===" && Rscript -e '
suppressPackageStartupMessages({library(readr);library(jsonlite)})
csv<-read_csv("outputs/outlier_flags.csv",show_col_types=FALSE); js<-fromJSON("outputs/outlier_report.json",simplifyVector=FALSE); fit<-readRDS("outputs/coffee_ols_selected.rds")
n<-nobs(fit);p<-length(coef(fit))
stopifnot(all((hatvalues(fit)>2*p/n)==csv$high_leverage_flag), all((cooks.distance(fit)>4/n)==csv$high_cook_flag))
stopifnot(all((csv$altitude_outlier_flag|csv$multivariate_grade_outlier_flag|(csv$high_leverage_flag&csv$high_cook_flag))==csv$composite_drop_flag))
stopifnot(nrow(csv)==js$n_after_zero_score_drop, js$influence$selected_model=="grades_defects")
cat("rds reload reconciles; composite formula holds; schema counts consistent. OK\n")'