R/np.condensity.R

In np: Nonparametric Kernel Smoothing Methods for Mixed Data Types

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

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

npcdens.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(npcdens, 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$condens <- napredict(ev$omit, ev$condens)
    ev$conderr <- napredict(ev$omit, ev$conderr)
    if (!is.null(ev$condens.raw))
      ev$condens.raw <- napredict(ev$omit, ev$condens.raw)

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

    return(ev)
  }

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


npcdens.conbandwidth <- 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("project"),
                                 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_condens_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 = "npcdens")
  
  ## 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 = "npcdens")
  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 = "npcdens"
  )
  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("npcdens 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 = "npcdens",
        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 = "npcdens"
    )
  }

  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 = "npcdens"
  )
  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_DENS,
      int_do_tree = .npcdensbw_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 density",
    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")
  )

  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(npcdenshat, 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 <- condensity(bws = bws,
                    xeval = txeval,
                    yeval = tyeval,
                    condens = myout$condens, conderr = myout$conderr,
                    congrad = myout$congrad, congerr = myout$congerr,
                    ll = myout$log_likelihood,
                    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_condens_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 = "npcdens()"
  )

}


npcdens.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, "conbandwidth")
  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(npcdensbw)
    bws.named <- FALSE
  } else {
    sc.bw <- sc
    sc.bw[[1]] <- quote(npcdensbw)
  }

  ## 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 density 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(npcdens, c(call.args, list(...)))
}