R/print_NBKP.R

In NBKP: Beta Kernel Process Modeling with Negative Binomial Response

#' @rdname print.NBKP
#' @title Print Methods for NBKP Objects
#'
#' @description Provides formatted console output for fitted NBKP model
#' objects and their predictions. The following specialized methods are supported:
#' \itemize{
#'   \item \code{print.NBKP} – display fitted NBKP model objects.
#'   \item \code{print.predict_NBKP} – display posterior predictive results.
#' }
#'
#' @param x An object of class \code{"NBKP"} or \code{"predict_NBKP"}.
#' @param ... Additional arguments passed to the generic \code{print} method
#'   (currently unused; included for S3 consistency).
#'
#' @return Invisibly returns the input object. Called for the side effect of
#'   printing human-readable summaries to the console.
#'
#' @seealso \code{\link{fit_NBKP}} for model fitting;
#'   \code{\link{predict.NBKP}} for posterior prediction.
#'
#' @references Zhao J, Qing K, Xu J (2025). \emph{BKP: An R Package for Beta
#'   Kernel Process Modeling}. arXiv. https://doi.org/10.48550/arxiv.2508.10447.
#'
#' @keywords NBKP
#'
#' @examples
#' \donttest{
#' # -------------------------- 1D Example --------------------------
#' set.seed(123)
#'
#' # Define true mean function
#' true_mu_fun <- function(x) {
#'   exp(sin(x) + 0.5)
#' }
#'
#' n <- 30
#' Xbounds <- matrix(c(-2, 2), nrow=1)
#' X <- tgp::lhs(n = n, rect = Xbounds)
#' true_mu <- true_mu_fun(X)
#' y <- rnbinom(n, size = 1, mu = true_mu)
#'
#' # Fit NBKP model
#' model1 <- fit_NBKP(X, y, Xbounds=Xbounds)
#' print(model1) # fitted object
#'
#' pred1 <- predict(model1)
#' print(pred1) # predictions
#'
#' # -------------------------- 2D Example --------------------------
#' set.seed(123)
#'
#' # Define 2D latent function and mean transformation
#' true_mu_fun <- function(X) {
#'   if(is.null(nrow(X))) X <- matrix(X, nrow=1)
#'   x1 <- 4*X[,1] - 2
#'   x2 <- 4*X[,2] - 2
#'   f <- sin(2*pi*x1) * cos(2*pi*x2)
#'   return(exp(f))
#' }
#'
#' n <- 100
#' Xbounds <- matrix(c(0, 0, 1, 1), nrow = 2)
#' X <- tgp::lhs(n = n, rect = Xbounds)
#' true_mu <- true_mu_fun(X)
#' y <- rnbinom(n, size = 0.5, mu = true_mu)
#'
#' # Fit NBKP model
#' model2 <- fit_NBKP(X, y, Xbounds=Xbounds)
#' print(model2)
#'
#' pred2 <- predict(model2)
#' print(pred2)
#' }
#'
#' @export
#' @method print NBKP
print.NBKP <- function(x, ...) {
  cat("\n      Negative Binomial Kernel Process (NBKP) Model    \n\n")
  cat(sprintf("Number of observations (n): %d\n", nrow(x$X)))
  cat(sprintf("Input dimensionality (d):   %d\n", ncol(x$X)))
  cat(sprintf("Kernel type:                %s\n", x$kernel))
  cat(sprintf("Optimized kernel parameters: %s\n",
              paste(sprintf("%.4f", x$theta_opt), collapse = ", ")))
  if (!is.na(x$loss_min)) {
    cat(sprintf("Minimum achieved loss:      %.5f\n", x$loss_min))
  }
  cat(sprintf("Loss function:              %s\n", x$loss))
  cat(sprintf("Prior type:                 %s\n", x$prior))
  if (x$prior == "fixed" || x$prior == "adaptive") {
    cat(sprintf("r0: %.3f\n", x$r0))
  }
  if (x$prior == "fixed") {
    cat(sprintf("mu0: %.3f\n", x$mu0))
  }
  cat(sprintf("Dispersion parameter (phi): %.3f\n", x$phi))
  invisible(x)
}

#' @rdname print.NBKP
#' @export
#' @method print predict_NBKP
print.predict_NBKP <- function(x, ...) {
  n <- length(x$mean)

  # Determine prediction input
  if (is.null(x$Xnew)) {
    cat("Prediction results on training data (X).\n")
    cat("Total number of training points:", n, "\n")
    X_disp <- x$X
  } else {
    cat("Prediction results on new data (Xnew).\n")
    cat("Total number of prediction points:", n, "\n")
    X_disp <- x$Xnew
  }

  d <- ncol(X_disp)

  # Determine how many rows to preview
  k <- min(6, n)
  if (n > k) {
    if (is.null(x$Xnew)) {
      cat("\nPreview of predictions for training data (first", k, "of", n, "points):\n")
    } else {
      cat("\nPreview of predictions for new data (first", k, "of", n, "points):\n")
    }
  } else {
    if (is.null(x$Xnew)) {
      cat("\nPredictions for all training data points:\n")
    } else {
      cat("\nPredictions for all new data points:\n")
    }
  }

  # Format X_disp for display
  X_preview <- head(X_disp, k)
  if (d == 1) {
    X_preview <- data.frame(x = round(X_preview, 4))
  } else if (d == 2) {
    X_preview <- as.data.frame(round(X_preview, 4))
    names(X_preview) <- c("x1", "x2")
  } else {
    # Only display first and last columns with ... in between
    X_preview_vals <- round(X_preview[, c(1, d), 4])
    X_preview <- as.data.frame(X_preview_vals)
    names(X_preview) <- c("x1", paste0("x", d))
    # Add a ... column
    X_preview$`...` <- rep("...", nrow(X_preview))
    # Reorder columns: x1, ..., xd
    X_preview <- X_preview[, c("x1", "...", paste0("x", d))]
  }

  # Construct results table
  pred_summary <- data.frame(
    mean = round(head(x$mean, k), 4),
    variance = round(head(x$variance, k), 4),
    lower = round(head(x$lower, k), 4),
    upper = round(head(x$upper, k), 4)
  )

  # Update CI column names
  ci_low <- round((1 - x$CI_level)/2 * 100, 2)
  ci_high <- round((1 + x$CI_level)/2 * 100, 2)
  names(pred_summary)[3:4] <- paste0(c(ci_low, ci_high), "% quantile")

  # Combine X preview and prediction
  res <- cbind(X_preview, pred_summary)
  print(res, row.names = FALSE)

  if (n > k) cat("...\n")
  invisible(x)
}