## ----knitr-config, include = FALSE--------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)

## ----setup--------------------------------------------------------------------
library(gDRcore)
library(gDRtestData)
library(gDRutils)
library(SummarizedExperiment)
library(BumpyMatrix)
library(data.table)

## ----load-data----------------------------------------------------------------
mae <- gDRutils::get_synthetic_data("finalMAE_small.qs2")
sa_se <- mae[["single-agent"]]
sa_se

## ----inspect-averaged---------------------------------------------------------
avg_cell <- BumpyMatrix::unsplitAsDataFrame(
  assay(sa_se, "Averaged"),
  row.field = "row", column.field = "column"
)
head(avg_cell[avg_cell$row == avg_cell$row[1] &
              avg_cell$column == avg_cell$column[1], ])

## ----minimal-fn---------------------------------------------------------------
# The simplest possible fit_fn: compute mean and SD of the response
summary_fn <- function(avg_dt) {
  list(
    x_mean = mean(avg_dt$x, na.rm = TRUE),
    x_sd = sd(avg_dt$x, na.rm = TRUE),
    n = NROW(avg_dt)
  )
}

## ----sa-basic-----------------------------------------------------------------
sa_out <- apply_fit(
  sa_se,
  fit_fn = summary_fn,
  data_type = "single-agent",
  output_assay = "custom_summary",
  fit_source = "demo"
)
assayNames(sa_out)

## ----sa-inspect---------------------------------------------------------------
summary_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(sa_out, "custom_summary"),
  row.field = "row", column.field = "column"
)
head(summary_df)

## ----sa-metrics---------------------------------------------------------------
# SE straight from the pipeline — already has a native "Metrics" assay
# (fit_source = "gDR")
"Metrics" %in% assayNames(sa_se)

# Apply a custom fit to the same assay, coexisting alongside gDR rows
custom_hill <- apply_fit(
  sa_se,
  fit_fn = fit_drug_response_metrics,
  data_type = "single-agent",
  output_assay = "Metrics",
  fit_source = "custom_hill"   # distinct key keeps native "gDR" rows intact
)

metrics_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(custom_hill, "Metrics"),
  row.field = "row", column.field = "column"
)
unique(metrics_df$fit_source)   # both "gDR" and "custom_hill"

## ----sa-idempotent------------------------------------------------------------
n_before <- NROW(BumpyMatrix::unsplitAsDataFrame(
  assay(sa_out, "custom_summary"),
  row.field = "row", column.field = "column"
))

# Call again — same fit_source, same data
sa_out2 <- apply_fit(
  sa_out,
  fit_fn = summary_fn,
  data_type = "single-agent",
  output_assay = "custom_summary",
  fit_source = "demo"
)
n_after <- NROW(BumpyMatrix::unsplitAsDataFrame(
  assay(sa_out2, "custom_summary"),
  row.field = "row", column.field = "column"
))

stopifnot(n_before == n_after)   # no duplicate rows
message("Row count before: ", n_before, " — after: ", n_after, " (no change)")

## ----sa-coexist---------------------------------------------------------------
extra_fn <- function(avg_dt) list(x_max = max(avg_dt$x, na.rm = TRUE))

sa_two <- sa_out |>
  apply_fit(extra_fn, "single-agent",
                   output_assay = "custom_summary",
                   fit_source = "extremes")

sources <- unique(BumpyMatrix::unsplitAsDataFrame(
  assay(sa_two, "custom_summary"),
  row.field = "row", column.field = "column"
)$fit_source)
message("fit_source values in assay: ", paste(sources, collapse = ", "))

## ----hill-params, eval = FALSE------------------------------------------------
# # Loosen the significance threshold or force the sigmoidal fit
# fit_drug_response_metrics(avg_dt, pcutoff = 0.1)
# fit_drug_response_metrics(avg_dt, force_fit = TRUE)
# # Different range for x_AOC_range computation
# fit_drug_response_metrics(avg_dt, range_conc = c(1e-3, 10))

