diff --git a/R/build_network.R b/R/build_network.R
index 772c7a1..6e606a0 100644
--- a/R/build_network.R
+++ b/R/build_network.R
@@ -1,31 +1,144 @@
-source(here::here("R", "graphon_distribution.R"))
+# 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
+    X_matrix,
+    v,
+    a,
+    phi,
+    rho_n,
+    Fv = pnorm
 ) {
-  xi <- vapply(v, FUN = function(x){pgraphon(x, a, Fv, X_matrix)}, numeric(1))
-  print(xi)
-  dim_adj <- length(v)
+  ## -------------------------- 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.")
+  }
   
-  rand_mat <- matrix(runif(dim_adj^2), dim_adj, dim_adj)
-  upper_mask <- upper.tri(rand_mat, diag=FALSE)
-  rand_mat[!upper_mask] <- 0
-  rand_mat <- rand_mat + t(rand_mat) + diag(nrow=dim_adj)
-  graphon_function <- function(x,y) { rho_n * phi(x,y)}
-  probabilities <- outer(xi, xi, FUN=graphon_function)
+  # 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(rnorm(4), nrow = 4, ncol=1)
-a <- 0.5
-v <- rnorm(4)
-
-adj <- compute_adj_matrix(X, v, a, phi = function(x,y) {x * y}, 0.5)
+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
+