R/np.condistribution.R

In np: Nonparametric Kernel Smoothing Methods for Mixed Data Types

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

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

npcdist.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 (!missing(data) && !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) {
      npValidateNewdataFormula(newdata, tt, include.response = TRUE)
      umf.args <- list(formula = tt, data = newdata)
      umf <- do.call(stats::model.frame, umf.args, envir = parent.frame())
      emf <- umf
      
      eydat <- emf[, bws$variableNames[["response"]], drop = FALSE]
      exdat <- emf[, bws$variableNames[["terms"]], drop = FALSE]
    }

    cd.args <- list(txdat = txdat, tydat = tydat)
    if (has.eval) {
      cd.args$exdat <- exdat
      cd.args$eydat <- eydat
    }
    cd.args$bws <- bws
    ev <- do.call(npcdist, c(cd.args, list(...)))

    ev$omit <- attr(umf,"na.action")
    ev$rows.omit <- as.vector(ev$omit)
    ev$nobs.omit <- length(ev$rows.omit)
    train.omit <- as.vector(attr(tmf, "na.action"))
    eval.omit <- if (has.eval) as.vector(attr(umf, "na.action")) else integer(0)
    ev$train.rows.omit <- if (length(train.omit)) train.omit else NA
    ev$train.nobs.omit <- length(train.omit)
    ev$eval.rows.omit <- if (length(eval.omit)) eval.omit else NA
    ev$eval.nobs.omit <- length(eval.omit)

    ev$condist <- napredict(ev$omit, ev$condist)
    ev$conderr <- napredict(ev$omit, ev$conderr)
    if (!is.null(ev$condist.raw))
      ev$condist.raw <- napredict(ev$omit, ev$condist.raw)

    if(ev$gradients){
        ev$congrad <- napredict(ev$omit, ev$congrad)
        ev$congerr <- napredict(ev$omit, ev$congerr)
    }

    return(ev)
  }

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


