## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  fig.width  = 7,
  fig.height = 5,
  warning    = FALSE,
  message    = FALSE,
  eval       = FALSE
)

## ----simulate-synthetic-------------------------------------------------------
# library(dicepro)
# set.seed(2101L)
# 
# sim_data <- simulation(
#   loi        = "gauss",
#   scenario   = "hierarchical",
#   nSample    = 30,
#   nGenes     = 200,
#   nCellsType = 10,
#   sigma_bio  = 0.07,
#   sigma_tech = 0.07,
#   seed       = 2101L
# )
# 
# cat("Reference  :", dim(sim_data$W), "\n")
# cat("Proportions:", dim(sim_data$p), "\n")
# cat("Bulk       :", dim(sim_data$B), "\n")
# cat("Row sums   :", round(range(rowSums(sim_data$p)), 4), "\n")

## ----sanity-check-synthetic---------------------------------------------------
# bulk_clean <- as.matrix(sim_data$W) %*% t(as.matrix(sim_data$p))
# 
# plot(
#   bulk_clean[seq_len(500)],
#   as.matrix(sim_data$B)[seq_len(500)],
#   xlab = "Clean bulk (first 500 entries)",
#   ylab = "Noisy bulk",
#   pch  = 19, cex = 0.4, col = "#2c7bb6",
#   main = "Noise model: clean vs noisy bulk"
# )
# abline(0, 1, col = "firebrick", lwd = 1.5)

## ----run-dicepro-synthetic----------------------------------------------------
# out <- dicepro(
#   reference             = as.matrix(sim_data$W)[, -c(1,5,10)],
#   bulk                  = as.matrix(sim_data$B),
#   methodDeconv          = "FARDEEP",
#   W_prime               = 0,
#   bulkName              = "SimBulk",
#   refName               = "SimRef",
#   hp_max_evals          = 500,
#   algo_select           = "random",
#   output_path           = tempdir(),
#   hspaceTechniqueChoose = "gamma_dominant",
#   normalize             = FALSE
# )

## ----best-hp-synthetic--------------------------------------------------------
# cat("Best lambda :", out$hyperparameters$lambda, "\n")
# cat("Best gamma  :", out$hyperparameters$gamma,  "\n")
# cat("Loss        :", out$metrics$loss,           "\n")
# cat("Constraint  :", out$metrics$constraint,     "\n")

## ----plot-hyperopt-synthetic, fig.height=9------------------------------------
# out$plot_hyperopt

## ----plot-pareto-synthetic----------------------------------------------------
# out$plot

## ----compare-props-synthetic--------------------------------------------------
# true_prop <- as.matrix(sim_data$p)
# pred_prop <- as.matrix(out$H)
# 
# common_ct   <- intersect(colnames(pred_prop), colnames(true_prop))
# true_common <- true_prop[, common_ct, drop = FALSE]
# pred_common <- pred_prop[, common_ct, drop = FALSE]
# 
# r_overall <- cor(as.vector(true_common), as.vector(pred_common))
# cat(sprintf("Overall Pearson r: %.3f\n", r_overall))
# 
# plot(
#   as.vector(true_common),
#   as.vector(pred_common),
#   xlab = "True proportions",
#   ylab = "Predicted proportions",
#   pch  = 19, cex = 0.5, col = "#2c7bb6aa",
#   main = sprintf("True vs Predicted  (r = %.3f)", r_overall)
# )
# abline(0, 1, col = "firebrick", lwd = 1.5)

## ----per-ct-cor-synthetic, fig.height=4---------------------------------------
# ct_cors <- vapply(common_ct, function(ct) {
#   cor(true_common[, ct], pred_common[, ct])
# }, numeric(1L))
# 
# par(mar = c(5, 10, 3, 1))
# barplot(
#   sort(ct_cors),
#   horiz = TRUE, las = 1,
#   col   = ifelse(sort(ct_cors) > 0.7, "#2c7bb6", "#d7191c"),
#   xlab  = "Pearson r",
#   main  = "Per-cell-type correlation (synthetic)",
#   xlim  = c(-0.2, 1)
# )
# abline(v = 0.7, lty = 2, col = "gray40")

## ----perf-metrics-synthetic---------------------------------------------------
# perf <- makeTable1Tool(pred_mat = pred_common, obs_mat = true_common)
# knitr::kable(perf$Perf, digits = 3,
#              caption = "Performance metrics -- fully synthetic data")

## ----simulate-bluecode--------------------------------------------------------
# library(dicepro)
# 
# sim_bc <- simulation_bluecode(
#   nSample    = 30,
#   sigma_bio  = 0.15,
#   sigma_tech = 0.02,
#   seed       = 2101L
# )
# 
# cat("Reference  :", dim(sim_bc$W), "\n")  # nGenes x 34
# cat("Proportions:", dim(sim_bc$p), "\n")  # 30 x 34
# cat("Bulk       :", dim(sim_bc$B), "\n")  # nGenes x 30
# cat("Row sums   :", round(range(rowSums(sim_bc$p)), 4), "\n")
# 
# # Real cell-type names from BlueCode
# head(colnames(sim_bc$p))

