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?