Le 04/12/12 21:27, Honglang Wang a écrit :
Since the "ker" function defined in R is very simple and vectorized, I
thought it was fine to call this function in C++. And I tried to code
this function inside C++, and I found it became ever slower: from
47.310s to 49.536s.
Hard to say what is going on. Perhaps you are allocating too many things
with :
arma::vec tp = t-t0;
arma::uvec q1 = arma::find(arma::abs(tp)<h);
arma::uvec q2 = arma::find(arma::abs(tp)>=h);
arma::vec kk = t;
kk.elem(q1) = ((1-arma::pow(tp.elem(q1)/h,2))/h)*0.75;
kk.elem(q2).zeros();
You are allocating tp, q1, q2 and kk.
Perhaps, there are ways not to do that.
arma brings you nice syntax, but sometimes this has a cost.
Sometimes, well written R code performs well and the speedup you gain
with just translating into C++ might not be huge.
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::export]]
List betahat(double t0, NumericMatrix Xr, NumericMatrix yr,
NumericVector tr, double h, int m) {
int n = Xr.nrow(), p = Xr.ncol();
arma::mat X(Xr.begin(), n, p, false);
arma::mat y(yr.begin(), n, 1, false);
arma::colvec t(tr.begin(), tr.size(), false);
arma::mat T = X;
T.each_col() %= (t-t0)/h;
arma::mat D = arma::join_rows(X,T);
arma::vec tp = t-t0;
arma::uvec q1 = arma::find(arma::abs(tp)<h);
arma::uvec q2 = arma::find(arma::abs(tp)>=h);
arma::vec kk = t;
kk.elem(q1) = ((1-arma::pow(tp.elem(q1)/h,2))/h)*0.75;
kk.elem(q2).zeros();
arma::mat W = (arma::diagmat(kk))/m;
arma::mat Inv = arma::trans(D)*W*D;
arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y;
arma::colvec betahat0(betahat.begin(),betahat.size()/2,false);
return List::create(Named("betahat") = betahat0);
}
Best wishes!
Honglang Wang
Office C402 Wells Hall
Department of Statistics and Probability
Michigan State University
1579 I Spartan Village, East Lansing, MI 48823
wangh...@msu.edu <mailto:wangh...@msu.edu>
On Tue, Dec 4, 2012 at 2:19 PM, Romain Francois
<rom...@r-enthusiasts.com <mailto:rom...@r-enthusiasts.com>> wrote:
You are calling back to R using Function. This is expensive.
What is ker ? Can you implement it in C++. This is a wild guess, but
that is I think where the bottleneck is.
Romain
Le 04/12/12 20:14, Honglang Wang a écrit :
Dear All,
I have tried out the first example by using RcppArmadillo, but I
am not
sure whether the code is efficient or not. And I did the
comparison of
the computation time.
1) R code using for loop in R: 87.22s
2) R code using apply: 77.86s
3) RcppArmadillo by using for loop in C++: 53.102s
4) RcppArmadillo together with apply in R: 47.310s
It is kind of not so big increase. I am wondering whether I used an
inefficient way for the C++ coding:
// [[Rcpp::depends(RcppArmadillo)__]]
#include <RcppArmadillo.h>
using namespace Rcpp;
// [[Rcpp::export]]
List betahat(Function ker, double t0, NumericMatrix Xr,
NumericMatrix
yr, NumericVector tr, double h, int m) {
int n = Xr.nrow(), p = Xr.ncol();
arma::mat X(Xr.begin(), n, p, false);
arma::mat y(yr.begin(), n, 1, false);
arma::colvec t(tr.begin(), tr.size(), false);
arma::mat T = X;
T.each_col() %= (t-t0)/h;
arma::mat D = arma::join_rows(X,T);
arma::vec kk =as<arma::vec>(ker(tr-t0,h));
arma::mat W = (arma::diagmat(kk))/m;
arma::mat Inv = arma::trans(D)*W*D;
arma::vec betahat = arma::inv(Inv)*arma::trans(D)*__W*y;
arma::colvec betahat0(betahat.begin(),__betahat.size()/2,false);
return List::create(Named("betahat") = betahat0);
}
Anyone could help me with how to increase the efficiency of the
coding?
Thanks.
Best wishes!
Honglang Wang
Office C402 Wells Hall
Department of Statistics and Probability
Michigan State University
1579 I Spartan Village, East Lansing, MI 48823
wangh...@msu.edu <mailto:wangh...@msu.edu>
<mailto:wangh...@msu.edu <mailto:wangh...@msu.edu>>
--
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30 <tel:%2B33%280%29%206%2028%2091%2030%2030>
R Graph Gallery: http://gallery.r-enthusiasts.__com
<http://gallery.r-enthusiasts.com>
`- http://bit.ly/SweN1Z : SuperStorm Sandy
blog: http://romainfrancois.blog.__free.fr
<http://romainfrancois.blog.free.fr>
|- http://bit.ly/RE6sYH : OOP with Rcpp modules
`- http://bit.ly/Thw7IK : Rcpp modules more flexible
_________________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-__project.org
<mailto:Rcpp-devel@lists.r-forge.r-project.org>
https://lists.r-forge.r-__project.org/cgi-bin/mailman/__listinfo/rcpp-devel
<https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel>
--
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30
R Graph Gallery: http://gallery.r-enthusiasts.com
`- http://bit.ly/SweN1Z : SuperStorm Sandy
blog: http://romainfrancois.blog.free.fr
|- http://bit.ly/RE6sYH : OOP with Rcpp modules
`- http://bit.ly/Thw7IK : Rcpp modules more flexible
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
