R/np.qregression.R

In np: Nonparametric Kernel Smoothing Methods for Mixed Data Types

npqreg <-
  function(bws, ...){
    args <- list(...)

    if (!missing(bws)){
      if (is.recursive(bws)){
        if (!is.null(bws$formula) && is.null(args$txdat))
          UseMethod("npqreg",bws$formula)
        else if (!is.null(bws$call) && is.null(args$txdat))
          UseMethod("npqreg",bws$call)
        else if (!is.call(bws))
          UseMethod("npqreg",bws)
        else
          UseMethod("npqreg",NULL)
      } else {
        UseMethod("npqreg", NULL)
      }
    } else {
      UseMethod("npqreg", NULL)
    }
  }

.npqreg.fit.control.names <- c("data", "newdata", "exdat", "tau", "gradients", "tol", "small", "itmax")
.npqreg.removed.solver.controls <- c("ftol",
                                     "lbc.dir", "dfc.dir", "cfac.dir", "initc.dir",
                                     "lbd.dir", "hbd.dir", "dfac.dir", "initd.dir")

.npqreg_validate_tau <- function(tau) {
  if (!is.numeric(tau) || !length(tau) || anyNA(tau) ||
      any(!is.finite(tau)) || any(tau <= 0) || any(tau >= 1))
    stop("'tau' must contain numeric values in (0,1)")
  as.double(tau)
}

.npqreg_tau_labels <- function(tau) {
  paste0("tau=", format(tau, trim = TRUE, scientific = FALSE))
}

.npqreg_napredict_eval <- function(omit, x) {
  if (!length(omit))
    return(x)
  omit <- as.integer(omit)
  keep <- seq_len(NROW(x) + length(omit))[-omit]
  if (is.null(dim(x))) {
    out <- x[rep(NA_integer_, length(x) + length(omit))]
    out[keep] <- x
    attr(out, "na.action") <- NULL
    return(out)
  }
  if (is.data.frame(x)) {
    out <- x[rep(NA_integer_, nrow(x) + length(omit)), , drop = FALSE]
    out[keep, ] <- x
    attr(out, "na.action") <- NULL
    return(out)
  }
  if (length(dim(x)) <= 2L) {
    out <- x[rep(NA_integer_, nrow(x) + length(omit)), , drop = FALSE]
    out[keep, ] <- x
    attr(out, "na.action") <- NULL
    return(out)
  }
  d <- dim(x)
  dn <- dimnames(x)
  new.dim <- c(d[1L] + length(omit), d[-1L])
  if (!is.null(dn)) {
    dn[[1L]] <- NULL
    if (length(dn) != length(new.dim) ||
        any(vapply(seq_along(dn), function(i) {
          !is.null(dn[[i]]) && length(dn[[i]]) != new.dim[[i]]
        }, logical(1L))))
      dn <- NULL
  }
  out <- array(NA_real_,
               dim = new.dim,
               dimnames = dn)
  out[keep, , ] <- x
  out
}

.npqreg_validate_newdata_terms <- function(newdata, xnames) {
  nd <- toFrame(newdata)
  missing.names <- setdiff(xnames, names(nd))
  if (length(missing.names))
    stop(sprintf(
      "newdata must contain columns: %s",
      paste(shQuote(xnames), collapse = ", ")
    ), call. = FALSE)
  invisible(TRUE)
}

.npqreg_fit_dots <- function(dots, allow.bandwidth.controls = FALSE) {
  dot.names <- names(dots)
  if (is.null(dot.names))
    dot.names <- rep("", length(dots))

  stale <- intersect(dot.names[nzchar(dot.names)], .npqreg.removed.solver.controls)
  if (length(stale) && !allow.bandwidth.controls) {
    stop(sprintf(
      "'%s' %s no longer accepted by npqreg; the canonical one-dimensional quantile extractor is controlled by 'tol', 'small', and 'itmax'",
      paste(stale, collapse = "', '"),
      if (length(stale) == 1L) "is" else "are"
    ))
  }

  if (!allow.bandwidth.controls) {
    bad <- dot.names == "" | !(dot.names %in% .npqreg.fit.control.names)
    if (any(bad))
      .np_reject_unused_dots(dots[bad], "npqreg")
  }

  keep <- (!nzchar(dot.names)) | (dot.names %in% .npqreg.fit.control.names)
  dots[keep]
}

