[Rcpp-devel] use of auxiliary functions

[Douglas Bates]

On Fri, Aug 6, 2010 at 3:14 PM, baptiste auguie
<baptiste.auguie at googlemail.com> wrote:
> Dear all,

> In a package using RcppArmadillo I define several functions that use
> the same piece of c++ code. For instance, several functions use a
> rotation matrix constructed from 3 angles. I am considering the
> following idea: 1) define a routine that performs the calculations; 2)
> create a wrapper in case I want to access it at R level; 3) simply
> call this routine in other c++ functions. The obvious advantage is
> that I minimize code duplication and errors associated with it. It is
> perhaps best illustrated with the code,

Code reuse is a good idea.

>  arma::mat euler(double phi, double theta, double psi)
>  {
>    arma::mat Rot(3,3);
>    Rot(0,0) = cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi);
>    Rot(0,1) = cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi);
>    Rot(0,2) = sin(psi)*sin(theta);
>
>    Rot(1,0) = -sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi);
>    Rot(1,1) = -sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi);
>    Rot(1,2) = cos(psi)*sin(theta);
>
>    Rot(2,0) = sin(phi)*sin(theta);
>    Rot(2,1) = -cos(phi)*sin(theta);
>    Rot(2,2) = cos(theta);
>    return (Rot);
>  }

It is probably not a big deal now but you are recalculating the same
trigonometric functions several times in there.  Why not calculate the
cos and sin of the each of the angles just once and then form the
elements of the rotation matrix as linear combinations?  (Once upon a
time we used to keep track of every floating point calculation and try
to work out ways to minimize that number.)

> // this is the wrapper for access at R level
>  SEXP rotation(SEXP Iphi, SEXP Itheta, SEXP Ipsi)
>  {
>    try {
>
>      // conversion of input variables
>      double phi  =  as<double>(Iphi), theta  =  as<double>(Itheta),
> psi  =  as<double>(Ipsi);
>      arma::mat Rot= euler(theta, phi, psi);
>      List res ;
>      res["R"] = Rot;
>
>      return res ;
>
>    } catch( std::exception &ex ) {
>      forward_exception_to_r( ex );
>    } catch(...) {
>      ::Rf_error( "c++ exception (unknown reason)" );
>    }
>    return R_NilValue; // -Wall
>  }

Myself I would use Romain's RCPP_FUNCTION_3 macro and write that as

RCPP_FUNCTION_3(NumericMatrix, rotation, NumericVector phi,
NumericVector theta, NumericVector psi) {
  return wrap(euler(*phi.begin(), *theta.begin(), *psi.begin()));
}

It would be even easier if you defined euler to return a NumericMatrix
and converted to the arma matrix only when needed.

You are passing the result back as a list of one element which may
make sense in context but generally is an unnecessary wrapper.

> // and perhaps
> // other SEXP functions using euler()
> // ...
>
> I am wondering if there is any downside to this approach, in terms of
> efficiency. In particular, when rotation() calls euler(), is it making
> copies of the variables in memory? Should I use pointers instead? (I
> have forgotten the most basic things about them).
> I have made a quick, inconclusive timing comparison of both
> techniques; on this particular example the cost would be negligible (I
> have other routines that use large matrices, however).
>
>
> Regards,
>
> baptiste
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
>