R/TOwen.R

In snreg: Regression with Skew-Normally Distributed Error Term

#' Owen's T Function via C Backend
#'
#' @title Compute Owen's T Function \code{T(h, a)}
#'
#' @description
#' \code{TOwen1} computes an Owen's \eqn{T}-function variant (or a related
#' special function) for vectors \code{h} and \code{a} based on the \code{t} function in \url{https://people.sc.fsu.edu/~jburkardt/c_src/owen/owen.html}. Non-finite inputs (in \code{h} or \code{a}) produce \code{NA}
#' at corresponding positions, while finite pairs are computed in C in a
#' vectorized fashion.
#'
#' @details
#' Owen's \eqn{T} function is commonly defined as
#' \deqn{T(h, a) \;=\; \frac{1}{2\pi} \int_{0}^{a} \frac{\exp\!\left(-\tfrac{1}{2}h^2 (1+t^2)\right)}{1+t^2} \, dt,}
#' for real \eqn{h} and \eqn{a}.
#'
#' The function accepts vector inputs and:
#' \itemize{
#'   \item Computes results only for entries where both \code{h} and \code{a} are finite.
#'   \item Returns \code{NA} where either \code{h} or \code{a} is non-finite.
#'   \item Optionally passes a \code{threads} hint to the C backend (ignored if not supported).
#' }
#'
#' @param h
#' numeric vector of \eqn{h} arguments.
#'
#' @param a
#' numeric vector of \eqn{a} arguments. Must be either the same length as \code{h}
#' or of length 1 (will be recycled by standard R rules).
#'
#' @param threads
#' integer. Number of threads to request from the C implementation (if supported).
#' Default is \code{1}.
#'
#' @return
#' A numeric vector of length \code{length(h)} containing \eqn{T(h_i, a_i)}.
#' Elements where either \code{h_i} or \code{a_i} is not finite are \code{NA}.
#' The returned object is given class \code{"snreg"} for downstream compatibility
#' with your package’s print/summary helpers.
#'
#' @examples
#' library(snreg)
#' 
#' # Basic usage. Vectorized 'a'
#' h <- c(-1, 0, 1, 2)
#' a <- 0.5
#' TOwen(h, a)
#' 
#' # Vectorized 'a' with non-finite entries; non-finite entries yield NA
#' a2 <- c(0.2, NA, 1, Inf)
#' TOwen(h, a2)
#'
#' @seealso
#' \code{\link{pnorm}}, \code{\link{dnorm}}
#'
#' @export
TOwen <- function(h, a, threads = 1) {
  n <- length(h)
  
  # Recycle 'a' if needed (standard R recycling rules)
  if (length(a) != n) {
    a <- rep(a, length.out = n)
  }
  
  # Indices to compute: both finite
  toInclude <- is.finite(h) & is.finite(a)
  
  if (sum(toInclude) < n) {
    out <- rep(NA_real_, n)
    n1  <- sum(toInclude)
    if (n1 > 0) {
      h1 <- as.double(h[toInclude])
      a1 <- as.double(a[toInclude])
      out[toInclude] <- .C(
        "TOwen",
        as.integer(n1),
        h1,
        a1,
        TOwenValues = double(n1),
        as.integer(threads)
      )$TOwenValues
    }
  } else {
    out <- .C(
      "TOwen",
      as.integer(n),
      as.double(h),
      as.double(a),
      TOwenValues = double(n),
      as.integer(threads)
    )$TOwenValues
  }
  
  class(out) <- "snreg"
  out
}

# OLD before AI

# TOwen <- function(h,a,threads = 1){
#   n <- length(h)
#   # cat.print(h)
#   # cat.print(a)
#   # cat.print(threads)
#   toInclude <- is.finite(h) & is.finite(a) 
#   # toInclude <- rep(TRUE,n)
#   if(sum(toInclude) < n){
#     tymch1 <- rep(NA, n)
#     n1 <- sum(toInclude)
#     h1 <- h[toInclude]
#     a1 <- a[toInclude]
#     tymch1[toInclude] <- .C("TOwen", as.integer(n1), as.double(h1), as.double(a1), TOwenValues = double(n1), as.integer(threads)) $ TOwenValues
#   } else {
#     tymch1 <- .C("TOwen", as.integer(n), as.double(h), as.double(a), TOwenValues = double(n), as.integer(threads)) $ TOwenValues
#   }
#   class(tymch1) <- "snreg"
#   return(tymch1)
# }