.npqreg_reject_gradient_order_dots <- function(dots) {
  dot.names <- names(dots)
  if (is.null(dot.names))
    dot.names <- rep("", length(dots))

  bad <- dot.names %in% c("gradient.order", "gradient_order")
  if (any(bad))
    .np_reject_unused_dots(dots[which(bad)[1L]], "npqreg")

  invisible(TRUE)
}

.npqreg_validate_itmax <- function(itmax) {
  if (!is.numeric(itmax) || length(itmax) != 1L || is.na(itmax) ||
      !is.finite(itmax) || itmax < 1 || itmax != floor(itmax) ||
      itmax > .Machine$integer.max)
    stop("'itmax' must be a positive integer")
  as.integer(itmax)
}

.npqreg_quantile_clamp_attr <- "npqreg.clamp"

.npqreg_quantile_clamp <- function(quantile) {
  clamp <- attr(quantile, .npqreg_quantile_clamp_attr, exact = TRUE)
  if (is.null(clamp) || length(clamp) != length(quantile))
    return(rep.int("none", length(quantile)))
  as.character(clamp)
}

.npqreg_mark_clamped_delta <- function(delta, clamp) {
  clamp <- as.character(clamp)
  bad <- which(!is.na(clamp) & clamp != "none")
  if (!length(bad))
    return(delta)

  if (!is.null(delta$quanterr) && length(delta$quanterr) >= max(bad))
    delta$quanterr[bad] <- NA_real_
  if (is.matrix(delta$quantgrad) && nrow(delta$quantgrad) >= max(bad))
    delta$quantgrad[bad, ] <- NA_real_
  if (is.matrix(delta$quantgerr) && nrow(delta$quantgerr) >= max(bad))
    delta$quantgerr[bad, ] <- NA_real_
  delta
}

.npqreg_assert_monotone_cdf_kernel <- function(bws) {
  order <- bws$cykerorder
  if (is.null(order))
    return(invisible(TRUE))
  order <- suppressWarnings(as.integer(order[1L]))
  if (!is.finite(order) || order != 2L)
    stop("npqreg requires cykerorder = 2 for conditional-quantile inversion; higher-order dependent-variable CDF kernels can be nonmonotone",
         call. = FALSE)
  invisible(TRUE)
}

.npqreg_strip_fit_controls_from_bw_call <- function(call) {
  for (nm in c("tau", "gradients", "tol", "small", "itmax", "newdata", "exdat")) {
    if (nm %in% names(call))
      call[[nm]] <- NULL
  }
  call
}

.npqreg_quantile_delta_from_conditional <- function(bws,
                                                    xdat,
                                                    ydat,
                                                    exdat,
                                                    quantile,
                                                    gradients = FALSE) {
  xdat <- toFrame(xdat)
  ydat <- toFrame(ydat)
  exdat <- toFrame(exdat)
  quantile <- as.double(quantile)
  gradients <- npValidateScalarLogical(gradients, "gradients")

  if (length(quantile) != nrow(exdat))
    stop("quantile delta helper requires one quantile per evaluation row")
  if (ncol(ydat) != 1L)
    stop("quantile delta helper requires a single response")

  eydat <- stats::setNames(data.frame(quantile), names(ydat)[1L])
  cdf.obj <- .np_conditional_eval_selected(
    bws = bws,
    xdat = xdat,
    ydat = ydat,
    exdat = exdat,
    eydat = eydat,
    cdf = TRUE,
    gradients = gradients
  )
  dens.obj <- .np_conditional_eval_selected(
    bws = bws,
    xdat = xdat,
    ydat = ydat,
    exdat = exdat,
    eydat = eydat,
    cdf = FALSE,
    gradients = FALSE
  )

  dens <- as.double(dens.obj$condens)
  quanterr <- as.double(cdf.obj$conderr) / NZD(dens)
  quanterr[!is.finite(quanterr) | quanterr < 0.0] <- NA_real_

  if (!gradients) {
    return(list(
      quanterr = quanterr,
      quantgrad = NA,
      quantgerr = NA,
      cdf = cdf.obj,
      dens = dens.obj
    ))
  }

  dens.mat <- matrix(NZD(dens),
                     nrow = nrow(cdf.obj$congrad),
                     ncol = ncol(cdf.obj$congrad))
  grad <- -cdf.obj$congrad / dens.mat
  grad[!is.finite(grad)] <- NA_real_

  gerr <- cdf.obj$congerr / dens.mat
  gerr[!is.finite(gerr) | gerr < 0.0] <- NA_real_

  list(
    quanterr = quanterr,
    quantgrad = grad,
    quantgerr = gerr,
    cdf = cdf.obj,
    dens = dens.obj
  )
}

