how to call a function created in Rcpp from another Rccp function

Question

I want to run an Rcpp function esat in the function h_evap. Both are in a common .cpp file and I execute it with sourceCPP in Rstudio. Here's the code. Both esat and h_evap are created and esat works fine. But h_evap gives me the output

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

I suspect the problem is in how I try to access esat from the global environment but can't figure out how to call esat to get the value of the output rather than a closure.

#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)
*/

An alternative is to use cppFunction. I've tried that and still get errors that are not clear to this Rcpp novice. Here's the code

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 compiles fine. h_evap returns an error message is not clear to me...

Answer 1

A slight rewrite of your file to avoid calling a C++ via an intermediate R function which is (generally) a bad idea and almost always an uncalled-for and heavy tax on performance.

As you defined a valid C++ function in the same file and before its use (so that you don't need a signature to declare it as e.g. a header file would do for you) can simply call it.

I also changed the loop index variable to get rid of one warning during compilation, and, while I was at it, removed using namespace Rcpp; and switched to explicit calls with namespace which is more explicit and a little 'safer' from surprises in larger code bases.

Edit: And as your loops are in fact invariant to the loop index, we can rewrite the code as vectorized calls whicg is shorter, simpler, faster, and easier to reason with. (And could, of course, be done from R too...)

Code

#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)
*/

Usage

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