npcdist.condbandwidth <-
  function(bws,
           txdat = stop("invoked without training data 'txdat'"),
           tydat = stop("invoked without training data 'tydat'"),
           exdat, eydat, gradients = FALSE,
           gradient.order = 1L,
           proper = FALSE,
           proper.method = c("isotonic"),
           proper.control = list(),
           ...){

    dots <- list(...)
    fit.start <- proc.time()[3]
    fit.progress.handoff <- isTRUE(dots$.np_fit_progress_handoff)
    gradients <- npValidateScalarLogical(gradients, "gradients")
    proper.args <- .np_condist_validate_proper_args(
      proper = proper,
      proper.method = proper.method,
      proper.control = proper.control
    )

    if (xor(missing(exdat),missing(eydat)))
      stop("evaluation data must be supplied for both 'exdat' and 'eydat'")

    no.exy = missing(exdat)

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

    if (!no.exy){
      exdat = toFrame(exdat)
      eydat = toFrame(eydat)

      if (! txdat %~% exdat )
        stop("'txdat' and 'exdat' are not similar data frames!")

      if (! tydat %~% eydat )
        stop("'tydat' and 'eydat' 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) != length(tydat))
      stop("length of bandwidth vector does not match number of columns of 'tydat'")

    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")

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

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

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

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

    eval.rows.omit <- integer(0)
    if (!no.exy){
      keep.eval <- rep_len(TRUE, nrow(exdat))
      eval.rows.omit <- attr(na.omit(data.frame(exdat, eydat)), "na.action")
      if (length(eval.rows.omit) > 0L)
        keep.eval[as.integer(eval.rows.omit)] <- FALSE

      if (!any(keep.eval))
        stop("Data has no rows without NAs")

      exdat <- exdat[keep.eval,,drop = FALSE]
      eydat <- eydat[keep.eval,,drop = FALSE]
    }


    tnrow = nrow(txdat)
    enrow = (if (no.exy) tnrow else nrow(exdat))

    ## 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.exy){
      exdat <- adjustLevels(exdat, bws$xdati, allowNewCells = TRUE)
      eydat <- adjustLevels(eydat, bws$ydati, allowNewCells = TRUE)
      npKernelBoundsCheckEval(exdat, bws$ixcon, bws$cxkerlb, bws$cxkerub, argprefix = "cxker")
      npKernelBoundsCheckEval(eydat, bws$iycon, bws$cykerlb, bws$cykerub, argprefix = "cyker")
    }

    proper.slice.context <- list(
      txdat = txdat,
      tydat = tydat,
      exdat = if (no.exy) NULL else exdat,
      eydat = if (no.exy) NULL else eydat
    )

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

    ## at this stage, data to be sent to the c routines must be converted to
    ## numeric type.
    
    tydat = toMatrix(tydat)

    tyuno = tydat[, bws$iyuno, drop = FALSE]
    tycon = tydat[, bws$iycon, drop = FALSE]
    tyord = tydat[, bws$iyord, drop = FALSE]


    txdat = toMatrix(txdat)

    txuno = txdat[, bws$ixuno, drop = FALSE]
    txcon = txdat[, bws$ixcon, drop = FALSE]
    txord = txdat[, bws$ixord, drop = FALSE]

    if (!no.exy){
      eydat = toMatrix(eydat)

      eyuno = eydat[, bws$iyuno, drop = FALSE]
      eycon = eydat[, bws$iycon, drop = FALSE]
      eyord = eydat[, bws$iyord, drop = FALSE]


      exdat = toMatrix(exdat)

      exuno = exdat[, bws$ixuno, drop = FALSE]
      excon = exdat[, bws$ixcon, drop = FALSE]
      exord = exdat[, bws$ixord, drop = FALSE]
    } else {
      eyuno = data.frame()
      eycon = data.frame()
      eyord = data.frame()

      exuno = data.frame()
      excon = data.frame()
      exord = data.frame()
    }

    reg.spec <- npConditionalRegEngineSpec(bws, where = "npcdist")
    reg.engine <- reg.spec$reg.engine
    basis.engine <- reg.spec$basis.engine
    degree.engine <- reg.spec$degree.engine
    bernstein.engine <- reg.spec$bernstein.engine
    glp.gradient.order <- npConditionalGradientOrder(
      bws = bws,
      reg.engine = reg.engine,
      gradient.order = gradient.order,
      where = "npcdist"
    )
    lp.degree0.lc.gradient <- isTRUE(gradients) &&
      npGlpDegree0FirstDerivativeLcOk(
        regtype.engine = reg.engine,
        degree.engine = degree.engine,
        gradient.order = glp.gradient.order,
        ncon = bws$xncon
      )
    glp.gradient.available <- NULL
    glp.gradient.partial <- FALSE
    if (isTRUE(gradients) &&
        identical(reg.engine, "lp") &&
        (bws$xncon > 0L)) {
      glp.gradient.available <- npGlpGradientAvailability(
        regtype.engine = reg.engine,
        degree.engine = degree.engine,
        gradient.order = glp.gradient.order,
        ncon = bws$xncon
      )
      if (!any(glp.gradient.available)) {
        stop("npcdist has no available derivative components for the requested gradient.order and fitted polynomial degree",
             call. = FALSE)
      }
      glp.gradient.partial <- !lp.degree0.lc.gradient &&
        any(!glp.gradient.available)
      if (glp.gradient.partial && any(!glp.gradient.available)) {
        npWarnGlpGradientPartialAvailability(
          where = "npcdist",
          degree.engine = degree.engine,
          gradient.order = glp.gradient.order,
          available = glp.gradient.available,
          con.names = colnames(txdat)[bws$ixcon]
        )
      }
    }
    if (isTRUE(gradients) &&
        identical(reg.engine, "lp") &&
        (bws$xncon > 0L) &&
        !lp.degree0.lc.gradient &&
        all(degree.engine == 0L)) {
      stop("regtype='lp' with degree=0 does not support derivatives; use gradients=FALSE for fitted/predicted values")
    }
    if (isTRUE(gradients) &&
        identical(reg.engine, "lp") &&
        !glp.gradient.partial) {
      npValidateGlpGradientDegree(
        regtype.engine = reg.engine,
        degree.engine = degree.engine,
        gradient.order = glp.gradient.order,
        ncon = bws$xncon,
        where = "npcdist"
      )
    }

    reg.c <- npRegtypeToC(
      regtype = if (identical(reg.engine, "lp") && !lp.degree0.lc.gradient) "lp" else "lc",
      degree = if (lp.degree0.lc.gradient) rep.int(0L, bws$xncon) else degree.engine,
      ncon = bws$xncon,
      context = "npcdist"
    )
    degree.c <- if (bws$xncon > 0L) {
      as.integer(if (is.null(reg.c$degree)) rep.int(0L, bws$xncon) else reg.c$degree)
    } else {
      integer(0)
    }
    basis.code <- as.integer(npLpBasisCode(basis.engine))
    do.compiled.gradients <- isTRUE(gradients) && !glp.gradient.partial

    myopti <- list(
        num_obs_train = tnrow,
        num_obs_eval = enrow,
        int_LARGE_SF = (if (bws$scaling) SF_NORMAL else SF_ARB),
        BANDWIDTH_den_extern = switch(bws$type,
            fixed = BW_FIXED,
            generalized_nn = BW_GEN_NN,
            adaptive_nn = BW_ADAP_NN),
        int_MINIMIZE_IO=if (isTRUE(getOption("np.messages"))) IO_MIN_FALSE else IO_MIN_TRUE,
        xkerneval = switch(bws$cxkertype,
            gaussian = CKER_GAUSS + bws$cxkerorder/2 - 1,
            epanechnikov = CKER_EPAN + bws$cxkerorder/2 - 1,
            uniform = CKER_UNI,
            "truncated gaussian" = CKER_TGAUSS),
        ykerneval = switch(bws$cykertype,
            gaussian = CKER_GAUSS + bws$cykerorder/2 - 1,
            epanechnikov = CKER_EPAN + bws$cykerorder/2 - 1,
            uniform = CKER_UNI,
            "truncated gaussian" = CKER_TGAUSS),
        uxkerneval = switch(bws$uxkertype,
            aitchisonaitken = UKER_AIT,
            liracine = UKER_LR),
        uykerneval = switch(bws$uykertype,
            aitchisonaitken = UKER_AIT,
            liracine = UKER_LR),
        oxkerneval = switch(bws$oxkertype,
            wangvanryzin = OKER_WANG,
            liracine = OKER_NLR,
        "racineliyan" = OKER_RLY),
        oykerneval = switch(bws$oykertype,
            wangvanryzin = OKER_WANG,
            liracine = OKER_NLR,
        "racineliyan" = OKER_RLY),
        num_yuno = bws$ynuno,
        num_yord = bws$ynord,
        num_ycon = bws$yncon,
        num_xuno = bws$xnuno,
        num_xord = bws$xnord,
        num_xcon = bws$xncon,
        no.exy = no.exy,
        gradients = do.compiled.gradients,
        ymcv.numRow = attr(bws$ymcv, "num.row"),
        xmcv.numRow = attr(bws$xmcv, "num.row"),
        densOrDist = NP_DO_DIST,
        int_do_tree = .npcdistbw_tree_code(
          bws = bws,
          ncon = bws$yncon + bws$xncon,
          ncat = bws$ynuno + bws$ynord + bws$xnuno + bws$xnord))

    cxker.bounds.c <- npKernelBoundsMarshal(bws$cxkerlb[bws$ixcon], bws$cxkerub[bws$ixcon])
    cyker.bounds.c <- npKernelBoundsMarshal(bws$cykerlb[bws$iycon], bws$cykerub[bws$iycon])

    myout <- .np_with_compiled_fit_progress(
      label = "Fitting conditional distribution",
      total = .np_condensdist_fit_total(bws = bws, tnrow = tnrow, enrow = enrow),
      handoff = fit.progress.handoff,
      handoff.detail = if (fit.progress.handoff) "starting" else NULL,
      .Call("C_np_density_conditional",
            as.double(tyuno), as.double(tyord), as.double(tycon),
            as.double(txuno), as.double(txord), as.double(txcon),
            as.double(eyuno), as.double(eyord), as.double(eycon),
            as.double(exuno), as.double(exord), as.double(excon),
            as.double(c(bws$xbw[bws$ixcon], bws$ybw[bws$iycon],
                        bws$ybw[bws$iyuno], bws$ybw[bws$iyord],
                        bws$xbw[bws$ixuno], bws$xbw[bws$ixord])),
            as.double(bws$ymcv), as.double(attr(bws$ymcv, "pad.num")),
            as.double(bws$xmcv), as.double(attr(bws$xmcv, "pad.num")),
            as.double(bws$nconfac), as.double(bws$ncatfac), as.double(bws$sdev),
            as.integer(myopti),
            as.integer(enrow),
            as.integer(bws$xndim),
            as.double(cxker.bounds.c$lb),
            as.double(cxker.bounds.c$ub),
            as.double(cyker.bounds.c$lb),
            as.double(cyker.bounds.c$ub),
            as.integer(reg.c$code),
            as.integer(degree.c),
            as.integer(bernstein.engine),
            basis.code,
            PACKAGE = "np")
    )
    names(myout)[1] <- "condist"

    if(gradients){
      if (!glp.gradient.partial) {
        myout$congrad = matrix(data=myout$congrad, nrow = enrow, ncol = bws$xndim, byrow = FALSE)
        rorder = numeric(bws$xndim)
        xidx <- seq_len(bws$xndim)
        rorder[c(xidx[bws$ixcon], xidx[bws$ixuno], xidx[bws$ixord])] <- xidx
        myout$congrad = myout$congrad[, rorder, drop = FALSE]

        myout$congerr = matrix(data=myout$congerr, nrow = enrow, ncol = bws$xndim, byrow = FALSE)
        myout$congerr = myout$congerr[, rorder, drop = FALSE]
      } else {
        myout$congrad <- matrix(NA_real_, nrow = enrow, ncol = bws$xndim)
        myout$congerr <- matrix(NA_real_, nrow = enrow, ncol = bws$xndim)
      }

      if (identical(reg.engine, "lp") && bws$xncon > 0L && !lp.degree0.lc.gradient) {
        cont.idx <- which(bws$ixcon)
        hat.gradient.idx <- if (glp.gradient.partial) {
          which(glp.gradient.available)
        } else {
          which(glp.gradient.order > 1L)
        }
        if (length(hat.gradient.idx)) {
          rhs <- rep.int(1.0, nrow(proper.slice.context$txdat))
          for (jj in hat.gradient.idx) {
            svec <- integer(bws$xncon)
            svec[jj] <- glp.gradient.order[jj]
            hat.args <- list(
              bws = bws,
              txdat = proper.slice.context$txdat,
              tydat = proper.slice.context$tydat,
              y = rhs,
              output = "apply",
              s = svec
            )
            if (!no.exy) {
              hat.args$exdat <- proper.slice.context$exdat
              hat.args$eydat <- proper.slice.context$eydat
            }
            myout$congrad[, cont.idx[jj]] <- as.vector(do.call(npcdisthat, hat.args))
            myout$congerr[, cont.idx[jj]] <- NA_real_
          }
        }
      }
    } else {
      myout$congrad = NA
      myout$congerr = NA
    }


    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)

    out <- condistribution(bws = bws,
                           xeval = txeval,
                           yeval = tyeval,
                           condist = myout$condist, conderr = myout$conderr,
                           congrad = myout$congrad, congerr = myout$congerr,
                           ntrain = tnrow, trainiseval = no.exy, gradients = gradients,
                           gradient.order = if (identical(reg.engine, "lp")) glp.gradient.order else NULL,
                           rows.omit = if (no.exy) train.rows.omit else eval.rows.omit,
                           train.rows.omit = train.rows.omit,
                           eval.rows.omit = if (no.exy) integer(0) else eval.rows.omit,
                           timing = bws$timing, total.time = total.time,
                           optim.time = optim.time, fit.time = fit.elapsed)
    out$nomad.time <- if (!is.null(bws$nomad.time) && is.finite(bws$nomad.time)) as.double(bws$nomad.time) else NA_real_
    out$powell.time <- if (!is.null(bws$powell.time) && is.finite(bws$powell.time)) as.double(bws$powell.time) else NA_real_

    .np_condist_finalize_proper_object(
      object = out,
      proper = proper.args$proper.requested,
      proper.method = proper.args$proper.method,
      proper.control = proper.args$proper.control,
      slice.context = proper.slice.context,
      where = "npcdist()"
    )

  }

npcdist.default <- function(bws, txdat, tydat, nomad = FALSE, ...){
  sc <- sys.call()
  sc.names <- names(sc)
  nomad <- npValidateNomadControl(nomad, "nomad")

  ## 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)
    bws.named <- FALSE
  } else {
    sc.bw <- sc
    sc.bw[[1]] <- quote(npcdistbw)
  }

  ## 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]
  }
  sc.bw$newdata <- NULL
  sc.bw$exdat <- NULL
  sc.bw$eydat <- NULL
  sc.bw$gradients <- NULL
  sc.bw$gradient.order <- NULL
  sc.bw$proper <- NULL
  sc.bw$proper.method <- NULL
  sc.bw$proper.control <- NULL
    
  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))
    }
  }
  if (!has.explicit.bws)
    call.args$.np_fit_progress_handoff <- TRUE
  do.call(npcdist, c(call.args, list(...)))
}