Search code examples
rrcpp

R use stats::optimize in Rcpp, simple example fails to compile


I want to use the R stats::optimize function in Rcpp because I haven't been able to find an Rcpp equivalent. The code below is my attempt at a simple example based on the Example in the optimize help, but fails. Here's the R function and results

f <- function (x) (x - .33)^2
xmin <- optimize(f, c(0, 1), tol = 0.0001)
xmin

This returns

$minimum
[1] 0.333

$objective
[1] 0

Here's the Rcpp code that fails when sourcing it.

    #include <Rcpp.h>
    const double tolerance = 1e-0;

   // [[Rcpp::export]]
    Rcpp::NumericVector f(Rcpp::NumericVector x) {
      return pow(x-0.33, 2);
    }
    
    Rcpp::List fTg_opt(const double optmin, const double optmax) {
      Rcpp::Environment base("package:stats");
      Rcpp::Function optimize_r = base["optimize"];    
      Rcpp::NumericVector interval = {optmin,optmax};
      return optimize_r(f, interval, tolerance);
    }

The Rstudio console has the following error messages.

> Rcpp::sourceCpp("R/cpp/testopt.cpp")
In file included from testopt.cpp:1:
In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp.h:27:
In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/RcppCommon.h:157:
In file included from /Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/traits/traits.h:45:
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/traits/is_convertible.h:35:10: error: function cannot return function type 'Rcpp::Vector<14> (Rcpp::Vector<14>)'
                static T MakeT() ;
                       ^
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/internal/wrap.h:770:75: note: in instantiation of template class 'Rcpp::traits::is_convertible<Rcpp::Vector<14> (Rcpp::Vector<14>), SEXPREC *>' requested here
            return wrap_dispatch_unknown(object, typename ::Rcpp::traits::is_convertible<T,SEXP>::type());
                                                                          ^
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/internal/wrap.h:787:20: note: in instantiation of function template specialization 'Rcpp::internal::wrap_dispatch_eigen<Rcpp::Vector<14> (Rcpp::Vector<14>)>' requested here
            return wrap_dispatch_eigen(object, typename traits::is_eigen_base<T>::type());
                   ^
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/internal/wrap.h:807:20: note: in instantiation of function template specialization 'Rcpp::internal::wrap_dispatch_unknown_importable<Rcpp::Vector<14> (Rcpp::Vector<14>)>' requested here
            return wrap_dispatch_unknown_importable(object, typename ::Rcpp::traits::is_importer<T>::type());
                   ^
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/internal/wrap_end.h:30:25: note: in instantiation of function template specialization 'Rcpp::internal::wrap_dispatch<Rcpp::Vector<14> (Rcpp::Vector<14>)>' requested here
          return internal::wrap_dispatch( object, typename ::Rcpp::traits::wrap_type_traits<T>::wrap_category() ) ;
                           ^
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/grow.h:44:26: note: in instantiation of function template specialization 'Rcpp::wrap<Rcpp::Vector<14> (Rcpp::Vector<14>)>' requested here
            return grow( wrap(head), tail ) ;
                         ^
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/grow.h:65:26: note: in instantiation of function template specialization 'Rcpp::internal::grow__dispatch<Rcpp::Vector<14> (Rcpp::Vector<14>)>' requested here
        return internal::grow__dispatch(typename traits::is_named<T>::type(), head, y);
                         ^
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/generated/grow__pairlist.h:45:9: note: in instantiation of function template specialization 'Rcpp::grow<Rcpp::Vector<14> (Rcpp::Vector<14>)>' requested here
        return grow( t1, grow( t2, grow( t3, R_NilValue ) ) ) ;
               ^
/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include/Rcpp/generated/Function__operator.h:45:20: note: in instantiation of function template specialization 'Rcpp::pairlist<Rcpp::Vector<14> (Rcpp::Vector<14>), Rcpp::Vector<14>, double>' requested here
            return invoke(pairlist(t1, t2, t3), R_GlobalEnv);
                          ^
testopt.cpp:13:20: note: in instantiation of function template specialization 'Rcpp::Function_Impl<PreserveStorage>::operator()<Rcpp::Vector<14> (Rcpp::Vector<14>), Rcpp::Vector<14>, double>' requested here
  return optimize_r(f, interval, tolerance);
                   ^
1 error generated.
make: *** [testopt.o] Error 1
clang++ -mmacosx-version-min=10.13 -std=gnu++14 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.1/Resources/library/Rcpp/include" -I"/Users/gcn/Documents/workspace/ISIMIPData/R/cpp" -I/usr/local/include   -fPIC  -Wall -g -O2  -c testopt.cpp -o testopt.o
Error in Rcpp::sourceCpp("R/cpp/testopt.cpp") : 
  Error 1 occurred building shared library.

Solution

  • One of your problems here is that you assume that becomes a function you submit to compilation under Rcpp::sourceCpp() is callable under its exported name.

    It is not. Try Rcpp::sourceCpp(..., verbose=TRUE), i.e. add that arguments, to see what is really called. Those you could pass around (using SEXP argunments and results, but they are unwieldy).

    To prove, here is a 'working but useless' version of your code. If we pass f() from R too, everything is callable.

    Morale: The interface still is SEXP .Call("name", SEXP a, SEXP b, ...) even if Rcpp hides that. No Free Lunch (TM). But as my comment hinted, there are optimization packages you can use with Rcpp.

    Code

    #include <Rcpp.h>
    
    // [[Rcpp::export]]
    Rcpp::List fTg_opt(Rcpp::Function f, const double optmin, const double optmax) {
        Rcpp::Environment base("package:stats");
        Rcpp::Function optimize_r = base["optimize"];
        Rcpp::NumericVector interval = {optmin,optmax};
        Rcpp::List res = optimize_r(f, interval);
        return res;
    }
    
    /*** R
    f <- function (x) (x - .33)^2
    xmin <- optimize(f, c(0, 1), tol = 0.0001)
    xmin
    
    fTg_opt(f, 0, 1)
    */
    

    Output

    > Rcpp::sourceCpp("~/git/stackoverflow/68674076/question.cpp")
    
    > f <- function (x) (x - .33)^2
    
    > xmin <- optimize(f, c(0, 1), tol = 0.0001)
    
    > xmin
    $minimum
    [1] 0.33
    
    $objective
    [1] 0
    
    
    > fTg_opt(f, 0, 1)
    $minimum
    [1] 0.33
    
    $objective
    [1] 0