.npqreg_selected_cdf_values <- function(bws,
                                        xdat,
                                        ydat,
                                        exdat,
                                        ycand) {
  ydat <- toFrame(ydat)
  yname <- names(ydat)[1L]
  eydat <- stats::setNames(data.frame(as.double(ycand)), yname)
  as.double(.np_conditional_eval_selected(
    bws = bws,
    xdat = xdat,
    ydat = ydat,
    exdat = exdat,
    eydat = eydat,
    cdf = TRUE,
    gradients = FALSE
  )$condist)
}

.npqreg_selected_cdf_cache_atom <- function(x) {
  if (inherits(x, "factor"))
    return(paste0("f:", as.integer(x)))
  if (inherits(x, "Date") || inherits(x, "POSIXt"))
    return(paste0("t:", sprintf("%a", as.double(x))))
  if (is.numeric(x) || is.integer(x) || is.logical(x))
    return(paste0("n:", sprintf("%a", as.double(x))))
  paste0("c:", as.character(x))
}

.npqreg_selected_cdf_cache_row_key <- function(exdat, i) {
  row <- exdat[i, , drop = FALSE]
  parts <- vapply(row, function(x) .npqreg_selected_cdf_cache_atom(x[[1L]]), character(1L))
  paste(parts, collapse = "\r")
}

.npqreg_selected_cdf_cache_row_keys <- function(exdat) {
  exdat <- as.data.frame(exdat)
  if (!nrow(exdat))
    return(character(0L))
  vapply(seq_len(nrow(exdat)), function(i) .npqreg_selected_cdf_cache_row_key(exdat, i), character(1L))
}

.npqreg_selected_cdf_cache_key <- function(row.key, ycand) {
  paste(c(row.key, paste0("y:", sprintf("%a", as.double(ycand)))), collapse = "\r")
}

.npqreg_selected_cdf_cache_new <- function(enabled) {
  cache <- new.env(parent = emptyenv(), hash = FALSE)
  cache$enabled <- isTRUE(enabled)
  cache$store <- new.env(parent = emptyenv(), hash = TRUE)
  cache$visits <- 0L
  cache$hits <- 0L
  cache$misses <- 0L
  cache$unique <- 0L
  cache
}

.npqreg_selected_cdf_cache_clear <- function(cache) {
  if (is.environment(cache)) {
    cache$enabled <- FALSE
    cache$store <- new.env(parent = emptyenv(), hash = TRUE)
  }
  invisible(NULL)
}

.npqreg_selected_cdf_cache_should_enable <- function(tau, exdat) {
  isTRUE(npObjectiveCacheEnabled()) &&
    (length(tau) > 1L || anyDuplicated(as.data.frame(exdat)) > 0L)
}

