# Load required helper (the file must define `pgraphon`)
if (requireNamespace("here", quietly = TRUE)) {
  source(here::here("R", "graphon_distribution.R"))
} else {
  stop("Package `here` is required to locate the helper file.")
}

# TODO: Check the naming consistency across the documentation

#' Compute an adjacency matrix from a graphon model
#'
#' @title Generate a the adjacency matrix using a graphon specification
#' @description
#' Given a covariate matrix `X_matrix`, a vector of intercepts `v`, a
#' coefficient vector `a`, a symmetric graphon kernel `phi`, and a scaling
#' factor `rho_n`, this function draws a random undirected adjacency matrix.
#' The construction follows the steps:
#'   1. Compute latent variables \eqn{\xi_i = F_v(a^\top X_i + v_i)}.
#'   2. Generate a symmetric matrix of i.i.d. Uniform$[0,1]$ random numbers
#'      with ones on the diagonal.
#'   3. Evaluate the graphon probabilities \eqn{p_{ij}= \rho_n \, \phi(\xi_i,\xi_j)}.
#'   4. Set \eqn{A_{ij}=1} if the uniform draw is smaller than the probability,
#'      otherwise \eqn{0}.  The resulting matrix is symmetric with zeros on the
#'      diagonal.
#'
#' @param X_matrix A numeric matrix of dimension *n \times d* (covariates).  
#'   Must not contain `NA` or `NaN` values.
#' @param v A numeric vector of length *n* (intercepts).  
#'   Must be the same length as the number of rows of `X_matrix`.
#' @param a A numeric vector of length *d* (coefficients).  
#'   Must be the same length as the number of columns of `X_matrix`.
#' @param phi A binary function `phi(x, y)` returning a numeric value.  
#'   It should be symmetric (`phi(x, y) == phi(y, x)`) and bounded in $[0,1]$.
#' @param rho_n A numeric scalar in $[0,1]$ that scales the graphon.
#' @param Fv A cumulative distribution function (CDF) used to transform the
#'   linear predictor.  Defaults to the standard normal CDF `pnorm`.  
#'   Must accept a numeric vector and return a numeric vector of the same length.
#'
#' @return An *n* × *n* integer matrix (`0L`/`1L`) representing the adjacency
#'   matrix of an undirected graph.
#'
#' @examples
#' ## Simple example with a multiplicative graphon
#' set.seed(1)
#' X  <- matrix(seq(-1, 1, length.out = 5), ncol = 1)
#' a  <- 2
#' v  <- seq(0, 0.8, length.out = 5)
#' phi_fun <- function(x, y) x * y          # multiplicative kernel
#' adj <- compute_adj_matrix(
#'   X_matrix = X,
#'   v        = v,
#'   a        = a,
#'   phi      = phi_fun,
#'   rho_n    = 0.5,
#'   Fv       = punif
#' )
#' adj
#'
#' @export
compute_adj_matrix <- function(
    X_matrix,
    v,
    a,
    phi,
    rho_n,
    Fv = pnorm
) {
  ## -------------------------- Guard rails -------------------------- ##
  # Missing / NULL checks
  if (missing(X_matrix) || is.null(X_matrix)) {
    stop("`X_matrix` is missing or NULL.")
  }
  if (missing(v) || is.null(v)) {
    stop("`v` is missing or NULL.")
  }
  if (missing(a) || is.null(a)) {
    stop("`a` is missing or NULL.")
  }
  if (missing(phi) || is.null(phi)) {
    stop("`phi` is missing or NULL.")
  }
  if (missing(rho_n) || is.null(rho_n)) {
    stop("`rho_n` is missing or NULL.")
  }
  if (missing(Fv) || is.null(Fv)) {
    stop("`Fv` is missing or NULL.")
  }
  
  # Type / dimension checks
  if (!is.matrix(X_matrix) || !is.numeric(X_matrix)) {
    stop("`X_matrix` must be a numeric matrix.")
  }
  if (any(is.na(X_matrix) | is.nan(X_matrix))) {
    stop("`X_matrix` must not contain NA or NaN values.")
  }
  n <- nrow(X_matrix)
  d <- ncol(X_matrix)
  
  if (!is.numeric(v) || length(v) != n) {
    stop("`v` must be a numeric vector of length equal to nrow(`X_matrix`).")
  }
  if (!is.numeric(a) || length(a) != d) {
    stop("`a` must be a numeric vector of length equal to ncol(`X_matrix`).")
  }
  if (!is.function(phi)) {
    stop("`phi` must be a function of two numeric arguments.")
  }
  if (!is.numeric(rho_n) || length(rho_n) != 1 || rho_n < 0 || rho_n > 1) {
    stop("`rho_n` must be a numeric scalar in [0, 1].")
  }
  if (!is.function(Fv)) {
    stop("`Fv` must be a function (CDF) that accepts a numeric vector.")
  }
  
  ## -------------------------- Core computation -------------------------- ##
  # 1. Latent variables ξ_i = Fv(aᵀ X_i + v_i)
  lin_pred <- as.vector(X_matrix %*% a) + v
  xi <- pgraphon(y = lin_pred, a = a, Fv = Fv, X_matrix = X_matrix)
  
  # 2. Symmetric matrix of Uniform[0,1] draws with 1 on the diagonal
  #    (only the upper triangle is sampled, then mirrored)
  rand_mat <- matrix(0, n, n)
  upper_idx <- which(upper.tri(rand_mat, diag = FALSE), arr.ind = TRUE)
  rand_mat[upper_idx] <- runif(n * (n - 1) / 2)
  rand_mat <- rand_mat + t(rand_mat) + diag(1, n)
  
  # 3. Graphon probabilities p_ij = rho_n * phi(ξ_i, ξ_j)
  graphon_fun <- function(x, y) rho_n * phi(x, y)
  probabilities <- outer(xi, xi, FUN = graphon_fun)
  
  # 4. Adjacency matrix: A_ij = 1 if U_ij < p_ij, else 0
  adj_mat <- (rand_mat < probabilities) * 1L
  
  ## Return the adjacency matrix
  adj_mat
}

set.seed(1)
X <- matrix(c(-1, -0.5, 0, 0.5, 1), nrow = 5, ncol=1)
a <- 2.0
v <- c(0, 0.2, 0.4, 0.6, 0.8)
adj <- compute_adj_matrix(X, v, a, phi = function(x,y) {x * y}, 0.5, Fv = punif)
adj