[Rcpp-devel] rcpp overhead

[Douglas Bates]

On Sat, Mar 10, 2012 at 6:05 PM, Kaveh Vakili <kaveh.vakili at ulb.ac.be> wrote:
> Hi Steve,
>
> Timing:
>
> i use:
>
> int start_s=clock();
> ..
> ..
> int stop_s=clock();
> cout << "time: " << (stop_s-start_s)/double(CLOCKS_PER_SEC)*1000 << endl;
>
>
> Function:
> Whatever it is, it's coming from these three functions (when i uncoment them, the rest of the code runs in comparable speed).
>
> using namespace Rcpp;
> using namespace Eigen;
> using namespace std;
> using Eigen::MatrixXf;
> using Eigen::VectorXf;
> using Eigen::RowVectorXf;
>
>
> VectorXi SampleR(int& m,int& p){
>        int i,j,n=m;
>        VectorXi x(n);
>        VectorXi y(p);
>        x.setLinSpaced(n,0,n-1);
>        VectorXf urd = VectorXf::Random(p).array().abs();
>        for(i=0;i<p;i++){
>                j=n*urd(i);
>                y(i)=x(j);
>                --n;
>                x(j)=x(n);
>        }
>        return y;
> }
> VectorXf FindLine(MatrixXf& xSub,RowVectorXf& xSub_mean){
>        int h = xSub.rows();
>        int p = xSub.cols();
>        VectorXi RIndex = SampleR(h,p);
>        VectorXf beta(p);
>        Eigen::Matrix<float,16,16>A;
>        for(int i=1;i<p;i++)    A.block(i,0,1,p)=xSub.row(RIndex(i));
>        A.block(0,0,1,p) = xSub_mean;
>        beta = VectorXf::Ones(p);
>        beta = A.topLeftCorner(p,p).lu().solve(beta);
>        beta/=beta.norm();
>        return beta;
> }
> VectorXf OneDir(MatrixXf& x,MatrixXf& xSub,RowVectorXf& xSub_mean,int& h,
> VectorXf& proj){
> VectorXf beta(x.cols());
> beta = FindLine(xSub,xSub_mean);
> proj = ((x*beta).array()-xSub_mean.dot(beta)).square();
> std::nth_element(proj.data(),proj.data()+h,proj.data()+proj.size());
> return log(proj.head(h).mean());
> }
>
>
> What do you think --is there something geeky here?

As others have said, it is not exactly clear what you are doing here
but when you start comparing R and Eigen-based C++ code the natural
approach is to use the RcppEigen package, which also has a plugin for
inline.

I happened to post an example of the use of RcppEigen on
http://dmbates.blogspot.com today.  The example involves sampling from
a collection of multinomial distributions which seems to be related to
what you are doing.

>>On Sat, Mar 10, 2012 at 5:58 PM, Kaveh Vakili <kaveh.vakili at ulb.ac.be> wrote:
>>> Hi all,
>>>
>>> the same code when timed in cpp (i.e. without any interfacing with R)
>>> runs about 2 times faster than when called from R and timed from R's system.time(). In this type of overhead normal or is it a sign i'm doing something not optimally?
>>
>>You'll have to provide the code you are using for test (the plain c++
>>as well as the R/Rcpp hybrid) if you really want smoke this out.
>>
>>My guess is that you're doing something a bit wonky (and it might just
>>be in how you are timing it), but it's impossible to say.
>>
>>-steve
>>
>>--
>>Steve Lianoglou
>>Graduate Student: Computational Systems Biology
>> | Memorial Sloan-Kettering Cancer Center
>> | Weill Medical College of Cornell University
>>Contact Info: http://cbio.mskcc.org/~lianos/contact
>>
>>
>
> _______________________________________________
> 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