.npqreg_selected_cdf_values_cached <- function(bws,
                                               xdat,
                                               ydat,
                                               exdat,
                                               ycand,
                                               cdf.cache = NULL,
                                               row.keys = NULL,
                                               cdf.values = .npqreg_selected_cdf_values) {
  if (!is.environment(cdf.cache) || !isTRUE(cdf.cache$enabled))
    return(cdf.values(bws, xdat, ydat, exdat, ycand))

  exdat <- toFrame(exdat)
  n.eval <- nrow(exdat)
  ycand <- as.double(ycand)
  if (length(ycand) == 1L && n.eval > 1L)
    ycand <- rep.int(ycand, n.eval)
  if (length(ycand) != n.eval)
    return(cdf.values(bws, xdat, ydat, exdat, ycand))
  if (is.null(row.keys) || length(row.keys) != n.eval)
    row.keys <- .npqreg_selected_cdf_cache_row_keys(exdat)

  out <- rep.int(NA_real_, n.eval)
  keys <- character(n.eval)
  miss.first <- integer(0L)
  miss.keys <- character(0L)
  miss.map <- new.env(parent = emptyenv(), hash = TRUE)

  for (i in seq_len(n.eval)) {
    key <- .npqreg_selected_cdf_cache_key(row.keys[[i]], ycand[[i]])
    keys[[i]] <- key
    cdf.cache$visits <- cdf.cache$visits + 1L
    if (exists(key, envir = cdf.cache$store, inherits = FALSE)) {
      cdf.cache$hits <- cdf.cache$hits + 1L
      out[[i]] <- get(key, envir = cdf.cache$store, inherits = FALSE)
    } else {
      cdf.cache$misses <- cdf.cache$misses + 1L
      if (!exists(key, envir = miss.map, inherits = FALSE)) {
        assign(key, length(miss.first) + 1L, envir = miss.map)
        miss.first <- c(miss.first, i)
        miss.keys <- c(miss.keys, key)
      }
    }
  }

  if (length(miss.first)) {
    values <- cdf.values(
      bws = bws,
      xdat = xdat,
      ydat = ydat,
      exdat = exdat[miss.first, , drop = FALSE],
      ycand = ycand[miss.first]
    )
    for (j in seq_along(miss.first)) {
      assign(miss.keys[[j]], values[[j]], envir = cdf.cache$store)
      cdf.cache$unique <- cdf.cache$unique + 1L
    }
  }

  missing.out <- which(!is.finite(out) & is.na(out))
  if (length(missing.out)) {
    for (i in missing.out)
      out[[i]] <- get(keys[[i]], envir = cdf.cache$store, inherits = FALSE)
  }
  out
}

.npqreg_assert_selected_cdf_metadata <- function(bws) {
  .npqreg_assert_monotone_cdf_kernel(bws)
  reg.engine <- if (is.null(bws$regtype.engine)) {
    if (is.null(bws$regtype)) "lc" else as.character(bws$regtype)
  } else {
    as.character(bws$regtype.engine)
  }
  if (!identical(reg.engine, "lp"))
    return(invisible(TRUE))

  if (is.null(bws$degree.engine) && is.null(bws$degree))
    stop("selected LP conditional distribution metadata missing from bandwidth object: degree")
  invisible(TRUE)
}

