# =============================================================================
# Title:
# - False Discovery Rate (FDR) Control for Spatial Raster Data. 
#
# Authors:
# - Oliver Gutiérrez-Hernández, Department of Geography, University of Málaga (UMA) (Corresponding author: olivergh@uma.es).
# - Luis V. García, Institute of Natural Resources and Agrobiology of Seville (IRNAS), Spanish National Research Council (CSIC).
#
# R environment and package versions:
# - R version 4.5.1 (2025-10-22)
# - CRAN packages:
#     • terra 1.8-80 — Raster and spatial data handling
#     • cp4p 0.3.6 — Adaptive FDR estimation and control procedures
#     • colorspace 2.1-0 — Provides hcl.colors() and Cividis palette
# - Bioconductor packages:
#     • multtest 2.64.0 — Implementation of step-up procedures (BH, BKY)
#     • qvalue 2.40.0 — Estimation of π₀ and FDR-based q-values
#     • limma 3.64.1 — General multiple testing framework (dependency for cp4p)
# - Base R packages (implicit): stats, graphics, grDevices, utils (all version 4.5.1)
#
# ⚙️ Note:
# The utility function `install_if_missing()` automatically checks for the presence of
# required packages and installs the most recent available version from CRAN or Bioconductor
# if they are not already installed. This ensures that the workflow remains reproducible
# while staying up to date with stable package releases.
# =============================================================================

install_if_missing <- function(package, bioc = FALSE) {
    if (!requireNamespace(package, quietly = TRUE)) {
        if (bioc) {
            if (!requireNamespace("BiocManager", quietly = TRUE)) {
                install.packages("BiocManager", dependencies = TRUE)
            }
            BiocManager::install(package, ask = FALSE, update = FALSE)
        } else {
            install.packages(package, dependencies = TRUE)
        }
    }
}

# --- Install required packages ---
install_if_missing("terra")
install_if_missing("cp4p")
install_if_missing("colorspace")  # ✅ Required for hcl.colors() and Cividis palette
install_if_missing("multtest", bioc = TRUE)
install_if_missing("qvalue", bioc = TRUE)
install_if_missing("limma", bioc = TRUE)

# --- Load packages ---
library(terra)
library(cp4p)
library(colorspace)  # ✅ Added for color palette support

# --- Load raster ---
raster_path <- "C:/data/p-values.tif"
test_raster <- rast(raster_path)
raster_base_name <- tools::file_path_sans_ext(basename(raster_path))
output_dir <- dirname(raster_path)

# --- Backup copy ---
original_copy_path <- file.path(output_dir, paste0(raster_base_name, "_original.tif"))
writeRaster(test_raster, original_copy_path, overwrite = TRUE)

# --- Extract and validate values ---
p_values <- values(test_raster)
valid_indices <- !is.na(p_values)
num_valid_pixels <- sum(valid_indices)
num_invalid_pixels <- sum(is.na(p_values))

# --- Adjust p-values ---
adjusted_p_values_bky <- rep(NA, length(p_values))
adjusted_p_values_bh  <- rep(NA, length(p_values))
adjusted_p_values_bky[valid_indices] <- adjust.p(p_values[valid_indices], pi0.method = "bky", alpha = 0.05)$adjp[, 2]
adjusted_p_values_bh[valid_indices]  <- p.adjust(p_values[valid_indices], method = "BH")

# --- Binary maps ---
binary_original <- ifelse(p_values <= 0.05, 1, 0)
binary_bky <- ifelse(adjusted_p_values_bky <= 0.05, 1, 0)
binary_bh  <- ifelse(adjusted_p_values_bh <= 0.05, 1, 0)
binary_original[!valid_indices] <- NA
binary_bky[!valid_indices] <- NA
binary_bh[!valid_indices] <- NA

# --- Create raster layers ---
raster_adjusted_bky <- setValues(test_raster, adjusted_p_values_bky)
raster_adjusted_bh  <- setValues(test_raster, adjusted_p_values_bh)
raster_binary_original <- setValues(test_raster, binary_original)
raster_binary_bky <- setValues(test_raster, binary_bky)
raster_binary_bh  <- setValues(test_raster, binary_bh)

# --- Summary statistics ---
num_rejected_original <- sum(binary_original == 1, na.rm = TRUE)
num_rejected_bky <- sum(binary_bky == 1, na.rm = TRUE)
num_rejected_bh  <- sum(binary_bh == 1, na.rm = TRUE)
perc_original <- (num_rejected_original / num_valid_pixels) * 100
perc_bky <- (num_rejected_bky / num_valid_pixels) * 100
perc_bh  <- (num_rejected_bh / num_valid_pixels) * 100

cat("\n--- Summary Statistics ---\n")
cat("Valid pixels:", num_valid_pixels, "\n")
cat("Significant pixels (raw p-values < 0.05):", num_rejected_original, "(", round(perc_original, 2), "%)\n")
cat("Significant pixels (FDR-BH < 0.05):", num_rejected_bh, "(", round(perc_bh, 2), "%)\n")
cat("Significant pixels (Adaptive FDR-BKY < 0.05):", num_rejected_bky, "(", round(perc_bky, 2), "%)\n")

# ---------------------------------------------------------------------------
# FIRST FIGURE (added): Raw p-values map + histogram (Cividis palette)
# ---------------------------------------------------------------------------

# Define color palette (Cividis)
pal_raw <- hcl.colors(100, "Cividis")

op <- par(mfrow = c(2,1), mar = c(4,5,3,6))

# 1) Raw p-values map
plot(test_raster,
     main = expression("Raw " * italic(p) * "-values"),
     col = pal_raw, zlim = c(0,1))

