source(here::here("R", "graphon_distribution.R"))

compute_adj_matrix <- function(
  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)
  
  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)
  
  adj_mat <- (rand_mat < probabilities) * 1L
}

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)
adj