如何从 RCPP 调用 model.matrix 或等价物,可能在线程代码中?

How to call model.matrix or equivalent from RCPP, possibly in threaded code?

我们希望在具有许多结果不相互依赖的循环的算法中使用线程来加快速度。

在我们希望移植到 rcpp 的代码中,有一个对 model.matrix.

的调用

这似乎并不容易移植。

进一步调查(关于为我们的用例运行什么代码),揭示了 lm 对象的 S3 方法对变量做了一些准备工作,然后调用函数的默认版本,如在此代码的复制粘贴:

function (object, ...) 
{
    if (n_match <- match("x", names(object), 0L)) 
        object[[n_match]]
    else {
        data <- model.frame(object, xlev = object$xlevels, ...)
        if (exists(".GenericCallEnv", inherits = FALSE)) 
            NextMethod("model.matrix", data = data, contrasts.arg = object$contrasts)
        else {
            dots <- list(...)
            dots$data <- dots$contrasts.arg <- NULL
            do.call("model.matrix.default", c(list(object = object, 
                data = data, contrasts.arg = object$contrasts), 
                dots))
        }
    }
}

该函数的默认版本至少将其部分功能提供给已编译的 C 函数:

function (object, data = environment(object), contrasts.arg = NULL, 
    xlev = NULL, ...) {
    t <- if (missing(data)) 
        terms(object)
    else terms(object, data = data)
    if (is.null(attr(data, "terms"))) 
        data <- model.frame(object, data, xlev = xlev)
    else {
        reorder <- match(vapply(attr(t, "variables"), deparse2, 
            "")[-1L], names(data))
        if (anyNA(reorder)) 
            stop("model frame and formula mismatch in model.matrix()")
        if (!identical(reorder, seq_len(ncol(data)))) 
            data <- data[, reorder, drop = FALSE]
    }
    int <- attr(t, "response")
    if (length(data)) {
        contr.funs <- as.character(getOption("contrasts"))
        namD <- names(data)
        for (i in namD) if (is.character(data[[i]])) 
            data[[i]] <- factor(data[[i]])
        isF <- vapply(data, function(x) is.factor(x) || is.logical(x), 
            NA)
        isF[int] <- FALSE
        isOF <- vapply(data, is.ordered, NA)
        for (nn in namD[isF]) if (is.null(attr(data[[nn]], "contrasts"))) 
            contrasts(data[[nn]]) <- contr.funs[1 + isOF[nn]]
        if (!is.null(contrasts.arg)) {
            if (!is.list(contrasts.arg)) 
                warning("non-list contrasts argument ignored")
            else {
                if (is.null(namC <- names(contrasts.arg))) 
                  stop("'contrasts.arg' argument must be named")
                for (nn in namC) {
                  if (is.na(ni <- match(nn, namD))) 
                    warning(gettextf("variable '%s' is absent, its contrast will be ignored", 
                      nn), domain = NA)
                  else {
                    ca <- contrasts.arg[[nn]]
                    if (is.matrix(ca)) 
                      contrasts(data[[ni]], ncol(ca)) <- ca
                    else contrasts(data[[ni]]) <- contrasts.arg[[nn]]
                  }
                }
            }
        }
    }
    else {
        isF <- FALSE
        data[["x"]] <- raw(nrow(data))
    }
    ans <- .External2(C_modelmatrix, t, data)
    if (any(isF)) 
        attr(ans, "contrasts") <- lapply(data[isF], attr, 
            "contrasts")
    ans
}

是否有某种方式可以从 Rcpp 调用 C_modelmatrix,无论是单线程还是多线程?是否有任何库或包在 Rcpp 中做基本相同的事情,所以我不必在这里重新发明轮子?如果可以避免,我宁愿不必完全重新实现 model.matrix 所做的一切。

因为我们实际上没有功能代码,所以还没有任何可以展示的代码。

我们试图加快调用函数的相关部分 model.matrix,如下所示:("model.y is an lm",数据都是 model.frame(model.y) )

ymat.t <- model.matrix(terms(model.y), data=pred.data.t)
ymat.c <- model.matrix(terms(model.y), data=pred.data.c)

这实际上不是一个基于结果的问题,更像是一个基于 approach/methods 的问题

您可以从 C++ 中调用 model.matrix,但不能以多线程方式调用。

也会有开销,但如果在代码中间深处需要函数调用,那么为了方便起见,这样做是值得的。

示例:

// [[Rcpp::export]]
RObject call(RObject x, RObject y){
  Environment env = Environment::global_env();
  Function f = env["model.matrix"];
  RObject res = f(x,y);
  return res;
}