如何从另一个 Rccp 函数调用在 Rcpp 中创建的函数

Question

我想在函数 h_evap 中运行一个 Rcpp 函数 esat。两者都在一个通用的 .cpp 文件中，我在 Rstudio 中使用 sourceCPP 执行它。这是代码。 esat 和 h_evap 都被创建并且 esat 工作正常。但是 h_evap 给了我输出

> esat(42)
[1] 256.7082
> h_evap(42)
Error in h_evap(42) : 
  Not compatible with requested type: [type=closure; target=double].

我怀疑问题在于我如何尝试从全局环境访问 esat，但不知道如何调用 esat 来获取输出值而不是闭包。

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]

    NumericVector esat(NumericVector Tk) {
      NumericVector esat_out(Tk.size(), NAN);
      for (size_t i=0; i<Tk.size(); i++) {
        esat_out[i] = 6.1121 * Tk[i];
      }
      return esat_out;
    }
    
    // [[Rcpp::export]]
    NumericVector h_evap(NumericVector Tk) {
      Environment env = Environment::global_env();
      NumericVector f_esat = env["esat"];
      NumericVector h_evap_out(Tk.size(), NAN);
      for (size_t i=0; i<Tk.size(); i++) {
        
        h_evap_out[i] = (313.15 - Tk[i]);
        h_evap_out[i] = h_evap_out[i] + f_esat(Tk[i]);
      }

  return h_evap_out;
}

/*** R
h_evap(42)
*/

另一种方法是使用 cppFunction。我已经尝试过了，但仍然得到这个 Rcpp 新手不清楚的错误。这是代码

library(Rcpp)

cppFunction('NumericVector esat(NumericVector Tk) {
  NumericVector esat_out(Tk.size(), NAN);
  for (size_t i=0; i<Tk.size(); i++) {
    esat_out[i] = 6.1121 * Tk[i];
  }
  return esat_out;
}')

cppFunction('NumericVector h_evap(NumericVector Tk) {
   NumericVector h_evap_out(Tk.size(), NAN);
  for (size_t i=0; i<Tk.size(); i++) {
    h_evap_out[i] = esat(Tk[i]);
  }
  return h_evap_out;
}')

esat 编译得很好。 h_evap 返回一个错误信息我不清楚...

Answer 1

稍微重写您的文件以避免通过中间 R 函数调用 C++，这（通常）是一个坏主意，而且几乎总是对性能造成不必要的沉重负担。

正如您在同一文件中定义一个有效的 C++ 函数在使用它之前（这样您就不需要签名来将其声明为例如头文件会为你做的）可以简单地调用它。

我还更改了循环索引变量以在编译期间消除一个警告，并且，当我这样做时，删除了 using namespace Rcpp; 并切换到使用命名空间的显式调用，这更明确，更“安全”一点在更大的代码库中。

编辑： 由于您的循环实际上对循环索引是不变的， 我们可以将代码重写为 vectorized 调用，这样更短、更简单、更快且更易于推理。 （当然，也可以从 R 完成...）

代码

#include <Rcpp.h>

// [[Rcpp::export]]
Rcpp::NumericVector esat(Rcpp::NumericVector Tk) {
    Rcpp::NumericVector esat_out(Tk.size(), NAN);
    for (R_xlen_t i=0; i<Tk.size(); i++) {
        esat_out[i] = 6.1121 * Tk[i];
    }
    return esat_out;
}

// [[Rcpp::export]]
Rcpp::NumericVector h_evap(Rcpp::NumericVector Tk) {
    Rcpp::NumericVector h_evap_out(Tk.size(), NAN);
    Rcpp::NumericVector f_out = esat(Tk);
    for (R_xlen_t i=0; i<Tk.size(); i++) {
        h_evap_out[i] = (313.15 - Tk[i]);
        h_evap_out[i] = h_evap_out[i] + f_out[i];
    }
    return h_evap_out;
}

// [[Rcpp::export]]
Rcpp::NumericVector esatV(Rcpp::NumericVector Tk) {
    Rcpp::NumericVector esat_out = 6.1121 * Tk;
    return esat_out;
}

// [[Rcpp::export]]
Rcpp::NumericVector h_evapV(Rcpp::NumericVector Tk) {
  Rcpp::NumericVector f_out = esatV(Tk);
  Rcpp::NumericVector h_evap_out = 313.15 - Tk + f_out;
  return h_evap_out;
}

/*** R
esat(42)
h_evap(42)
esatV(42)
h_evapV(42)
*/

用法

> Rcpp::sourceCpp("~/git/stackoverflow/68605528/answer.cpp")

> esat(42)
[1] 256.708

> h_evap(42)
[1] 527.858

> esatV(42)
[1] 256.708

> h_evapV(42)
[1] 527.858
> 

Answer 2

通过对各种堆栈溢出 Q 和 As 的一些明智研究以及 Dirk 的一些关键提示，以下代码似乎可以满足我的需求。两个关键点 - 使用 Function 并将 esat 的输出写入 NumericVector 并将其读入 h_evap。我还在做包...

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]

NumericVector esat(NumericVector Tk) {
  NumericVector esat_out(Tk.size(), NAN);
  for (size_t i=0; i<Tk.size(); i++) {
    esat_out[i] = 6.1121 * Tk[i];
  }
  return esat_out;
}

// [[Rcpp::export]]
NumericVector h_evap(NumericVector Tk) {
  Environment env = Environment::global_env();
  Function f_esat = env["esat"];
  NumericVector h_evap_out(Tk.size(), NAN);
  NumericVector f_out = f_esat(Tk);
  for (size_t i=0; i<Tk.size(); i++) {

    h_evap_out[i] = (313.15 - Tk[i]);
    h_evap_out[i] = h_evap_out[i] + f_out[i];
  }

  return h_evap_out;
}

/*** R
esat(42)
h_evap(42)
*/