## ----hill-ref-----------------------------------------------------------------
hill_out <- apply_fit(
  sa_se,
  fit_fn = fit_drug_response_metrics,
  data_type = "single-agent",
  output_assay = "custom_hill",
  fit_source = "hill_ref"
)

hill_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(hill_out, "custom_hill"),
  row.field = "row", column.field = "column"
)
head(hill_df[, c("row", "column", "normalization_type",
                  "ec50", "xc50", "h", "r2", "fit_type")])

## ----summary-fn---------------------------------------------------------------
# Aggregate: mean xc50 and flag whether both norm types fitted successfully
hill_summary_fn <- function(fit_dt) {
  list(
    mean_xc50 = mean(fit_dt$xc50, na.rm = TRUE),
    mean_r2 = mean(fit_dt$r2, na.rm = TRUE),
    all_converged = all(fit_dt$fit_type == "DRC3pHillFitModelFixS0", na.rm = TRUE)
  )
}

hill_with_summary <- apply_fit(
  sa_se,
  fit_fn = fit_drug_response_metrics,
  data_type = "single-agent",
  output_assay = "custom_hill",
  summary_fn = hill_summary_fn,
  summary_assay = "custom_hill_summary",
  fit_source = "hill_ref"
)
assayNames(hill_with_summary)

## ----summary-fn-inspect-------------------------------------------------------
sumdf <- BumpyMatrix::unsplitAsDataFrame(
  assay(hill_with_summary, "custom_hill_summary"),
  row.field = "row", column.field = "column"
)
head(sumdf[, c("row", "column", "mean_xc50", "mean_r2", "all_converged")])

## ----combo-apply-combo-scores-------------------------------------------------
# Use the small synthetic combo dataset which has both Averaged and Metrics
combo_mae <- gDRutils::get_synthetic_data("finalMAE_combo_matrix_small")
combo_name <- gDRutils::get_supported_experiments("combo")
combo_se_full <- combo_mae[[combo_name]]

# combo_se_full already has Metrics from fit_SE.combinations
combo_scored <- apply_combo_scores(combo_se_full)
assayNames(combo_scored)

## ----combo-apply-combo-scores-inspect-----------------------------------------
scores_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(combo_scored, "scores"),
  row.field = "row", column.field = "column"
)
scores_df[, c("row", "column", "normalization_type", "bliss_score", "hsa_score")]

## ----combo-se-----------------------------------------------------------------
# Build a minimal synthetic combination SE (Averaged only, no Metrics needed)
combo_dt <- data.table::CJ(
  row = c("DrugA", "DrugB"),
  column = "CellLine1",
  normalization_type = c("GR", "RV"),
  Concentration = c(0, 0.1, 1.0),
  Concentration_2 = c(0, 0.1, 1.0)
)
set.seed(42L)
combo_dt[, x := pmax(0.05,
  1 - 0.3 * Concentration / (Concentration + 0.5) -
      0.2 * Concentration_2 / (Concentration_2 + 0.5) +
      rnorm(.N, 0, 0.03))]

data_cols <- setdiff(names(combo_dt), c("row", "column"))
combo_bumpy <- BumpyMatrix::splitAsBumpyMatrix(
  combo_dt[, data_cols, with = FALSE],
  row = combo_dt$row, col = combo_dt$column
)
combo_se <- SummarizedExperiment(assays = list(Averaged = combo_bumpy))

## ----combo-bliss--------------------------------------------------------------
bliss_out <- apply_fit(
  combo_se,
  fit_fn = bliss_fit_fn,
  data_type = "combination",
  output_assay = "custom_bliss",
  fit_source = "bliss"
)

bliss_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(bliss_out, "custom_bliss"),
  row.field = "row", column.field = "column"
)
bliss_df[, c("row", "column", "normalization_type",
             "bliss_score", "bliss_excess_mean", "n_combo_points")]

