create_matrix_Q <- function(
    qgraphon,
    a,
    K,
    Fv,
    matrix_X = NULL,
    scaled = FALSE
) {
  ## 1.1 Check inputs ==========================================================
  if (!is.numeric(a) || !is.vector(a)) stop("'a' must be a numeric vector")
  if (!is.numeric(K) || length(K) != 1 || K <= 0) stop("'K' must be a positive integer")

  if (!is.function(Fv)) stop("'F_v' must be a function")
  if (!is.matrix(matrix_X)) stop("matrix_X must be a matrix")
  if (!is.logical(scaled)) stop("`scaled` must be a logical!")
  if (ncol(matrix_X) != length(a)) {
    stop("Number of columns of `matrix_X` (", ncol(matrix_X), ") must equal length(a) (", length(a), ")")
  }
  
  ## 1.3 Compute the graphon quantiles =========================================
  k <- seq(0, K) / K
  n <- nrow(matrix_X)
  # here there is an automatic switch included, if fX is not null and we have a
  # scalar case, then qpgrahon automatically switches to the analytical 
  # expression. The intended use is for small values of n
  graphon_quantiles <- qgraphon(k)
  
  ## 1.4 Build the matrix Q ====================================================
  inner_products = as.vector(matrix_X %*% a)
  
  # outer(y, x, "-")  gives a matrix with entry (j,i) = y[j] - x[i]
  #   then we apply the CDF `F_v` to the whole matrix at once.
  #   finally we take the difference of successive rows (j) to obtain the
  #   increments required by equation (3.1).
  cdf_mat <- Fv(outer(graphon_quantiles, inner_products, "-")) # (K +1) x n matrix
  Q <- diff(cdf_mat, lag=1)                           # operates along rows
  
  if (scaled) { Q <- 1 / sqrt(n) * Q }
  Q
}

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

n <- 100
K <- 3
a <- c(2.0, -0.5)
X <- matrix(rnorm(2 * n), nrow = n, ncol = 2)
Fv <- function(x) {dnorm(x, mean=0, sd=1)}

qgraphon <- make_distribution_func(a=a, Fv=Fv, X_matrix=X)
Q <- create_matrix_Q(qgraphon, a, K, Fv, X)
Q