.npqreg_invert_selected_cdf <- function(bws,
                                        xdat,
                                        ydat,
                                        exdat,
                                        tau,
                                        tol,
                                        small,
                                        itmax,
                                        cdf.cache = NULL,
                                        cdf.row.keys = NULL,
                                        cdf.values = .npqreg_selected_cdf_values) {
  .npqreg_assert_selected_cdf_metadata(bws)

  xdat <- toFrame(xdat)
  ydat <- toFrame(ydat)
  exdat <- toFrame(exdat)
  itmax <- .npqreg_validate_itmax(itmax)
  y <- as.double(ydat[[1L]])
  y <- y[is.finite(y)]
  if (!length(y))
    stop("npqreg selected-CDF inversion requires finite response support")

  n.eval <- nrow(exdat)
  y.min <- min(y)
  y.max <- max(y)
  if (!is.finite(y.min) || !is.finite(y.max))
    stop("npqreg selected-CDF inversion found non-finite response support")
  if (identical(y.min, y.max)) {
    out <- rep.int(y.min, n.eval)
    attr(out, .npqreg_quantile_clamp_attr) <- rep.int("constant", n.eval)
    return(out)
  }

  lo <- rep.int(y.min, n.eval)
  hi <- rep.int(y.max, n.eval)

  flo <- .npqreg_selected_cdf_values_cached(bws, xdat, ydat, exdat, lo, cdf.cache = cdf.cache,
                                            row.keys = cdf.row.keys, cdf.values = cdf.values)
  fhi <- .npqreg_selected_cdf_values_cached(bws, xdat, ydat, exdat, hi, cdf.cache = cdf.cache,
                                            row.keys = cdf.row.keys, cdf.values = cdf.values)
  if (any(!is.finite(flo)) || any(!is.finite(fhi)))
    stop("npqreg selected-CDF inversion encountered non-finite bracket values")

  done.low <- flo >= tau
  done.high <- fhi < tau
  active <- !(done.low | done.high)

  maxiter <- itmax
  iter <- 0L
  while (any(active) && iter < maxiter) {
    iter <- iter + 1L
    mid <- (lo[active] + hi[active]) / 2.0
    fmid <- .npqreg_selected_cdf_values_cached(
      bws = bws,
      xdat = xdat,
      ydat = ydat,
      exdat = exdat[active, , drop = FALSE],
      ycand = mid,
      cdf.cache = cdf.cache,
      row.keys = if (is.null(cdf.row.keys)) NULL else cdf.row.keys[active],
      cdf.values = cdf.values
    )
    if (any(!is.finite(fmid)))
      stop("npqreg selected-CDF inversion encountered non-finite refinement values")

    active.idx <- which(active)
    upper <- fmid >= tau
    hi[active.idx[upper]] <- mid[upper]
    lo[active.idx[!upper]] <- mid[!upper]

    width <- hi[active.idx] - lo[active.idx]
    scale <- pmax(abs(hi[active.idx]), abs(lo[active.idx]), 1.0)
    active[active.idx] <- width > (tol * scale + small)
  }

  if (any(active))
    stop("npqreg selected-CDF inversion failed to converge within 'itmax'")

  out <- (lo + hi) / 2.0
  out[done.low] <- y.min
  out[done.high] <- y.max
  clamp <- rep.int("none", n.eval)
  clamp[done.low] <- "lower"
  clamp[done.high] <- "upper"
  attr(out, .npqreg_quantile_clamp_attr) <- clamp
  out
}

npqreg.formula <-
  function(bws, data = NULL, newdata = NULL, ...){

    tt <- terms(bws)
    m <- match(c("formula", "data", "subset", "na.action"),
               names(bws$call), nomatch = 0)
    tmf <- bws$call[c(1,m)]
    tmf[[1]] <- as.name("model.frame")
    tmf[["formula"]] <- tt
    if (!is.null(data))
      tmf[["data"]] <- substitute(data)
    mf.args <- as.list(tmf)[-1L]
    umf <- tmf <- do.call(stats::model.frame, mf.args, envir = environment(tt))

    tydat <- tmf[, bws$variableNames[["response"]], drop = FALSE]
    txdat <- tmf[, bws$variableNames[["terms"]], drop = FALSE]

    has.eval <- !is.null(newdata)
    if (has.eval) {
      .npqreg_validate_newdata_terms(newdata, bws$variableNames[["terms"]])
      tt <- drop.terms(tt, match(bws$variableNames$response, attr(tt, 'term.labels')))
      umf.args <- list(formula = tt, data = newdata)
      umf <- do.call(stats::model.frame, umf.args, envir = parent.frame())
      emf <- umf
      exdat <- emf[, bws$variableNames[["terms"]], drop = FALSE]
    }

    q.args <- list(txdat = txdat, tydat = tydat)
    if (has.eval)
      q.args$exdat <- exdat
    q.args$bws <- bws
    tbw <- do.call(npqreg, c(q.args, .npqreg_fit_dots(list(...))))

    tbw$omit <- attr(umf,"na.action")
    tbw$rows.omit <- as.vector(tbw$omit)
    tbw$nobs.omit <- length(tbw$rows.omit)

    tbw$quantile <- .npqreg_napredict_eval(tbw$omit, tbw$quantile)
    tbw$quanterr <- .npqreg_napredict_eval(tbw$omit, tbw$quanterr)

    if(tbw$gradients){
        tbw$quantgrad <- .npqreg_napredict_eval(tbw$omit, tbw$quantgrad)
        tbw$quantgerr <- .npqreg_napredict_eval(tbw$omit, tbw$quantgerr)
    }

    return(tbw)
  }

