.libPaths(c(file.path(getwd(), ".local/R"), .libPaths()))
stopifnot(requireNamespace("MatrixGenerics", quietly = TRUE))
stopifnot(as.character(packageVersion("MatrixGenerics")) == "1.24.0",
          as.character(packageVersion("matrixStats")) == "1.5.0")
dir.create("evidence/outputs", showWarnings = FALSE, recursive = TRUE)

# Parse the upstream file; evaluate exactly its original top-level assignment.
# No substitutions, wrappers inside the function, or rewritten function body.
source_path <- "vendor/DESeq2/R/core.R"
expressions <- parse(source_path, keep.source = TRUE)
is_target <- function(x) {
  is.call(x) && identical(x[[1]], as.name("<-")) &&
    identical(x[[2]], as.name("estimateSizeFactorsForMatrix"))
}
indices <- which(vapply(expressions, is_target, logical(1)))
stopifnot(length(indices) == 1L)
upstream <- new.env(parent = baseenv())
eval(expressions[[indices]], envir = upstream)
estimate <- upstream$estimateSizeFactorsForMatrix
writeLines(deparse(estimate), "evidence/executed-function.txt")
writeLines(deparse(stats::median.default), "evidence/r-median.txt")

deterministic_counts <- function(n, m) {
  matrix(as.double((seq_len(n * m) * 7919) %% 100003 + 1), nrow = n)
}
cases <- list(
  odd = matrix(c(10,30,12,100,80,20,50,35,220,25,40,110,90,60,70), nrow = 5),
  even = deterministic_counts(6, 3),
  single_gene = matrix(c(8, 20, 80), nrow = 1),
  single_sample = matrix(c(4, 18, 90, 1000, 22), ncol = 1),
  tied = matrix(rep(c(10,20,40), each = 6), ncol = 3),
  zero_rows = rbind(deterministic_counts(5, 3), c(0,4,8), c(0,0,0), c(3,0,9)),
  two_samples = deterministic_counts(31, 2),
  four_samples = deterministic_counts(32, 4),
  large_counts = deterministic_counts(17, 3) + 2^40,
  thousand_genes = deterministic_counts(1001, 8)
)
manifest <- list()
for (name in names(cases)) {
  counts <- cases[[name]]
  scales <- 2^(seq_len(ncol(counts)) - 1)
  shifted <- sweep(counts, 2, scales, "*")
  global <- counts * 7
  retained <- which(apply(counts, 1, function(row) all(is.finite(row) & row > 0)))
  stopifnot(length(retained) > 0, all(counts >= 0), all(counts == floor(counts)),
            all(shifted < 2^53), all(global < 2^53))
  base_sf <- estimate(counts)
  scaled_sf <- estimate(shifted)
  global_sf <- estimate(global)
  filtered_sf <- estimate(counts[retained, , drop = FALSE])
  stopifnot(all(is.finite(c(base_sf, scaled_sf, global_sf, filtered_sf))),
            all(c(base_sf, scaled_sf, global_sf, filtered_sf) > 0))
  for (kind in c("base", "scaled", "global", "filtered")) {
    values <- get(paste0(if (kind == "base") "base" else kind, "_sf"))
    writeBin(as.double(values), paste0("evidence/outputs/", name, "-", kind, ".bin"),
             size = 8L, endian = "little")
  }
  write.table(counts, paste0("evidence/outputs/", name, "-counts.tsv"),
              sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)
  manifest[[name]] <- data.frame(case = name, rows = nrow(counts),
    samples = ncol(counts), retained = length(retained),
    scales = paste(scales, collapse = ","))
}
write.table(do.call(rbind, manifest), "evidence/cases.tsv", sep = "\t",
            row.names = FALSE, quote = FALSE)

# Boundary observations are empirical, outside the proved input contract.
observe <- function(label, thunk) {
  result <- tryCatch(suppressWarnings(thunk()), error = function(e) e)
  text <- if (inherits(result, "error")) paste0("error: ", conditionMessage(result)) else
    paste(format(result, digits = 17), collapse = ",")
  data.frame(case = label, observed = text)
}
boundaries <- list(
  observe("all_rows_have_zero", function() estimate(matrix(c(0,2,4,0), nrow=2))),
  observe("all_NA", function() estimate(matrix(NA_real_, nrow=2, ncol=3))),
  observe("mixed_NA", function() estimate(rbind(c(NA,2,3), c(1,2,4)))),
  observe("all_zero", function() estimate(matrix(0, nrow=2, ncol=3))),
  observe("empty_rows", function() estimate(matrix(numeric(), nrow=0, ncol=3))),
  observe("negative", function() estimate(rbind(c(-1,2,3), c(1,2,4)))),
  observe("positive_infinity", function() estimate(rbind(c(Inf,2,3), c(1,2,4)))),
  observe("empty_controls", function() estimate(cases$odd, controlGenes=integer())),
  observe("external_geoMeans", function() estimate(cases$odd, geoMeans=rep(1,5))),
  observe("poscounts", function() estimate(cases$zero_rows, type="poscounts"))
)
write.table(do.call(rbind, boundaries), "evidence/boundaries.tsv", sep="\t",
            row.names=FALSE, quote=TRUE)

# Distinguish exp(median(log ratios)) from median(ratios) for even length.
ev <- matrix(c(1,9,1,1), nrow=2)
ratios <- ev / exp(rowMeans(log(ev)))
up <- estimate(ev)
ordinary <- apply(ratios, 2, median)
stopifnot(abs(up[1] - ordinary[1]) > 0.1)
write.table(data.frame(sample=1:2, upstream=up, median_of_ratios=ordinary),
            "evidence/even-median-distinction.tsv", sep="\t", row.names=FALSE)
writeLines(capture.output(sessionInfo()), "evidence/r-session.txt")
cat("Executed unchanged pinned upstream function on", length(cases), "synthetic matrices.\n")