## ----combo-hss----------------------------------------------------------------
hss_out <- apply_fit(
  combo_se,
  fit_fn = hss_fit_fn,
  data_type = "combination",
  output_assay = "custom_hss",
  fit_source = "hss"
)

hss_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(hss_out, "custom_hss"),
  row.field = "row", column.field = "column"
)
hss_df[, c("row", "column", "normalization_type",
           "hss_score", "hss_excess_mean")]

## ----multi-fit----------------------------------------------------------------
combo_multi <- apply_fits(
  combo_se,
  fit_fns = list(
    custom_bliss = bliss_fit_fn,
    custom_hss = hss_fit_fn
  ),
  data_type = "combination",
  fit_source = "synergy_panel"
)
assayNames(combo_multi)

## ----shared-precompute, eval = FALSE------------------------------------------
# # Each top-level name maps to an output assay; inner lists are the rows
# bliss_and_hss_combined <- function(dt) {
#   # Expensive step done ONCE per cell
#   sa1 <- dt[dt$Concentration_2 == 0 & dt$Concentration > 0, ]
#   sa2 <- dt[dt$Concentration == 0   & dt$Concentration_2 > 0, ]
# 
#   list(
#     custom_bliss = list(bliss_score = mean(sa1$x) - mean(sa2$x)), # simplified
#     custom_hss = list(hss_score = min(c(sa1$x, sa2$x)))
#   )
# }
# 
# apply_fits(
#   combo_se,
#   fit_fns = list(custom_bliss = bliss_and_hss_combined,
#                     custom_hss = bliss_and_hss_combined),
#   data_type = "combination",
#   fit_source = "shared"
# )

## ----pipe-chain---------------------------------------------------------------
result_se <- combo_se |>
  apply_fit(bliss_fit_fn, "combination",
                   output_assay = "custom_bliss",
                   fit_source = "bliss") |>
  apply_fit(hss_fit_fn, "combination",
                   output_assay = "custom_hss",
                   fit_source = "hss") |>
  apply_fit(
    function(dt) list(n_obs = NROW(dt)),
    "combination",
    output_assay = "combo_diagnostics",
    fit_source = "qc"
  )

assayNames(result_se)

## ----error-warn---------------------------------------------------------------
buggy_fn <- function(dt) {
  if (dt$normalization_type[1] == "GR") stop("GR not supported")
  list(x_rv = mean(dt$x, na.rm = TRUE))
}

se_partial <- withCallingHandlers(
  apply_fit(
    sa_se, buggy_fn, "single-agent",
    output_assay = "rv_only",
    fit_source = "rv_fn",
    on_error = "warn"
  ),
  warning = function(w) {
    message("[caught] ", conditionMessage(w))
    invokeRestart("muffleWarning")
  }
)
# Only RV rows are written; GR cells were skipped with a warning
rv_df <- BumpyMatrix::unsplitAsDataFrame(
  assay(se_partial, "rv_only"),
  row.field = "row", column.field = "column"
)
unique(rv_df$normalization_type)

## ----quick-ref, eval = FALSE--------------------------------------------------
# # Single fit → one output assay
# apply_fit(
#   se,
#   fit_fn,
#   data_type = "single-agent", # or "combination", "time-course"
#   slicing_cols = NULL, # NULL → data_type default
#   slicing_values = NULL, # NULL → all unique values found
#   input_assay = NULL, # NULL → data_type default ("Averaged")
#   output_assay, # REQUIRED — your assay name
#   summary_fn = NULL, # optional cell-level aggregator
#   summary_assay = NULL,
#   merge = "merge", # or "replace"
#   on_error = "warn", # or "stop"
#   fit_source                         # REQUIRED — upsert key tag
# )
# 
# # Multiple fits → one BumpyMatrix pass
# apply_fits(
#   se,
#   fit_fns, # named list: name = output assay, value = fit function
#   data_type = "single-agent",
#   fit_source,
#   ...
# )

## ----session-info-------------------------------------------------------------
sessionInfo()