# 2) Histogram of raw p-values (10 intervals, same as map scale)
breaks_seq <- seq(0, 1, length.out = 11)
hist_vals <- hist(p_values[valid_indices],
                  breaks = breaks_seq,
                  plot = FALSE)

bar_colors <- pal_raw[as.numeric(cut(hist_vals$mids,
                                     breaks = seq(0, 1, length.out = length(pal_raw) + 1)))]

barplot(hist_vals$counts,
        names.arg = round(hist_vals$mids, 2),
        xlab = expression(italic(p)*"-values"),
        ylab = "Frequency",
        main = expression("Histogram of raw " * italic(p) * "-values"),
        col = bar_colors,
        border = "white",
        space = 0)

par(op)

# --- Threshold plots (FDR-BH and Adaptive FDR-BKY) ---
ordered_p_values <- sort(p_values[valid_indices])
m <- length(ordered_p_values)
q <- 0.05
thresh_bh <- (1:m) * q / m
res_bky <- adjust.p(ordered_p_values, pi0.method = "bky", alpha = q)
pi0_est <- res_bky$pi0
thresh_bky <- (1:m) * q / (m * pi0_est)

par(mfrow = c(1, 2), mar = c(5, 5, 4, 2) + 0.1, cex.lab = 1.2, cex.axis = 1.1, cex.main = 1.3)

plot(1:m, ordered_p_values, pch = 19, col = ifelse(ordered_p_values <= thresh_bh, 'blue', 'grey'),
     xlab = expression("Rank of " * italic(p) * "-values"), ylab = expression(italic(p)*"-values"),
     main = expression("FDR-BH threshold (" * italic(q) * " = 0.05)"), ylim = c(0, 1))
lines(1:m, thresh_bh, col = 'skyblue', lwd = 2)
grid()

plot(1:m, ordered_p_values, pch = 19, col = ifelse(ordered_p_values <= thresh_bky, 'blue', 'grey'),
     xlab = expression("Rank of " * italic(p) * "-values"), ylab = expression(italic(p)*"-values"),
     main = expression("Adaptive FDR-BKY threshold (" * italic(q) * " = 0.05)"),
     ylim = c(0, max(ordered_p_values, na.rm = TRUE)))
lines(1:m, thresh_bky, col = 'skyblue', lwd = 2)
grid()

# --- Bar plots ---
draw_bar <- function(sig, total, main_title) {
    bar_colors <- c("gray95", "blue")
    bar_labels <- rev(c("Discoveries", "Not discoveries"))
    bar_vals <- rev(c(sig, total - sig))
    bar_perc <- round(bar_vals / total * 100, 1)
    bar_text <- paste0(bar_vals, " (", bar_perc, "%)")
    bp <- barplot(bar_vals, names.arg = bar_labels, col = bar_colors,
                  main = main_title, xlab = "No. of tests", horiz = TRUE, las = 1, xlim = c(0, total),
                  xaxs = "i")
    text(x = bar_vals + total * 0.02, y = bp, labels = bar_text, pos = 4, cex = 1.2, col = "black")
}

par(mfrow = c(3, 2), mar = c(5, 5, 4, 2) + 0.1, cex.lab = 1.2, cex.axis = 1.1, cex.main = 1.3)
plot(raster_binary_original, main = expression("Significant (raw " * italic(p) * "-values, α = 0.05)"), col = c("gray95", "blue"), legend = FALSE)
draw_bar(num_rejected_original, num_valid_pixels, expression("Significant (raw " * italic(p) * "-values, α = 0.05)"))
plot(raster_binary_bh, main = expression("Discoveries (FDR-BH, " * italic(q) * " = 0.05)"), col = c("gray95", "blue"), legend = FALSE)
draw_bar(num_rejected_bh, num_valid_pixels, expression("Discoveries (FDR-BH, " * italic(q) * " = 0.05)"))
plot(raster_binary_bky, main = expression("Discoveries (Adaptive FDR-BKY, " * italic(q) * " = 0.05)"), col = c("gray95", "blue"), legend = FALSE)
draw_bar(num_rejected_bky, num_valid_pixels, expression("Discoveries (Adaptive FDR-BKY, " * italic(q) * " = 0.05)"))

# --- File output ---
writeRaster(raster_adjusted_bh,  file.path(output_dir, paste0(raster_base_name, "_adjusted_FDR-BH.tif")), overwrite = TRUE)
writeRaster(raster_adjusted_bky, file.path(output_dir, paste0(raster_base_name, "_adjusted_FDR-adaptive.tif")), overwrite = TRUE)
writeRaster(raster_binary_original, file.path(output_dir, paste0(raster_base_name, "_binary_original.tif")), overwrite = TRUE)
writeRaster(raster_binary_bh,  file.path(output_dir, paste0(raster_base_name, "_binary_FDR-BH.tif")), overwrite = TRUE)
writeRaster(raster_binary_bky, file.path(output_dir, paste0(raster_base_name, "_binary_FDR-adaptive.tif")), overwrite = TRUE)

cat("\n--- File Outputs ---\n")
cat("Original raster copy:", original_copy_path, "\n")
cat("Adjusted p-values raster (FDR-BH):", file.path(output_dir, paste0(raster_base_name, "_adjusted_FDR-BH.tif")), "\n")
cat("Adjusted p-values raster (Adaptive FDR):", file.path(output_dir, paste0(raster_base_name, "_adjusted_FDR-adaptive.tif")), "\n")
cat("Binary raster (raw p-values):", file.path(output_dir, paste0(raster_base_name, "_binary_original.tif")), "\n")
cat("Binary raster (FDR-BH):", file.path(output_dir, paste0(raster_base_name, "_binary_FDR-BH.tif")), "\n")
cat("Binary raster (Adaptive FDR):", file.path(output_dir, paste0(raster_base_name, "_binary_FDR-adaptive.tif")), "\n")