## ----compartment-overview-----------------------------------------------------
# # Cell-type to compartment mapping (mirrors .bluecode_cell_groups)
# compartment_map <- c(
#   rep("Immune",      9),
#   rep("Stromal",     8),
#   rep("Endothelial", 3),
#   rep("Epithelial",  5),
#   rep("Muscle",      9)
# )
# names(compartment_map) <- colnames(sim_bc$p)
# 
# # Aggregate proportions by compartment for each sample
# comp_props <- t(apply(sim_bc$p, 1, function(row) {
#   tapply(row, compartment_map[names(row)], sum)
# }))
# 
# boxplot(
#   comp_props,
#   col  = c("#4393c3", "#74c476", "#fd8d3c", "#9ecae1", "#fb6a4a"),
#   ylab = "Compartment proportion",
#   main = "Simulated compartment proportions (BlueCode)",
#   las  = 2
# )

## ----sanity-check-bluecode----------------------------------------------------
# bulk_clean_bc <- as.matrix(sim_bc$W) %*% t(as.matrix(sim_bc$p))
# 
# plot(
#   bulk_clean_bc[seq_len(500)],
#   as.matrix(sim_bc$B)[seq_len(500)],
#   xlab = "Clean bulk (first 500 entries)",
#   ylab = "Noisy bulk",
#   pch  = 19, cex = 0.4, col = "#74c476",
#   main = "Noise model: clean vs noisy bulk (BlueCode)"
# )
# abline(0, 1, col = "firebrick", lwd = 1.5)

## ----run-dicepro-bluecode-----------------------------------------------------
# out_bc <- dicepro(
#   reference             = as.matrix(sim_bc$W),
#   bulk                  = as.matrix(sim_bc$B),
#   methodDeconv          = "FARDEEP",
#   W_prime               = 0,
#   bulkName              = "BlueCodeBulk",
#   refName               = "BlueCode",
#   hp_max_evals          = 100,
#   algo_select           = "random",
#   output_path           = tempdir(),
#   hspaceTechniqueChoose = "gamma_dominant",
#   normalize             = FALSE
# )

## ----best-hp-bluecode---------------------------------------------------------
# cat("Best lambda :", out_bc$hyperparameters$lambda, "\n")
# cat("Best gamma  :", out_bc$hyperparameters$gamma,  "\n")
# cat("Loss        :", out_bc$metrics$loss,            "\n")
# cat("Constraint  :", out_bc$metrics$constraint,      "\n")

## ----plot-hyperopt-bluecode, fig.height=9-------------------------------------
# out_bc$plot_hyperopt

## ----plot-pareto-bluecode-----------------------------------------------------
# out_bc$plot

## ----compare-props-bluecode---------------------------------------------------
# true_prop_bc <- as.matrix(sim_bc$p)
# pred_prop_bc <- as.matrix(out_bc$H)
# 
# common_ct_bc   <- intersect(colnames(pred_prop_bc), colnames(true_prop_bc))
# true_common_bc <- true_prop_bc[, common_ct_bc, drop = FALSE]
# pred_common_bc <- pred_prop_bc[, common_ct_bc, drop = FALSE]
# 
# r_overall_bc <- cor(as.vector(true_common_bc), as.vector(pred_common_bc))
# cat(sprintf("Overall Pearson r: %.3f\n", r_overall_bc))
# 
# plot(
#   as.vector(true_common_bc),
#   as.vector(pred_common_bc),
#   xlab = "True proportions",
#   ylab = "Predicted proportions",
#   pch  = 19, cex = 0.5, col = "#74c47699",
#   main = sprintf("True vs Predicted -- BlueCode  (r = %.3f)", r_overall_bc)
# )
# abline(0, 1, col = "firebrick", lwd = 1.5)

## ----per-ct-cor-bluecode, fig.height=6----------------------------------------
# ct_cors_bc <- vapply(common_ct_bc, function(ct) {
#   cor(true_common_bc[, ct], pred_common_bc[, ct])
# }, numeric(1L))
# 
# par(mar = c(5, 14, 3, 1))
# barplot(
#   sort(ct_cors_bc),
#   horiz = TRUE, las = 1,
#   col   = ifelse(sort(ct_cors_bc) > 0.7, "#2c7bb6", "#d7191c"),
#   xlab  = "Pearson r",
#   main  = "Per-cell-type correlation (BlueCode)",
#   xlim  = c(-0.2, 1)
# )
# abline(v = 0.7, lty = 2, col = "gray40")

## ----perf-metrics-bluecode----------------------------------------------------
# perf_bc <- makeTable1Tool(pred_mat = pred_common_bc, obs_mat = true_common_bc)
# knitr::kable(perf_bc$Perf, digits = 3,
#              caption = "Performance metrics -- BlueCode simulation")

## ----compare-strategies-------------------------------------------------------
# cat(sprintf(
#   "Overall Pearson r\n  Synthetic  : %.3f\n  BlueCode   : %.3f\n",
#   r_overall, r_overall_bc
# ))

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