npqreg.call <-
  function(bws, ...) {
    npqreg(txdat = .np_eval_bws_call_arg(bws, "xdat"),
           tydat = .np_eval_bws_call_arg(bws, "ydat"),
           bws = bws, ...)
  }

npqreg.conbandwidth <-
  function(bws, ...){
    stop("incorrect bandwidth type: expected conditional distribution bandwidths instead of conditional density bandwidths")
  }

npqreg.condbandwidth <-
  function(bws,
           txdat = stop("training data 'txdat' missing"),
           tydat = stop("training data 'tydat' missing"),
           exdat,
           tau = 0.5,
           gradients = FALSE,
           tol = 1.490116e-04,
           small = 1.490116e-05, itmax = 10000,
           ...){

    fit.start <- proc.time()[3]
    fit.dots <- .npqreg_fit_dots(list(...))
    if (length(fit.dots))
      stop(sprintf("unused npqreg fit argument '%s'", names(fit.dots)[1L]))
    gradients <- npValidateScalarLogical(gradients, "gradients")
    itmax <- .npqreg_validate_itmax(itmax)
    if (!is.numeric(tol) || length(tol) != 1L || is.na(tol) ||
        !is.finite(tol) || tol <= 0)
      stop("'tol' must be a positive finite numeric scalar")
    if (!is.numeric(small) || length(small) != 1L || is.na(small) ||
        !is.finite(small) || small <= 0)
      stop("'small' must be a positive finite numeric scalar")
    tol <- as.double(tol)
    small <- as.double(small)

    no.ex = missing(exdat)

    txdat = toFrame(txdat)
    tydat = toFrame(tydat)

    tau <- .npqreg_validate_tau(tau)
    ntau <- length(tau)
    tau.labels <- .npqreg_tau_labels(tau)

    if (dim(tydat)[2] != 1)
      stop("'tydat' has more than one column")

    if (!no.ex){
      exdat = toFrame(exdat)
      
      if (! txdat %~% exdat )
        stop("'txdat' and 'exdat' are not similar data frames!")
    }

    if (length(bws$xbw) != length(txdat))
      stop("length of bandwidth vector does not match number of columns of 'txdat'")

    if (length(bws$ybw) != 1)
      stop("length of bandwidth vector does not match number of columns of 'tydat'")

    if (any(bws$iyord) || any(bws$iyuno) || coarseclass(tydat[,1]) != "numeric")
      stop("'tydat' is not continuous")

    if ((any(bws$ixcon) &&
         !all(vapply(txdat[, bws$ixcon, drop = FALSE], inherits, logical(1), c("integer", "numeric")))) ||
        (any(bws$ixord) &&
         !all(vapply(txdat[, bws$ixord, drop = FALSE], inherits, logical(1), "ordered"))) ||
        (any(bws$ixuno) &&
         !all(vapply(txdat[, bws$ixuno, drop = FALSE], inherits, logical(1), "factor"))))
      stop("supplied bandwidths do not match 'txdat' in type")

    ## catch and destroy NA's
    keep.rows <- rep_len(TRUE, nrow(txdat))
    rows.omit <- attr(na.omit(data.frame(txdat, tydat)), "na.action")
    if (length(rows.omit) > 0L)
      keep.rows[as.integer(rows.omit)] <- FALSE

    if (!any(keep.rows))
      stop("Training data has no rows without NAs")

    txdat <- txdat[keep.rows,,drop = FALSE]
    tydat <- tydat[keep.rows,,drop = FALSE]

    eval.omit <- NULL
    if (!no.ex){
      keep.eval <- rep_len(TRUE, nrow(exdat))
      eval.omit <- attr(na.omit(exdat), "na.action")
      if (length(eval.omit) > 0L)
        keep.eval[as.integer(eval.omit)] <- FALSE
      exdat <- exdat[keep.eval,,drop = FALSE]
    }
    
    tnrow = dim(txdat)[1]
    enrow = (if (no.ex) tnrow else dim(exdat)[1])

    ## re-assign levels in training and evaluation data to ensure correct
    ## conversion to numeric type.
    
    txdat <- adjustLevels(txdat, bws$xdati)
    tydat <- adjustLevels(tydat, bws$ydati)
    
    if (!no.ex){
      exdat <- adjustLevels(exdat, bws$xdati, allowNewCells = TRUE)
    }

    ## grab the evaluation data before it is converted to numeric
    if(no.ex){
      txeval <- txdat
    } else {
      txeval <- exdat
    }
    txdat.df <- txdat
    tydat.df <- tydat
    if (!no.ex)
      exdat.df <- exdat

    cdf.cache <- .npqreg_selected_cdf_cache_new(
      .npqreg_selected_cdf_cache_should_enable(tau = tau, exdat = txeval)
    )
    cdf.row.keys <- if (isTRUE(cdf.cache$enabled))
      .npqreg_selected_cdf_cache_row_keys(txeval)
    else NULL
    on.exit(.npqreg_selected_cdf_cache_clear(cdf.cache), add = TRUE)

    fit_one_tau <- function(tau_i) {
      yq <- .npqreg_invert_selected_cdf(
        bws = bws,
        xdat = txdat.df,
        ydat = tydat.df,
        exdat = txeval,
        tau = tau_i,
        tol = tol,
        small = small,
        itmax = itmax,
        cdf.cache = cdf.cache,
        cdf.row.keys = cdf.row.keys
      )
      qclamp <- .npqreg_quantile_clamp(yq)
      qdelta <- .npqreg_quantile_delta_from_conditional(
        bws = bws,
        xdat = txdat.df,
        ydat = tydat.df,
        exdat = txeval,
        quantile = yq,
        gradients = gradients
      )
      qdelta <- .npqreg_mark_clamped_delta(qdelta, qclamp)
      list(
        yq = yq,
        yqerr = qdelta$quanterr,
        yqgrad = if (gradients) qdelta$quantgrad else NA,
        yqgerr = if (gradients) qdelta$quantgerr else NA
      )
    }

    tau.out <- lapply(tau, fit_one_tau)

    if (ntau == 1L) {
      myout <- tau.out[[1L]]
    } else {
      myout <- list(
        yq = do.call(cbind, lapply(tau.out, `[[`, "yq")),
        yqerr = do.call(cbind, lapply(tau.out, `[[`, "yqerr")),
        yqgrad = NA,
        yqgerr = NA
      )
      colnames(myout$yq) <- tau.labels
      colnames(myout$yqerr) <- tau.labels
      if (gradients) {
        p <- ncol(tau.out[[1L]]$yqgrad)
        grad.names <- colnames(tau.out[[1L]]$yqgrad)
        myout$yqgrad <- array(NA_real_,
                              dim = c(enrow, p, ntau),
                              dimnames = list(NULL, grad.names, tau.labels))
        myout$yqgerr <- array(NA_real_,
                              dim = c(enrow, p, ntau),
                              dimnames = list(NULL, grad.names, tau.labels))
        for (j in seq_len(ntau)) {
          myout$yqgrad[, , j] <- tau.out[[j]]$yqgrad
          myout$yqgerr[, , j] <- tau.out[[j]]$yqgerr
        }
      }
    }

    if (!no.ex && length(eval.omit) > 0L) {
      myout$yq <- .npqreg_napredict_eval(eval.omit, myout$yq)
      myout$yqerr <- .npqreg_napredict_eval(eval.omit, myout$yqerr)
      if (gradients) {
        myout$yqgrad <- .npqreg_napredict_eval(eval.omit, myout$yqgrad)
        myout$yqgerr <- .npqreg_napredict_eval(eval.omit, myout$yqgerr)
      }
      txeval <- .npqreg_napredict_eval(eval.omit, txeval)
    }


    fit.elapsed <- proc.time()[3] - fit.start
    optim.time <- if (!is.null(bws$total.time) && is.finite(bws$total.time)) as.double(bws$total.time) else NA_real_
    total.time <- fit.elapsed + (if (is.na(optim.time)) 0.0 else optim.time)

    qregression(bws = bws,
                xeval = txeval,
                tau = tau,
                quantile = myout$yq,
                quanterr = myout$yqerr,
                quantgrad = myout$yqgrad,
                quantgerr = myout$yqgerr,
                ntrain = tnrow,
                trainiseval = no.ex,
                gradients = gradients,
                timing = bws$timing, total.time = total.time,
                optim.time = optim.time, fit.time = fit.elapsed)
  }


