#' Construct the matrix Q from graphon quantiles and covariates
#'
#' Constructs the matrix \(Q\) whose entries are given by differences of
#' \(F_v\) evaluated at graphon quantiles shifted by the inner products
#' \(X_j^\top a\). Specifically, for \(k = 1,\ldots,K\) and
#' \(j = 1,\ldots,n\),
#'
#' \deqn{
#' Q_{kj} =
#' F_v\left(\hat F_a^{-1}\left(\frac{k}{K}\right)
#' - X_j^\top a\right)
#' -
#' F_v\left(\hat F_a^{-1}\left(\frac{k-1}{K}\right)
#' - X_j^\top a\right).
#' }
#'
#' Here, \code{qgraphon} is a function that returns the graphon quantiles
#' \(\hat F_a^{-1}(u)\), \code{Fv} is the distribution function \(F_v\),
#' \code{a} is the parameter vector, and \code{matrix_X} contains the
#' covariate vectors \(X_j\) as rows.
#'
#' If \code{scaled = TRUE}, the resulting matrix is multiplied by
#' \(1/\sqrt{n}\).
#'
#' @param qgraphon A function that computes the graphon quantile function.
#'   It must accept a numeric vector of probabilities in \([0,1]\) and return
#'   the corresponding quantiles.
#' @param a A numeric parameter vector. Its length must equal the number of
#'   columns of \code{matrix_X}.
#' @param K A positive integer specifying the number of intervals used to
#'   construct the matrix \(Q\).
#' @param Fv A function representing the distribution function \(F_v\).
#'   It must accept numeric input and return values of the same length.
#' @param matrix_X A numeric matrix whose rows contain the covariate vectors
#'   \(X_j\). The number of columns must equal \code{length(a)}.
#' @param scaled Logical indicating whether the resulting matrix should be
#'   scaled by \(1/\sqrt{n}\), where \(n\) is the number of rows of
#'   \code{matrix_X}. Defaults to \code{FALSE}.
#'
#' @return A numeric \(K \times n\) matrix. The \((k,j)\)-th entry is
#'   \deqn{
#'   F_v\left(\hat F_a^{-1}(k/K) - X_j^\top a\right)
#'   -
#'   F_v\left(\hat F_a^{-1}((k-1)/K) - X_j^\top a\right).
#'   }
#'   If \code{scaled = TRUE}, the matrix is multiplied by \(1/\sqrt{n}\).
#'
#' @examples
#' n <- 100
#' K <- 3
#' a <- c(2.0, -0.5)
#'
#' X <- matrix(
#'   rnorm(2 * n),
#'   nrow = n,
#'   ncol = 2
#' )
#'
#' Fv <- function(x) {
#'   pnorm(x, mean = 0, sd = 1)
#' }
#'
#' qgraphon <- make_distribution_func(
#'   a = a,
#'   Fv = Fv,
#'   X_matrix = X
#' )
#'
#' Q <- create_matrix_Q(
#'   qgraphon = qgraphon,
#'   a = a,
#'   K = K,
#'   Fv = Fv,
#'   matrix_X = X
#' )
#'
#' dim(Q)
#'
#' @export
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 (!is.null(matrix_X) && 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