npqreg.default <- function(bws, txdat, tydat, nomad = FALSE, ...){
  nomad <- npValidateNomadControl(nomad, "nomad")
  early.dots <- list(...)
  .npqreg_reject_gradient_order_dots(early.dots)
  sc <- sys.call()
  sc.names <- names(sc)

  ## here we check to see if the function was called with tdat =
  ## if it was, we need to catch that and map it to dat =
  ## otherwise the call is passed unadulterated to npudensbw

  bws.named <- any(sc.names == "bws")
  txdat.named <- any(sc.names == "txdat")
  tydat.named <- any(sc.names == "tydat")

  no.bws <- missing(bws)
  no.txdat <- missing(txdat)
  no.tydat <- missing(tydat)
  has.explicit.bws <- (!no.bws) && isa(bws, "condbandwidth")
  bws.formula <- (!no.bws) && inherits(bws, "formula")

  if (bws.named && no.txdat && no.tydat && bws.formula) {
    sc$`bws` <- NULL
    sc$formula <- bws
    sc.bw <- sc
    sc.bw[[1]] <- quote(npcdistbw)
    sc.bw <- .npqreg_strip_fit_controls_from_bw_call(sc.bw)
    bws.named <- FALSE
  } else if (bws.formula) {
    sc.bw <- sc
    sc.bw$`bws` <- NULL
    sc.bw[[1]] <- quote(npcdistbw)
    sc.bw <- .npqreg_strip_fit_controls_from_bw_call(sc.bw)
    bws.named <- FALSE
  } else {
    sc.bw <- sc
    sc.bw[[1]] <- quote(npcdistbw)
    sc.bw <- .npqreg_strip_fit_controls_from_bw_call(sc.bw)
  }

  ## if bws was passed in explicitly, do not compute bandwidths
    
  if(txdat.named)
    txdat <- toFrame(txdat)

  if(tydat.named)
    tydat <- toFrame(tydat)

  if(bws.named){
    sc.bw$bandwidth.compute <- FALSE
  }

  ostxy <- c('txdat','tydat')
  nstxy <- c('xdat','ydat')
  
  m.txy <- match(ostxy, names(sc.bw), nomatch = 0)

  if(any(m.txy > 0)) {
    names(sc.bw)[m.txy] <- nstxy[m.txy > 0]
  }

  use.outer.bandwidth.progress <- !.np_bw_call_uses_nomad_degree_search(
    sc.bw,
    caller_env = parent.frame()
  )

  tbw <- if (!has.explicit.bws) {
    if (use.outer.bandwidth.progress) {
      .np_progress_select_bandwidth_enhanced(
        "Selecting conditional distribution bandwidth",
        .np_eval_bw_call(sc.bw, caller_env = parent.frame())
      )
    } else {
      .np_eval_bw_call(sc.bw, caller_env = parent.frame())
    }
  } else {
    .np_eval_bw_call(sc.bw, caller_env = parent.frame())
  }

  call.args <- list(bws = tbw)
  if (no.bws) {
    call.args$txdat <- txdat
    call.args$tydat <- tydat
  } else {
    if (txdat.named) call.args$txdat <- txdat
    if (tydat.named) call.args$tydat <- tydat
    if ((!bws.named) && (!txdat.named) && (!no.tydat) && (!tydat.named)) {
      call.args <- c(call.args, list(tydat))
    }
  }
  dots <- list(...)
  if (has.explicit.bws)
    fit.dots <- .npqreg_fit_dots(dots)
  else
    fit.dots <- .npqreg_fit_dots(dots, allow.bandwidth.controls = TRUE)
  do.call(npqreg, c(call.args, fit.dots))
}