Re: [Rcpp-devel] [RcppArmadillo] Result of Rcpp Wrap() for Sparse Matrix

[Douglas Bates]

On Wed, Jun 14, 2017 at 9:06 AM Serguei Sokol 
wrote:

> Le 14/06/2017 à 15:21, Douglas Bates a écrit :
> >
> >
> > On Wed, Jun 14, 2017 at 3:59 AM Serguei Sokol  <mailto:serguei.so...@gmail.com>> wrote:
> >
> > Le 13/06/2017 à 18:24, Douglas Bates a écrit :
> >  > On Tue, Jun 13, 2017 at 10:56 AM Binxiang Ni <
> binxian...@gmail.com <mailto:binxian...@gmail.com>  binxian...@gmail.com
> > <mailto:binxian...@gmail.com>>> wrote:
> >  >
> >  > Hi,
> >  >
> >  > I am working on fixing sparse matrix conversion for
> RcppArmadillo. Now a problem comes up to me: what kind of sparse matrix is
> expected to pass from
> >  > Armadillo to R? That is, what should the result of wrap() be?
> dgCMatrix(if logical, lgCMatrix or ngCMatrix)  or their original type?
> >  >
> >  >
> >  > What do you mean by "their original type"?
> >  >
> >  > It seems that the correspondence is
> >  > Armadillo   Matrix package
> >  > sp_mat   <=> dgCMatrix
> >  > sp_cx_mat <=> zgCMatrix
> >  > sp_imat  <=> igCMatrix
> > I would also consider the format used in a package slam.
> > It simply stores the indexes and non-zero values in a triplet
> (i,j,v).
> >
> >
> > That is the format of the dgTMatrix class from the Matrix package for R
> but not, as far as I can tell, in Armadillo.  A brief glance at the
> Armadillo
> > documentation indicates that sparse matrices are always in the
> compressed sparse column (CSC) format.
> Indeed, but nothing prevents Binxiang to develop a wrap() that will convert
> armadillo format to one or many of R formats, right?
>

Why? Is there a reason for doing type conversion from the dgCMatrix format
to another format in an Rcpp wrap function instead of with the existing
functions from the Matrix package?

Bear in mind that dgCMatrix is an efficient format both in terms of the
amount of memory required  (that's the "compressed" part of the name) and
in terms of performing operations with the matrix.  Most operations on
sparse matrices stored in the triplet format start by creating a CSC or CSR
(compressed sparse row) form of the matrix anyway.
___
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

Re: [Rcpp-devel] [RcppArmadillo] Result of Rcpp Wrap() for Sparse Matrix

[Douglas Bates]

On Wed, Jun 14, 2017 at 3:59 AM Serguei Sokol 
wrote:

> Le 13/06/2017 à 18:24, Douglas Bates a écrit :
> > On Tue, Jun 13, 2017 at 10:56 AM Binxiang Ni  <mailto:binxian...@gmail.com>> wrote:
> >
> > Hi,
> >
> > I am working on fixing sparse matrix conversion for RcppArmadillo.
> Now a problem comes up to me: what kind of sparse matrix is expected to
> pass from
> > Armadillo to R? That is, what should the result of wrap() be?
> dgCMatrix(if logical, lgCMatrix or ngCMatrix)  or their original type?
> >
> >
> > What do you mean by "their original type"?
> >
> > It seems that the correspondence is
> > Armadillo   Matrix package
> > sp_mat   <=> dgCMatrix
> > sp_cx_mat <=> zgCMatrix
> > sp_imat  <=> igCMatrix
> I would also consider the format used in a package slam.
> It simply stores the indexes and non-zero values in a triplet (i,j,v).
>

That is the format of the dgTMatrix class from the Matrix package for R but
not, as far as I can tell, in Armadillo.  A brief glance at the Armadillo
documentation indicates that sparse matrices are always in the compressed
sparse column (CSC) format.

I would point out that the sparse matrix facilities in Eigen and RcppEigen
are much more extensive than those in Armadillo.
___
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

Re: [Rcpp-devel] [RcppArmadillo] Result of Rcpp Wrap() for Sparse Matrix

[Douglas Bates]

On Tue, Jun 13, 2017 at 10:56 AM Binxiang Ni  wrote:

> Hi,
>
> I am working on fixing sparse matrix conversion for RcppArmadillo. Now a
> problem comes up to me: what kind of sparse matrix is expected to pass from
> Armadillo to R? That is, what should the result of wrap() be? dgCMatrix(if
> logical, lgCMatrix or ngCMatrix)  or their original type?
>

What do you mean by "their original type"?

It seems that the correspondence is
Armadillo   Matrix package
sp_mat   <=> dgCMatrix
sp_cx_mat <=> zgCMatrix
sp_imat  <=> igCMatrix

After that you are going to need to convert the values vector in the
Armadillo representation to a datatype available in R.  (In fact there is
always an implicit conversion because the rowind and colptr vectors are
unsigned integers in Armadillo and integers in R.  However that change is
just a cast of a pointer.)

I wouldn't worry about the logical values because LOGICAL in R is just a
32-bit integer.  ngCMatrix in R is a pattern matrix which only stores the
positions of the non-zeros, not their values.  These are used in the
symbolic phase of many sparse matrix operations, particularly
decompositions.  The symbolic phase only uses the pattern of nonzeros.  It
doesn't appear that there is a corresponding type in Armadillo.
___
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

Re: [Rcpp-devel] symmatu does not give symmetric matrix

[Douglas Bates]

First, this is not a good way of constructing a matrix to test the Cholesky
decomposition because the result is rank-deficient. It is a 2 by 2 matrix
of rank 1.  In fact it is the very definition of a rank-1 matrix - a column
vector multiplied by its transpose.

I am surprised at the result because I can't think of a way in computer
arithmetic that the off-diagonals would be different.  They should both be

prod(x)

The dput function will print the values out to full precision so try

dput(prod(x))
dput(x %*% t(x))

to see the full precision results.

Thirdly, R has two functions, crossprod and tcrossprod, to compute X'X and
XX', respectively.  They compute one triangle of the result and copy it
into the other, so the result is guaranteed to be symmetric. Try

dput(tcrossprod(x))

Finally, what are you going to do with the Cholesky decomposition?  If you
are doing this simply to get the decomposition, compiling C++ code, even
with all the help provided by Rcpp, is doing things the hard way. Just use
the chol function in R or the cholesky decomposition in the Matrix package
if you want more control.

On Wed, Jul 27, 2016 at 7:40 AM Dirk Eddelbuettel  wrote:

>
> Sorry for the delay in responding.  Mail delivery was stuck at R-Forge for
> a
> few days it seems.
>
> On 21 July 2016 at 10:46, Jendoubi, Takoua wrote:
> | Dear all,
> |
> | I am using RcppArmadillo to deal with some matrix computations.
> Specifically I
> | need to cholesky factorization of some symmetric matrices.
> |
> | I am generating random vectors using Rcpp and using them to construct
> symmetric
> | matrices.
> |
> | I have an error stating that my matrix is not symmetric although it
> definitely
> | should be. Here is the example I am working with (in R):
> |
> | >x
> | [1] -1.6683320 -0.8597148
> | >x%*%t(x)
> |  [,1]  [,2]
> | [1,] 2.783332 1.4342896
> | [2,] 1.434290 0.7391095
> |
> | Apparently, it is a rounding-off error. Is there any way to ensure that
> x%*%t
> | (x) gives an exactly symmetric matrix to use for cholesky factorization?
>
> I don't know of an automatic way.  R surely has no native type for
> symmetric
> matrices (at least not at the SEXP level as it really only has vectors with
> dim attributes).  I would probably start with some helper functions at the
> R
> or C++ level.  Some may have better tricks...
>
> Dirk
>
>
> --
> http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org
> ___
> 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
>
___
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

Re: [Rcpp-devel] RcppEigen: Windows binary from CRAN crashes R, but not when installing from source.

[Douglas Bates]

On Thu, Oct 16, 2014 at 11:35 AM, Kevin Ushey  wrote:

> I think John's advice is spot on here. The issue is only seen when
> `NDEBUG` is not defined.
>
> I can reproduce the crash (assertion failure) by ensuring I have
>
> CXXFLAGS=-UNDEBUG
>
> in my ~/.R/Makevars. Note that:
>
> 1. An assertion failure from Eigen implies you are doing something
> that you should not be doing, and
> 2. R by default sets -DNDEBUG whenever compiling by default, so I am
> surprised that you are not seeing it (does your package have a custom
> Makefile / Makevars or something to that effect? Or do you have your
> own custom Makevars somewhere?)
>
> Anyway, let's assume Eigen is right -- this means you're multiplying
> two non-conforming matrices, and hence your matrix product is
> undefined. Which makes sense, since you're now trying to multiply two
> non-conforming matrices. And if you want a scalar * matrix
> multiplication then you need to be using a different function.
>
> Note that this is exactly what Eigen's assertion was telling you here!
>

Which brings up the issue of the default Makefile for R packages overriding
other settings and enforcing -DNDEBUG.  I know that I would prefer to learn
that there were identifiable inconsistencies in the underlying code, even
if it meant that the failed assertion caused R to crash.  Those assertions
are there for a reason.  The decision to enforce a "fail faulty" scheme,
where known errors are silently allowed to go undetected, seems peculiar to
me.


>
> On Thu, Oct 16, 2014 at 8:58 AM, John Buonagurio
>  wrote:
> > Hi Henrik,
> >
> > You can just add #define EIGEN_NO_DEBUG to make it clear that you want
> to specifically disable Eigen assertions, without relying on what sets or
> does not set the DNDEBUG flag. It's in the Eigen documentation:
> >
> > http://eigen.tuxfamily.org/dox/TopicPreprocessorDirectives.html
> >
> > John
> >
> >> -Original Message-
> >> From: rcpp-devel-boun...@lists.r-forge.r-project.org [mailto:
> rcpp-devel-
> >> boun...@lists.r-forge.r-project.org] On Behalf Of Henrik Singmann
> >> Sent: Thursday, October 16, 2014 10:39 AM
> >> To: Dirk Eddelbuettel
> >> Cc: rcpp-de...@r-forge.wu-wien.ac.at
> >> Subject: Re: [Rcpp-devel] RcppEigen: Windows binary from CRAN crashes
> R, but
> >> not when installing from source.
> >>
> >> 
> >> Hi Dirk,
> >>
> >> I am sorry to address this again, as I am in no position to argue, but
> the problem
> >> is still the following:
> >> The code runs perfectly as your example showed but will fail when using
> either
> >> devtools or winbuilder or CRAN when *inside a package*. The problem of
> >> producing a minimally reproducible example of your liking is that
> either of those
> >> conditions cannot be met: I need the package to reproduce it.
> >>
> >> I also agree that it is an issue of multiplying non-conformable
> vectors. But as can
> >> be easily shown using the code you had sent, when compiling it on owns
> own
> >> (i.e., not using the conditions reproducing the crash described above)
> this only
> >> leads to a result of 0. However, when using the binary from either of
> the above
> >> it crashes R with an assertion error.
> >> It is also important to note that this did not happen in the past.
> >>
> >> So I probably just need to set the correct compiler flags or disable
> DNDEBUG or
> >> such a thing, but I did not manage to do this in such a way to prohibit
> the crash
> >> under the above described circumstances.
> >> I would really like to receive any advice on how to avoid this crash
> (assertion
> >> error) when using the binary compiled on CRAN, this is in the end the
> critical
> >> issue. *How can I disable the assertion error compiling in my code?*
> Just adding
> >> "#undef NDEBUG" at the beginning didn't work.
> >>
> >>
> >> On an unrelated note, the issue of within package or outside of a
> package also
> >> concerns the question of "using Eigen::..." versus prefacing all calls
> to Eigen
> >> parts with Eigen:: directly.
> >> While the code you had sent works great when using sourceCpp(), I didn't
> >> manage to get it to work in a package (even after wildly using
> >> compileAttributes). I had to replace all calls of e.g., VectorXd with
> >> Eigen::VectorXd. Is there a trick of how to do this inside a package?
> >>
> >> Btw, I agree that using the RcppAttributes is great. I hadn't used it,
> because, you
> >> know, "never touch a running system." But as it failed now, it is
> perhaps time for
> >> a change.
> >>
> >> Thanks again,
> >> Henrik
> >>
> >>
> >> Am 16.10.2014 um 16:04 schrieb Dirk Eddelbuettel:
> >> > On 16 October 2014 at 08:35, Dirk Eddelbuettel wrote:
> >> > |
> >> > |
> >> > | On 16 October 2014 at 15:07, Henrik Singmann wrote:
> >> > | | Hi Dirk and Kevin,
> >> > | |
> >> > | | I have now rebuild the package using the code Dirk send me (i.e.,
> using
> >> attributes) and the code still reliably crashes my R on Linux when
> using devtools
> >> (independent of RStudio), but not 

[Rcpp-devel] Yet another instance of "function 'dataptr' not provided ..."

[Douglas Bates]

I must have been away from writing R/Rcpp code for too long.

I started off trying to reproduce a calculation that is, literally, a
one-liner in Julia.  See

http://nbviewer.ipython.org/gist/dmbates/9746197

Now admittedly the calculation of the sums of the n choose k possible
subsets of size k from a vector or length n is aided by the fact that there
is a combinations iterator in Julia.  Nonetheless it is pretty amazing that
this can be written as

combsums(v::Vector, k) = [sum(c) for c in combinations(v,k)]

using this iterator and a comprehension.

I tried to write the combinations iterator in C++ but eventually tied
myself in knots so I decided to back off and write a function that does the
equivalent of the R function sample() and use that to generate a random
sample from the distribution of sums.

I have the C++ code working using Rcpp attributes but now I must generate a
package.  I have to be missing something obvious because my attempt at

http://github.com/dmbates/randomizationTest

produces the dread "function 'dataptr' not provided ..." message when I try
to invoke the R function randomsums.

Which of the many vignettes or manuals should I start reading?

(I can't believe I have spent two days going through innumerable
contortions to try to achieve the effect of a one-liner.)
___
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

Re: [Rcpp-devel] RcppEigen: Avoiding accessing elements with coeff(i, j) in sparse matrix

[Douglas Bates]

Two other small points.  I am assuming that the set of indices in the
second argument is sorted.  Also, it would be best to declare the function
as

bool foo(const MSpMat X, const Eigen::Map idx) as in the
enclosed


On Tue, Mar 25, 2014 at 10:45 AM, Douglas Bates  wrote:

> The enclosed works on this example.  The logic is to check each column in
> the index set for all the elements of the index set, except the one on the
> diagonal, being in the set of row indices. I have added some diagnostic
> output to demonstrate the flow.
>
>
>
>
> On Tue, Mar 25, 2014 at 9:48 AM, Douglas Bates wrote:
>
>> On Tue, Mar 25, 2014 at 6:42 AM, Søren Højsgaard wrote:
>>
>>> Dear all,
>>> I have a large sparse adjacency matrix X for an undirected graph. For a
>>> subset 'idx' of the vertices I want to find out if the subgraph defined by
>>> this subset is complete (i.e. has edges between all variables). So in R,
>>> one would check if X[idx,idx] has only ones outside the diagonal, but this
>>> is slow and I want to do it with Rcpp. Here is one attempt where I simply
>>> add up the elements above (or below) the diagonal, and to access the
>>> elements of X I use coeff which is relatively slow (because X is a sparse
>>> matrix).
>>>
>>> #include 
>>> typedef Eigen::MappedSparseMatrix MSpMat;
>>> typedef Eigen::SparseVector SpVec;
>>> typedef SpVec::InnerIterator InIter;
>>>
>>> //[[Rcpp::export]]
>>> double foo (MSpMat X, Eigen::VectorXi idx){
>>>   SpVec sidx = idx.sparseView();
>>>   double out=0;
>>>   int i2, j2;
>>>   for (InIter i_(sidx); i_; ++i_){
>>> i2 = i_.index();
>>> for (InIter j_(sidx); j_; ++j_){
>>>   j2 = j_.index();
>>>   if (i2>j2)
>>> out += X.coeff( i2, j2);
>>> }
>>>   }
>>>   return out;
>>> }
>>>
>>> /*** R
>>> library(Matrix)
>>> M1 <- new("dgCMatrix"
>>> , i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L,
>>> 5L, 3L, 4L)
>>> , p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
>>> , Dim = c(6L, 6L)
>>> , Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c",
>>> "d", "e", "f"))
>>> , x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
>>> , factors = list()
>>> )
>>> M1
>>>
>>> foo(M1, c(2,3,4))
>>>
>>> */
>>> Can anyone see a better way of doing this?
>>>
>>> I was thinking about whether the sparse matrix iterators could be a
>>> solution, but I can't a hold on it.
>>
>>
>> Sparse matrix inner and outer iterators are definitely the way to go.
>>  I'll write a version using them.
>>
>>
>>> I was also thinking about using sparse matrices in RcppArmadillo, but I
>>> guess the problem (speed) is the same (haven't tried!).
>>>
>>
>>
>>> Any advice will be greatly appreciated.
>>> Cheers
>>> Søren
>>>
>>>
>>>
>>> ___
>>> 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
>>>
>>>
>
// [[Rcpp::depends(RcppEigen)]]
#include 
typedef Eigen::MappedSparseMatrix MSpMat;
typedef MSpMat::InnerIterator InIter;


//[[Rcpp::export]]
bool foo (const MSpMat X, const Eigen::Map idx){
int n = X.cols();
if (X.rows() != n) throw std::invalid_argument("Sparse matrix X must be square");
for (int i = 0; i < idx.size(); ++i) {
int jj = idx[i] - 1;// 0-based column index
InIter it(X,jj);
Rcpp::Rcout << "i = " << i << ", jj = " << jj << std::endl;
for (int k = 0; k < idx.size(); ++k) {
int kk = idx[k]-1;  // 0-based row index to search for
if (kk == jj) continue;  // skip positions on the diagonal
Rcpp::Rcout << "  kk = " << kk << ", it.row =";
bool foundit = false;
for (; it; ++it) {
Rcpp::Rcout << " " << it.row();
if (it.row() == kk) {
foundit = true;
++it;
break;
}
if (it.row() > kk) return false;
}
if (!foundit) return false;
Rcpp::Rcout << std::endl;
}
}
return true;
}

/*** R
library(Matrix)
M1 <- new("dgCMatrix"
, i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L, 5L, 3L, 4L)
, p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
, Dim = c(6L, 6L)
, Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c", "d", "e", "f"))
, x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
, factors = list()
)
M1

foo(M1, c(2L,3L,4L))
***/
___
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

Re: [Rcpp-devel] RcppEigen: Avoiding accessing elements with coeff(i, j) in sparse matrix

[Douglas Bates]

The enclosed works on this example.  The logic is to check each column in
the index set for all the elements of the index set, except the one on the
diagonal, being in the set of row indices. I have added some diagnostic
output to demonstrate the flow.




On Tue, Mar 25, 2014 at 9:48 AM, Douglas Bates  wrote:

> On Tue, Mar 25, 2014 at 6:42 AM, Søren Højsgaard wrote:
>
>> Dear all,
>> I have a large sparse adjacency matrix X for an undirected graph. For a
>> subset 'idx' of the vertices I want to find out if the subgraph defined by
>> this subset is complete (i.e. has edges between all variables). So in R,
>> one would check if X[idx,idx] has only ones outside the diagonal, but this
>> is slow and I want to do it with Rcpp. Here is one attempt where I simply
>> add up the elements above (or below) the diagonal, and to access the
>> elements of X I use coeff which is relatively slow (because X is a sparse
>> matrix).
>>
>> #include 
>> typedef Eigen::MappedSparseMatrix MSpMat;
>> typedef Eigen::SparseVector SpVec;
>> typedef SpVec::InnerIterator InIter;
>>
>> //[[Rcpp::export]]
>> double foo (MSpMat X, Eigen::VectorXi idx){
>>   SpVec sidx = idx.sparseView();
>>   double out=0;
>>   int i2, j2;
>>   for (InIter i_(sidx); i_; ++i_){
>> i2 = i_.index();
>> for (InIter j_(sidx); j_; ++j_){
>>   j2 = j_.index();
>>   if (i2>j2)
>> out += X.coeff( i2, j2);
>> }
>>   }
>>   return out;
>> }
>>
>> /*** R
>> library(Matrix)
>> M1 <- new("dgCMatrix"
>> , i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L,
>> 5L, 3L, 4L)
>> , p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
>> , Dim = c(6L, 6L)
>> , Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c",
>> "d", "e", "f"))
>> , x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
>> , factors = list()
>> )
>> M1
>>
>> foo(M1, c(2,3,4))
>>
>> */
>> Can anyone see a better way of doing this?
>>
>> I was thinking about whether the sparse matrix iterators could be a
>> solution, but I can't a hold on it.
>
>
> Sparse matrix inner and outer iterators are definitely the way to go.
>  I'll write a version using them.
>
>
>> I was also thinking about using sparse matrices in RcppArmadillo, but I
>> guess the problem (speed) is the same (haven't tried!).
>>
>
>
>> Any advice will be greatly appreciated.
>> Cheers
>> Søren
>>
>>
>>
>> ___
>> 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
>>
>>
// [[Rcpp::depends(RcppEigen)]]
#include 
typedef Eigen::MappedSparseMatrix MSpMat;
typedef MSpMat::InnerIterator InIter;


//[[Rcpp::export]]
bool foo (MSpMat X, Eigen::VectorXi idx){
int n = X.cols();
if (X.rows() != n) throw std::invalid_argument("Sparse matrix X must be square");
for (int i = 0; i < idx.size(); ++i) {
int jj = idx[i] - 1;// 0-based column index
InIter it(X,jj);
Rcpp::Rcout << "i = " << i << ", jj = " << jj << std::endl;
for (int k = 0; k < idx.size(); ++k) {
int kk = idx[k]-1;  // 0-based row index to search for
if (kk == jj) continue;  // skip positions on the diagonal
Rcpp::Rcout << "  kk = " << kk << ", it.row =";
bool foundit = false;
for (; it; ++it) {
Rcpp::Rcout << " " << it.row();
if (it.row() == kk) {
foundit = true;
++it;
break;
}
if (it.row() > kk) return false;
}
if (!foundit) return false;
Rcpp::Rcout << std::endl;
}
}
return true;
}

/*** R
library(Matrix)
M1 <- new("dgCMatrix"
, i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L, 5L, 3L, 4L)
, p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
, Dim = c(6L, 6L)
, Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c", "d", "e", "f"))
, x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
, factors = list()
)
M1

foo(M1, c(2,3,4))
***/
___
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

Re: [Rcpp-devel] RcppEigen: Avoiding accessing elements with coeff(i, j) in sparse matrix

[Douglas Bates]

On Tue, Mar 25, 2014 at 6:42 AM, Søren Højsgaard  wrote:

> Dear all,
> I have a large sparse adjacency matrix X for an undirected graph. For a
> subset 'idx' of the vertices I want to find out if the subgraph defined by
> this subset is complete (i.e. has edges between all variables). So in R,
> one would check if X[idx,idx] has only ones outside the diagonal, but this
> is slow and I want to do it with Rcpp. Here is one attempt where I simply
> add up the elements above (or below) the diagonal, and to access the
> elements of X I use coeff which is relatively slow (because X is a sparse
> matrix).
>
> #include 
> typedef Eigen::MappedSparseMatrix MSpMat;
> typedef Eigen::SparseVector SpVec;
> typedef SpVec::InnerIterator InIter;
>
> //[[Rcpp::export]]
> double foo (MSpMat X, Eigen::VectorXi idx){
>   SpVec sidx = idx.sparseView();
>   double out=0;
>   int i2, j2;
>   for (InIter i_(sidx); i_; ++i_){
> i2 = i_.index();
> for (InIter j_(sidx); j_; ++j_){
>   j2 = j_.index();
>   if (i2>j2)
> out += X.coeff( i2, j2);
> }
>   }
>   return out;
> }
>
> /*** R
> library(Matrix)
> M1 <- new("dgCMatrix"
> , i = c(1L, 2L, 3L, 0L, 2L, 3L, 0L, 1L, 3L, 0L, 1L, 2L, 4L, 5L, 3L,
> 5L, 3L, 4L)
> , p = c(0L, 3L, 6L, 9L, 14L, 16L, 18L)
> , Dim = c(6L, 6L)
> , Dimnames = list(c("a", "b", "c", "d", "e", "f"), c("a", "b", "c",
> "d", "e", "f"))
> , x = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
> , factors = list()
> )
> M1
>
> foo(M1, c(2,3,4))
>
> */
> Can anyone see a better way of doing this?
>
> I was thinking about whether the sparse matrix iterators could be a
> solution, but I can't a hold on it.


Sparse matrix inner and outer iterators are definitely the way to go.  I'll
write a version using them.


> I was also thinking about using sparse matrices in RcppArmadillo, but I
> guess the problem (speed) is the same (haven't tried!).
>


> Any advice will be greatly appreciated.
> Cheers
> Søren
>
>
>
> ___
> 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
>
>
___
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

Re: [Rcpp-devel] Sparse matrix operations

[Douglas Bates]

On Mon, Feb 3, 2014 at 11:48 AM, Dirk Eddelbuettel  wrote:

>
> On 3 February 2014 at 10:23, French, Joshua wrote:
> | Soren and Doug,
> |
> | Thanks for the info about the sparse matrices.  I'll give it a go and
> seen what
> | happens.  If it works out, then perhaps I'll be able to provide a nice
> | RcppGallery example.
>
> Yes -- there isn't much yet in RcppArmadillo but this is expected to
> grow. And CRAN now has rARPACK so we get cheaply to the object code too.
>
> One last comment I should have made earlier: "dense" data from R comes over
> cheaply as a SEXP; "sparse" will always require a copy.  Keep that in mind,
> and on the margin don't believe anything any of us say but keep profiling
> and
> measuring :)
>

Not really.  RcppEigen offers conversion of SEXPs to Eigen::SparseMatrix
and Eigen::MappedSparseMatrix types.  The second uses the storage from R
without copying and hence should generally be declared with the const
modifier.  I think that wonderful vignetter by Bates and Eddelbuettel
describes this :-)
___
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

Re: [Rcpp-devel] Sparse matrix operations

[Douglas Bates]

On Fri, Jan 31, 2014 at 1:14 PM, French, Joshua
wrote:

> Hello everyone,
>
> For those of you who have used the sparse matrix capabilities of
> Armadillo/RcppArmadillo, how seamless are the operations in moving between
> sparse and dense matrices?  I know there are some tricks for getting
> sparseMatrix objects into SpMat objects (Dirk has a nice Rcpp Gallery post
> on this), but is it just as easy when it comes to doing Linear Algebra? My
> google ninja skills didn't seem to produce any helpful information.
>
> E.g., if A is a dense matrix and B is a sparse matrix in Armadillo, it is
> pretty much seamless to do calculations like A + B, A*B, solve(B, A), etc.?
>
>

I'm not that familiar with the sparse matrix capabilities of Armadillo but
I do know those of Eigen fairly well.  Operations on sparse matrices
require more care and thought than do those on dense matrices.  It is
possible to use "one size fits all" algorithms but they can be slow or
inaccurate compared to more specialized approaches to operations like
solve(B,A).  For dense matrices you could use an LU factorization or a QR
factorization on a symmetric B and never notice the difference relative to
using a Cholesky factorization.  With sparse matrices you could notice the
difference.

Operations like A + B would not be terribly meaningful because the result
would be dense if A is dense so you may as well extract a dense version of
B and use that.  A * B is definitely available in Eigen and I imagine also
available with Armadillo.  A general solve(B, A) could be written using a
sparse LU but that may not be an effective way of solving such a system.
 Most sparse matrix operations have a symbolic phase and a numeric phase.
 In the symbolic phase the problem is analyzed to determine the positions
of the nonzeros in the result and possibly auxiliary information such as a
fill-reducing permutation.  Often if you need to solve several systems of
the same form but with different numerical values you can save and update
the factorization rather than starting from scratch each time.

This is why I prefer Eigen which has classes representing factorizations in
addition to the one-shot methods.

>
> I'm trying to decide whether it would be worth it to convert my code from
> R using sparseMatrix objects (Matrix package) to a C++ equivalent using
> RcppArmadillo.  On the same topic, would it be easier/more difficult if I
> used RcppEigen instead?
>
> Thanks,
>
> Joshua
> --
> Joshua French, Ph.D.
> Assistant Professor
> Department of Mathematical and Statistical Sciences
> University of Colorado Denver
> joshua.fre...@ucdenver.edu
> http://math.ucdenver.edu/~jfrench/
> Ph:  303-315-1709  Fax:  303-315-1701
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] How many copies of my object do I make?

[Douglas Bates]

Typo

On Tue, Dec 10, 2013 at 11:48 AM, Douglas Bates  wrote:

> You can avoid one copy operation by declaring X as
>
>  const MSpMat X(as(xx_));
>

xx_ should be XX_.
___
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

Re: [Rcpp-devel] How many copies of my object do I make?

[Douglas Bates]

You can avoid one copy operation by declaring X as

 const MSpMat X(as(xx_));

As it is you are allocating storage for a copy of the sparse matrix which
will not be modified and hence does not need to be copied.


On Tue, Dec 10, 2013 at 11:33 AM, Dirk Eddelbuettel  wrote:

>
> On 9 December 2013 at 22:14, Søren Højsgaard wrote:
> | Dear all,
> | Using RcppEigen I've created a function for converting a sparse matrix
> (a dgCMatrix) to a standard dense matrix:
> |
> | typedef Eigen::SparseMatrix SpMatd;
> | typedef Eigen::MappedSparseMatrix MSpMat;
> |
> | // [[Rcpp::export]]
> | SEXP do_dgCMatrix2matrix ( SEXP XX_ ){
> |   S4 DD(wrap(XX_));
> |   List dn = clone(List(DD.slot("Dimnames")));
> |   SpMatd X(as(XX_));
> |   Eigen::MatrixXd dMat;
> |   dMat = Eigen::MatrixXd( X );
> |   NumericMatrix Xout(wrap(dMat));
> |   Xout.attr("dimnames") = dn;
> |   return Xout;
> | }
> |
> | The function does what I expect it to, but I am curious about how many
> times I really create some sort of "copy" of my input matrix and when all I
> do is put a small layer of "vernish" over my input matrix?
> |
> | 1) dMat is dense matrix, and that must (I guess) refer to new "object"??
> |
> | 2) X is a sparse matrix, and I suppose it is a new object??
> |
> | 2') If I had declared it to be MSpMat, then it would really just be the
> input object with a few bells and whistles??
> |
> | 3) But how about DD and Xout? Are they new objects (in terms of memory
> usage)?
> |
> | Put differently, if my input matrix occupies B bytes of memory, how many
> times B have I then created along the way?
>
> Unsure. You really need to look.
>
> For sure when i) clone is called, ii) each time you pass from sparse to
> dense
> and back.  For the others, you need to dig deeper.
>
> Dirk
>
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> ___
> 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
>
>
___
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

[Rcpp-devel] Question on lme4 book

[Douglas Bates]

Yesterday Taylor Russ asked

What it the proper citation for the lme4 package and the Bates' book?

 Also, can lme4 datasets (e.g., Pastes, ScotsSec, InstEval etc.) be
 used for illustration in publications?  Can the authors grant
 permission or is the permission from the source needed?

 Many thanks for the package and the book.  When can I hold a
non-digital copy in my hands?

I inadvertently deleted the message and so must respond without maintaining
the thread.

The data sets can be used in other publications.  At least my understanding
is that the data themselves cannot be copyright (despite the "Microsoft
Patents 1's, 0's" headline in The Onion many years ago - for those of you
who don't know that The Onion is a satirical newspaper, that didn't really
occur).  It is only the representation of the data, such as a table in a
copyright publication, that can be copyright.  I suppose I should provide
the usual caveat, "I am (thankfully) not a lawyer".

The other lme4 authors may be able to respond to the question of citing the
lme4 package.  I regret to say that I don't know of a good way of citing
the book and that there won't be non-digital copies.

Partly this can be attributed to my personality - I'm good at starting
projects but not so good at finishing them.  However, finishing the book
would involve spending time maintaining and developing the lme4 package for
CRAN and I have completely lost my enthusiasm for doing so.

As many of you know, I am doing most of my work in the Julia language (
www.julialang.org) now.  R is wonderful and I enjoyed most of my time
working on R and R packages but there are inherent limitations to R,
particularly when trying to achieve good performance on fitting complex
models to large data sets, that make this difficult.  It would be
attractive to have a "pure R" implementation of mixed-models but I don't
see a way of making it run quickly and without using a lot of memory.  In
Julia I can build a package that achieves good performance without the need
to interface to code written in C, C++ or Fortran - in the sense that my
package doesn't need to require compilation of code outside of that
provided by the language itself.

It is not surprising that the design of R is starting to show its age.
 Although R has only been around for 15-18 years, its syntax and much of
the semantics are based on the design of "S3" which is 25-30 years old.

R packages can include code to be compiled along with the interface code
and there are many wonderful tools to facilitate this - such as the Rcpp
package, the devtools package and RStudio support for these packages.  I
used these in the compiled code underlying lme4_1.0.

But even though Dirk would describe the use of Rcpp as "seamless", in my
experience it is not, especially if you wish to have your package available
on CRAN.

Maintaining an Rcpp-based package on CRAN these days is a case of "no good
deed shall go unpunished" and "the flogging will continue until morale
improves".  I am the maintainer of the RcppEigen package which apparently
also makes me the maintainer of an Eigen port to Solaris.  When compilers
on Solaris report errors in code from Eigen I am supposed to fix them.
 This is difficult in that I don't have access to any computers running
Solaris, which is a proprietary operating system as far as I can tell, and
Eigen is a complex code base using what is called "template
meta-programming" in C++.  Making modifications to such code can be
difficult.  I can't claim to fully understand all the details in Eigen and
in Rcpp.  I am a user of these code bases, not a developer. The Eigen
authors themselves don't test their code under Solaris because they don't
have access to Solaris systems either and they don't regard Solaris as an
important platform for numerical computing.  The CRAN maintainers feel
differently, which puts me in a box.

There are days when I am tempted to say, "okay, if RcppEigen is not
suitable for CRAN then remove it" which would result in removal of all the
packages that depend on it, including lme4.  That may seem childish of me
but I really don't know what else to do.

So I have reached the point of saying "goodbye" to R, Rcpp and lme4 and
switching all of my development effort to Julia.  I'm sorry but others are
going to need to determine how to maintain lme4 to the satisfaction of the
CRAN maintainers or whether there should be an alternative distribution
mechanism for R packages.
___
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

Re: [Rcpp-devel] RcppEigen needs 1GB memory to compile

[Douglas Bates]

On Mon, Dec 2, 2013 at 2:20 PM, Dirk Eddelbuettel  wrote:

>
> On 2 December 2013 at 13:52, Douglas Bates wrote:
> | The important thing is to use the -g0 flag.  Even though RcppEigen is  a
>
> Right.
>
> And for the opencpu server deployment, you may want to edit the -g out
> /etc/R/Makeconf as well. As I recall, there was an r-devel thread in which
> the desire to override such 'earlier' settings in what R CMD ... uses via
> later settings, and as I recall Simon stated that it couldn not be done.
>
> | header-only package we include an example R function fastlm. If you
> leave the
> | symbols in the DLL file you get a massive library size whereas stripping
> the
> | symbols provides you with a much smaller file size.  And because
> packages that
> | use RcppEigen only use the headers, not the DLL file, it doesn't matter
> if the
> | symbols are stripped.
>
> Do you think we should move fastLm out?
>

And face the "Wrath of Achim" for removing a function that is described in
a JSS paper?
___
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

Re: [Rcpp-devel] RcppEigen needs 1GB memory to compile

[Douglas Bates]

The important thing is to use the -g0 flag.  Even though RcppEigen is  a
header-only package we include an example R function fastlm. If you leave
the symbols in the DLL file you get a massive library size whereas
stripping the symbols provides you with a much smaller file size.  And
because packages that use RcppEigen only use the headers, not the DLL file,
it doesn't matter if the symbols are stripped.


On Mon, Dec 2, 2013 at 1:31 PM, Dirk Eddelbuettel  wrote:

>
> On 2 December 2013 at 11:17, Jeroen Ooms wrote:
> | I noticed that RcppEigen fails to install on my servers because it
> requires
> | more than 1gb of memory to compile. Is this expected? Are there any
> flags or
>
> Yes. It also falls over on R-Forge.
>
> | options I could set to trade of some memory for cpu? A short simulation:
>
> I have had to resort to similar tricks on some smaller architecture when
> autobuilding demanding Debian packages.  Try -O0 -g0 for a start.
>
> Or if you need a smile, this still holds:
>
>   http://dilbert.com/strips/comic/1995-06-24/
>
> Dirk
>
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> ___
> 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
>
___
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

[Rcpp-devel] Header-only version of Rcpp provides a couple of gotcha's in RcppEigen

[Douglas Bates]

If I install the latest, header-only, git version of Rcpp (f117a70), which
is a great achievement by the way,  installation of RcppEigen encounters a
few minor problems.  One seems to be related to conflicting definitions of
trunc in Rmath.h and in fstream

In file included from fastLm.cpp:23:
In file included from ./fastLm.h:25:
In file included from ../inst/include/RcppEigen.h:25:
In file included from ../inst/include/RcppEigenForward.h:39:
In file included from ../inst/include/unsupported/Eigen/SparseExtra:22:
/usr/lib/gcc/x86_64-linux-gnu/4.8/../../../../include/c++/4.8/fstream:641:60:
error: no member named 'Rf_ftrunc' in 'std::ios_base'
In file included from  ios_base::openmode __mode =
ios_base::out|ios_base::trunc)
   ~~^
RcppEigen.cpp:22:
In file included from ../inst/include/RcppEigen.h:25:
/home/bates/build/R-devel/include/Rmath.h:In file included from
../inst/include/RcppEigenForward.h348::3915::
In file included from  note:
../inst/include/unsupported/Eigen/SparseExtra:expanded 22:
from
/usr/lib/gcc/x86_64-linux-gnu/4.8/../../../../include/c++/4.8/fstreammacro
:641'trunc'
:60: error: #define trunc   ftrunc
no ^
member /home/bates/build/R-devel/include/Rmath.hnamed :247'Rf_ftrunc' :17in
: 'std::ios_base'
note: expanded from macro 'ftrunc'
#define ftrunc  Rf_ftrunc
^
 ios_base::openmode __mode =
ios_base::out|ios_base::trunc)

The other looks like I was playing a bit fast and loose with converting an
attribute to a list

make: *** Waiting for unfinished jobs
fastLm.cpp:217:20: error: call to constructor of 'List' (aka 'Vector<19>')
is ambiguous
List  dimnames(NumericMatrix(Xs).attr("dimnames"));
  ^~~

There is one more problem, but I think that needs a namespace qualifier.

fastLm.cpp:240:6: error: use of undeclared identifier
'forward_exception_to_r'; did you mean 'Rcpp::forward_exception_to_r'?
forward_exception_to_r( ex );
^~
Rcpp::forward_exception_to_r
/home/bates/R/x86_64-unknown-linux-gnu-library/3.1/Rcpp/include/Rcpp/exceptions.h:172:17:
note: 'Rcpp::forward_exception_to_r' declared here
inline void forward_exception_to_r( const std::exception& ex){
^

Do these all look like issues that I should address in RcppEigen?
___
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

[Rcpp-devel] New release of RcppEigen based on Eigen 3.2.0

[Douglas Bates]

A new release of RcppEigen based on Eigen 3.2.0 has been uploaded to CRAN.
___
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

Re: [Rcpp-devel] How much speedup for matrix operations?

[Douglas Bates]

Another thing you can do in compiled code that is difficult to do in R is
to avoid the allocate/free sequences when calling a code segment thousands
of times.  R's functional semantics (i.e. the expectation that functions do
not modify their arguments) is great - until you really want to overwrite a
value and you can't.  With Eigen and Armadillo and even with raw BLAS calls
you can pass the array that will be the result as well as the arguments.


On Wed, Nov 6, 2013 at 12:19 PM, Simon Zehnder  wrote:

> As Romain already mentioned, there could be rather a memory issue -
> copying is one issue. Very often time gets lost in loads and writes with
> cache misses. How large are your matrices? Have you done a memory-profiling?
>
> Best
> Simon
>
> On 06 Nov 2013, at 19:04, Xavier Robin  wrote:
>
> > On 11/6/13 6:38 PM, Romain Francois wrote:
> >> This very much depends on the code but there is a good chance that
> RcppArmadillo will generate code making less data copies, etc ...
> >>
> >> Hard to say without seeing the code.
> >>
> >> Romain
> >>
> > Most of the code (or at least the slow, highly repeated parts) look like:
> >
> >>A <- t(c + t(W) %*% X)
> >>E <- (exp(A) * (1 - 1 / A) + 1 / A) / (exp(A) - 1)
> >>E[abs(A) < sqrt(.Machine$double.eps) * 2 ] <- 0.5
> >>
> >>B <- t(b + W %*% t(E))
> >>X2 <- 1 / (1 + exp(-B))
> >>
> >>A2 <- t(c + t(W) %*% X2)
> >>E2 <- (exp(A2) * (1 - 1 / A2) + 1 / A2) / (exp(A2) - 1)
> >>E2[abs(A2) < sqrt(.Machine$double.eps) * 2 ] <- 0.5
> >>
> >>delta <- (t(X) %*% E - t(X2) %*% E2)
> >>W <- W + delta
> >
> > Where b and c are vectors, W and X matrices. All this is encapsulated in
> a function, that is called a few thousand times in a for loop, with some
> sanity checks. (But it didn't appear to have much impact on the speed... if
> I remove the matrix operations so it does nothing, it executes nearly
> instantly). I understand from Dirk and Douglas that it probably isn't going
> to make a huge difference, though (not by orders).
> >
> > Thanks,
> > Xavier
> >
> > --
> > Xavier Robin, PhD
> > Cellular Signal Integration Group (C-SIG) - http://www.lindinglab.org
> > Center for Biological Sequence Analysis (CBS) - http://www.cbs.dtu.dk
> > Department of Systems Biology - Technical University of Denmark (DTU)
> > Anker Engelundsvej, Building 301, DK-2800 Lyngby, DENMARK.
> >
> > ___
> > 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
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] How much speedup for matrix operations?

[Douglas Bates]

By default Eigen does not use BLAS, which can be good or bad, depending on
the situation.  I notice that the second largest total time is spent in
t.default which may mean that you are using an operation like

t(X) %*% X

If so, you can save yourself time by using the crossprod or tcrossprod
functions.  For example, the expression above is more cleanly written as

crossprod(X)

or, for t(X) %*% Y,

crossprod(X, Y)

Eigen or Armadillo could help you avoid making unnecessary copied but if
your calculation does end up being dominated by matrix multiplications you
can't expect to gain much speed relative to R.  You may want to check the
type of BLAS you are using.  For Intel processors MKL is generally the
fastest (but proprietary) with OpenBLAS in second place.

I see that I am, as often happens, giving you similar advice to Dirk's
response.


On Wed, Nov 6, 2013 at 11:35 AM, Xavier Robin  wrote:

> Hi,
>
> I have a pure-R code that spends most of the time performing vector and
> matrix operations, as shown by the summary of Rprof:
>
>>self.time self.pct total.time total.pct
>> "%*%" 903.2477.67 903.24 77.67
>> "t.default"76.26 6.56  76.26 6.56
>> "-"36.60 3.15  36.60 3.15
>> "+"24.44 2.10  24.44 2.10
>> "/"24.22 2.08  24.22 2.08
>> "exp"  20.26 1.74  20.26 1.74
>> "predict.myClass"  17.68 1.52 503.82 43.32
>> "*"11.90 1.02  11.90 1.02
>> "t" 9.38 0.81 811.94 69.82
>> "update.batch"  8.04 0.69 654.68 56.30
>> ...
>>
> So mostly matrix %*% matrix multiplications, transpositions, vector +-/*
> matrix operations and exponentiations, representing >95% of the computation
> time.
> I have very few loops and if/else blocks.
>
> I want to speed up this code, and I am considering reimplementing it (or
> part of it) with RcppEigen or RcppArmadillo.
>
> However, I read that both Eigen and Amarillo use the underlying BLAS, like
> R.
> My question is, can I expect any significant speed-up from an Rcpp
> re-implementation in this case, given it is already mostly matrix algebra
> (which are supposed to be pretty efficient in R)?
>
> Thanks,
> Xavier
>
> --
> Xavier Robin, PhD
> Cellular Signal Integration Group (C-SIG) - http://www.lindinglab.org
> Center for Biological Sequence Analysis (CBS) - http://www.cbs.dtu.dk
> Department of Systems Biology - Technical University of Denmark (DTU)
> Anker Engelundsvej, Building 301, DK-2800 Lyngby, DENMARK.
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] RcppEigen not compiling under Solaris 10 INTEL

[Douglas Bates]

On Mon, Sep 30, 2013 at 4:42 PM, Douglas Bates  wrote:

> On Mon, Sep 30, 2013 at 4:16 PM, Dirk Eddelbuettel  wrote:
>
>>
>> Hi Martyn,
>>
>> On 30 September 2013 at 21:02, Martyn Plummer wrote:
>> | I've already been asked to look into this by Brian.  I can't say I've
>> rolled up my sleeves very far, but the problem appears to be in Eigen
>> itself.
>>
>> Cool!  Good to know.
>>
>> | I should clarify that I'm not volunteering to maintain RcppEigen, but I
>> will provide patches as and when.
>>
>> That too is appreciated!
>>
>> The Eigen code on CRAN has not been updated in a while. The version number
>> (0.3.1.2.1) reflect both upstream Eigen 3.1.2 (two patch releases behind
>> Eigen in the 3.1.* series which is now at 3.1.4) and the mandatory change
>> for
>> JSS once the paper was accepted.
>>
>> Changing 3.1.2 to 3.1.3/3.1.4 may help. Or it may not.  Maybe 3.2.0
>> upstream
>> will help.
>>
>> Changing to 3.1.4 or, better, 3.2, would be a good idea.  In particular,
> the 3.2 version provides Cholesky, LU and QR sparse factorizations that do
> not depend on SuiteSparse (the package of C libraries of which CHOLMOD is
> part).
>
> Unfortunately, the current maintainer played a little fast and loose by
> adding a couple of wrinkles to the Eigen code itself and using these
> wrinkles in lme4.
>
> I had been hoping to put off the switch to Eigen 3.2 for a while longer as
> I am facing some deadline pressures on another project.  But perhaps it
> would be best to create an experimental branch using Eigen 3.2.  Best case,
> the fixes for lme4 are easy.  Worst case ...
>
> Ben's message came through while I was writing this.  If you are
> interested look at the last link he gives
>
> http://permalink.gmane.org/gmane.comp.lib.eigen/991
>
> and you will see that the sort of error messages showing up in the
> RcppEigen log were present in 2010 and the reaction of the Eigen developers
> is that the compiler is complaining about code that follows the standard
> and is compilable on other systems.
>

It seems that a good place for someone with access to a Solaris system to
start is to download versions 3.1.4 and 3.2 of Eigen from
http://eigen.tuxfamily.org and see if the simple program

$ cat test.cc
#include "eigen/Eigen/Eigen"
using namespace Eigen;
int main() {
 return 0;
}
$
$ CC test.cc

compiles with one of those versions.  If that doesn't compile the way
forward will be difficult.
___
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

Re: [Rcpp-devel] RcppEigen not compiling under Solaris 10 INTEL

[Douglas Bates]

On Mon, Sep 30, 2013 at 4:16 PM, Dirk Eddelbuettel  wrote:

>
> Hi Martyn,
>
> On 30 September 2013 at 21:02, Martyn Plummer wrote:
> | I've already been asked to look into this by Brian.  I can't say I've
> rolled up my sleeves very far, but the problem appears to be in Eigen
> itself.
>
> Cool!  Good to know.
>
> | I should clarify that I'm not volunteering to maintain RcppEigen, but I
> will provide patches as and when.
>
> That too is appreciated!
>
> The Eigen code on CRAN has not been updated in a while. The version number
> (0.3.1.2.1) reflect both upstream Eigen 3.1.2 (two patch releases behind
> Eigen in the 3.1.* series which is now at 3.1.4) and the mandatory change
> for
> JSS once the paper was accepted.
>
> Changing 3.1.2 to 3.1.3/3.1.4 may help. Or it may not.  Maybe 3.2.0
> upstream
> will help.
>
> Changing to 3.1.4 or, better, 3.2, would be a good idea.  In particular,
the 3.2 version provides Cholesky, LU and QR sparse factorizations that do
not depend on SuiteSparse (the package of C libraries of which CHOLMOD is
part).

Unfortunately, the current maintainer played a little fast and loose by
adding a couple of wrinkles to the Eigen code itself and using these
wrinkles in lme4.

I had been hoping to put off the switch to Eigen 3.2 for a while longer as
I am facing some deadline pressures on another project.  But perhaps it
would be best to create an experimental branch using Eigen 3.2.  Best case,
the fixes for lme4 are easy.  Worst case ...

Ben's message came through while I was writing this.  If you are interested
look at the last link he gives

http://permalink.gmane.org/gmane.comp.lib.eigen/991

and you will see that the sort of error messages showing up in the
RcppEigen log were present in 2010 and the reaction of the Eigen developers
is that the compiler is complaining about code that follows the standard
and is compilable on other systems.


> All this is now of somewhat greater importance as lme4 on CRAN uses
> RcppEigen.  "Just make it work and freeze it" is one approach, updating to
> current versions is another.
>
> As I said, it's an open job for someone who wants to earn some karma.
>
> But in any event, I am glad that you accepted that challenge.  Whichever
> timeline you can work, it'll be appreciated.  We can set you up with
> R-Forge
> commit privileges for this (and thereby the rest of Rcpp) as well.
>
> Cheers, Dirk
>
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> ___
> 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
>
___
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

Re: [Rcpp-devel] RcppEigen not compiling under Solaris 10 INTEL

[Douglas Bates]

Change "either x86 or sparse" to "either x86 or sparc".  I write sparse
more than I write sparc and my fingers took over.


On Mon, Sep 30, 2013 at 4:20 PM, Douglas Bates  wrote:

> On Mon, Sep 30, 2013 at 3:56 PM, Dirk Eddelbuettel  wrote:
>
>>
>> On 30 September 2013 at 16:46, Sul, Young L wrote:
>> | I know that from
>> http://cran.fhcrc.org/web/checks/check_results_RcppEigen.html
>> | RcppEigen isn?t compilable until Solaris 10 yet. I?m wondering if
>> anyone has a
>> | sense if this will be fixed in the future?
>>
>> It's tricky:
>>
>>   -- packages Rcpp${FOO} fail for various values of ${FOO} also because
>> Rcpp
>>  itself may not build (and in this particular case solaris-x86 seems
>> to
>>  pass for Rcpp (Yay!) but solaris-sparc still fails (booh) so here it
>> is
>>  indeed an RcppEigen issue)
>>
>>   -- none of the core member of the Rcpp group has access to either
>> Solaris
>>  variant (though as I recall Martyn Plummer has a solaris x86 and
>>  provided Romain access to this in the past)
>>
>>   -- so we cannot fix Rcpp, nor can we fix packages building on top of
>> Rcpp
>>
>>   -- and as "we" don't use these machine, fixes are not of high priority
>> (but
>>  if someone could provide funding, Romain may be willing to tackle
>> this
>>  but you'd have to ask him)
>>
>> However, we have always accepted patches for this.
>>
>> With RcppEigen, the added difficulty is that the listed maintainer is no
>> longer all that active around R.
>>
>
> Actually the listed maintainer was writing a reply when your reply came
> through :-)
>
> I had been willing to bet that no one used RcppEigen on Solaris and its
> failing to build on Solaris was a non-issue, but Young Sui has proved me
> wrong.
>
> However, the problems Dirk described hold for me too.  I don't have access
> to a Solaris system, either x86 or sparse, and don't know of a way that I
> could access such a system.  Even if I could get access I doubt that I
> would be able to fix the parts complained about in the logs because they
> are in the Eigen code   If I recall correctly the Eigen authors said they
> were not going to try to support Solaris because - wait for it - they don't
> have access to a Solaris system.
>
> I recall a suggestion of using openSolaris under virtualbox or some other
> virtualization software but, according to the Wikipedia page, openSolaris
> no longer exists.
>
> So that is an open spot -- if someone ambitious enough to tackle this would
>> step forward and adopt RcppEigen, we'd all appreciate it.
>>
>
> Yes, I would be willing to pass RcppEigen onto another maintainer.
> However It happens that there will be a rather complicated shuffle
> involving the Matrix, RcppEigen and lme4 packages in the near future and
> probably I should be involved in that.
>
> The situation is that a recent release (and a minor version number at
> that) of the CHOLMOD software, on which the sparse Cholesky factorization
> in Matrix and in RcppEigen is based, added a field to the
> cholmod_factor_struct.  If we use the new CHOLMOD version in Matrix, which
> we would like to do, then the installed versions of other packages will
> start to fail when the new Matrix is installed.  It is going to be messy
> because there is no way to force an upgrade of a downstream package when a
> new Matrix is installed.
>
>
___
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

Re: [Rcpp-devel] RcppEigen not compiling under Solaris 10 INTEL

[Douglas Bates]

On Mon, Sep 30, 2013 at 3:56 PM, Dirk Eddelbuettel  wrote:

>
> On 30 September 2013 at 16:46, Sul, Young L wrote:
> | I know that from
> http://cran.fhcrc.org/web/checks/check_results_RcppEigen.html
> | RcppEigen isn?t compilable until Solaris 10 yet. I?m wondering if anyone
> has a
> | sense if this will be fixed in the future?
>
> It's tricky:
>
>   -- packages Rcpp${FOO} fail for various values of ${FOO} also because
> Rcpp
>  itself may not build (and in this particular case solaris-x86 seems to
>  pass for Rcpp (Yay!) but solaris-sparc still fails (booh) so here it
> is
>  indeed an RcppEigen issue)
>
>   -- none of the core member of the Rcpp group has access to either Solaris
>  variant (though as I recall Martyn Plummer has a solaris x86 and
>  provided Romain access to this in the past)
>
>   -- so we cannot fix Rcpp, nor can we fix packages building on top of Rcpp
>
>   -- and as "we" don't use these machine, fixes are not of high priority
> (but
>  if someone could provide funding, Romain may be willing to tackle this
>  but you'd have to ask him)
>
> However, we have always accepted patches for this.
>
> With RcppEigen, the added difficulty is that the listed maintainer is no
> longer all that active around R.
>

Actually the listed maintainer was writing a reply when your reply came
through :-)

I had been willing to bet that no one used RcppEigen on Solaris and its
failing to build on Solaris was a non-issue, but Young Sui has proved me
wrong.

However, the problems Dirk described hold for me too.  I don't have access
to a Solaris system, either x86 or sparse, and don't know of a way that I
could access such a system.  Even if I could get access I doubt that I
would be able to fix the parts complained about in the logs because they
are in the Eigen code   If I recall correctly the Eigen authors said they
were not going to try to support Solaris because - wait for it - they don't
have access to a Solaris system.

I recall a suggestion of using openSolaris under virtualbox or some other
virtualization software but, according to the Wikipedia page, openSolaris
no longer exists.

So that is an open spot -- if someone ambitious enough to tackle this would
> step forward and adopt RcppEigen, we'd all appreciate it.
>

Yes, I would be willing to pass RcppEigen onto another maintainer.
However It happens that there will be a rather complicated shuffle
involving the Matrix, RcppEigen and lme4 packages in the near future and
probably I should be involved in that.

The situation is that a recent release (and a minor version number at that)
of the CHOLMOD software, on which the sparse Cholesky factorization in
Matrix and in RcppEigen is based, added a field to the
cholmod_factor_struct.  If we use the new CHOLMOD version in Matrix, which
we would like to do, then the installed versions of other packages will
start to fail when the new Matrix is installed.  It is going to be messy
because there is no way to force an upgrade of a downstream package when a
new Matrix is installed.
___
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

Re: [Rcpp-devel] Symbol not found when including an external library

[Douglas Bates]

Do you know of the nloptr and nloptwrap R packages based on NLopt.  The
nloptr package installs the NLopt library.  Unfortunately it does not put
the header files in inst/include in the source package so they will be
available in the include subdirectory of the installed package.



On Wed, Sep 4, 2013 at 11:51 AM, Simon Zehnder  wrote:

> Dear Rcpp::Users, Rcpp::Devels,
>
> I am following the ongoing discussion 'How to use external C++ libraries
> in R packages withRcpp' (Makevars are still one of my weaknesses) and I
> have a related question.
>
> I am trying to use the C++ optimization library nlopt (
> http://ab-initio.mit.edu/wiki/index.php/NLopt) in my package that is
> based on RcppArmadillo. In a first step I installed nlopt in a subfolder of
> the 'src' folder of my package structure (so I have '/src/nlopt-2.3').
> Second, I included the header file into my source file right below
> '#include '. Then I set my Makevars file to:
>
> ## Use the R_HOME indirection to support installations of multiple R
> version
> NLOPT_VERSION = 2.3
>
> NLOPT_LIBS = -lm nlopt-${NLOPT_VERSION}/lib/libnlopt_cxx.a
>
> NLOPT_INCL = -I./nlopt-${NLOPT_VERSION}/include
>
> PKG_CPPFLAGS = ${NLOPT_INCL}
> PKG_CFLAGS = -pipe ${NLOPT_INCL}
> PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` ${NLOPT_LIBS}
>
> Calling 'R CMD INSTALL packageName' compiles without errors but produces
> the well-known error during loading the shared library:
>
> Error in dyn.load(file, DLLpath = DLLpath, ...) :
>   unable to load shared object
> '/Library/Frameworks/R.framework/Versions/3.0/Resources/library/finmix/libs/finmix.so':
>
> dlopen(/Library/Frameworks/R.framework/Versions/3.0/Resources/library/finmix/libs/finmix.so,
> 6): Symbol not found: __ZNSt8__detail15_List_node_base11_M_transferEPS0_S1_
>
> I find the missing symbol (using 'nm -A libnlopt_cxx.a | grep "transfer"'
> on the shell) in the static library of nlopt (libnlopt_cxx.a) to which I
> link in the Makevars, but I do not know yet, how I have to modify the
> Makevars to include also the static library of nlopt into the shared
> library of my package. Does anyone have a suggestion how to do it?
>
>
> Best
>
> Simon
>
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] returning array from C

[Douglas Bates]

On Tue, Jun 11, 2013 at 10:51 AM, Steve Jaffe  wrote:

> Is there a way to 'wrap', say, an array of double allocated on the heap in
> C/C++ and return it to R without copying, ie as the data inside a REALSXP?
>

I don't know of any way of doing that.  R is very possessive about memory
and expects to be able to garbage-collect any memory allocated in an
SEXPREC.


> (Let's say for example that this array allocation is in a 3rd-party
> library which can't be modified)
>
> Looking at the 'wrap' code it appears that copying is done in all cases
> (except for the 'external pointer').
>

Yes.  The only way to avoid copying is if the external library allows
objects to be constructed from pre-allocated memory.  This is used in the
RcppArmadillo and RcppEigen packages.


>
> Clearly the memory ownership issue is one that would need to be handled.
> Handing over ownership to R would seem the simplest approach. But it would
> be more interesting to be able to integrate with something like
> boost::shared_array so that the lifetime would be correctly managed across
> both the C and R environments.
>
> I've read through Rcpp documentation and "Writing R Extensions" without
> coming up with a good answer, but I'm new to this and I suspect someone
> familiar with R internals (unlike me) would immediately see what was
> involved.
>
> Thanks for your help.
>
> Steve
>
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] Segfault error during simulation in Rcpp

[Douglas Bates]

The segfaults seem to be related to PutRNGstate and I don't see that you
have declared an instance of the class that causes the RNGstate to be
accessed and restored (I have forgotten the name of the class but it should
be fairly easy to find in the examples).  When you use random number
generators from within Rcpp you either need to call getRNGstate/putRNGstate
yourself or to simply declare an instance of this class.



On Mon, May 6, 2013 at 6:15 AM, Matteo Fasiolo wrote:

> Hi All,
>
>  I am trying to simulate from the following model:
>
> Y_t = Pois( phi * N_t );
> N_t = r * N_{t-1} * exp( -N_{t-1}^theta + e_t )
> e_t ~ Norm(0, sigma)
>
> so I have written a Rcpp function with prototype:
>
> genRickerCpp(int days, int nSimul, int nBurn, NumericMatrix params)
>
> where:
>
> days = length of each simulated path
> nSimul = number of simulated paths
> nBurn = number of simulations I'm going to discard before storing
> each path
> params = matrix of input parameters that can be either 1 by 4 (in which
>   case all the paths are simulated using the same parameters)
>   or nSimul by 4 (in which case each path uses a different
>   vector of parameters)
>
> This is my code:
>
> #include 
>
> using namespace Rcpp;
>
> // [[Rcpp::export]]
> NumericMatrix genRickerCpp(int days, int nSimul, int nBurn, NumericMatrix
> params)
> {
>   int nParams = params.ncol();
>   int totDays = nBurn + days;
>   bool multiParams = false;
>
>   if(nParams != 4) stop("Wrong number of parameters");
>   if(params.nrow() > 1) { multiParams = true; }
>   if(multiParams == true && params.nrow() != nSimul)
>   stop("Number of parameters vectors is different from the number of
> simulations");
>
>   double r = exp(params(0, 0));
>   double theta = exp(params(0, 1));
>   double sigma = exp(params(0, 2));
>   double phi = exp(params(0, 3));
>
>   NumericVector procNoise( rnorm( totDays * nSimul ) );
>   NumericVector initState( runif( nSimul ) );
>   NumericMatrix output( nSimul, days );
>
>   NumericVector::iterator noiseIter = procNoise.begin();
>   NumericVector::iterator initIter = initState.begin();
>
>   double currState;
>
>   for(int iRow = 0; iRow < nSimul; iRow++, initIter++)
>   {
>
>if( multiParams == true )
>{
> r = exp(params(iRow, 0));
> theta = exp(params(iRow, 1));
> sigma = exp(params(iRow, 2));
> phi = exp(params(iRow, 3));
>}
>
>currState = *initIter;
>
>for(int iCol = 1; iCol <= nBurn; iCol++, noiseIter++){
>  currState = r * currState * exp( - pow( currState, theta ) +
> *noiseIter * sigma );
>}
>
>output(iRow, 0) = rpois(1, phi * currState)[0];
>
>for(int iCol = 1; iCol < days; iCol++, noiseIter++){
>  currState = r * currState * exp( - pow( currState, theta ) +
> *noiseIter * sigma );
>  output(iRow, iCol) = rpois(1, phi * currState)[0];
>}
>
>   }
>
>   return output;
>
> }
>
>
> the function seems to work well, I tried to compare the output the an
> equivalent R function
> and I get the same results. The problem is that if I run it a lot of times:
>
> library(Rcpp)
> sourceCpp("~/Desktop/genRickerCpp.cpp")
>
> for(ii in 1:10^6){
> data <- genRickerCpp(days = 1, nSimul = 1, nBurn = 1,
>  params = matrix(log(c(r = exp(3.8), theta = 1, sigma
> = 0.3, phi = 10)), 1, 4))
> data <- as.numeric(data)
> }
>
> occasionally R crashes with error:
>
>  *** caught segfault ***
> address 0x28, cause 'memory not mapped'
>
> the strange thing is that in most cases I can call it 10^6 times without
> any error.
> I tried to go through the code in gdb, but I didn't see anything wrong. I
> also ran the previous
> R code in valgrind, and there I get the following errors while the code is
> running:
>
> ==3031== Invalid read of size 1
> ==3031==at 0x4EF7BF1: ??? (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4EF7CE5: ??? (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4EF8894: Rf_duplicate (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4EA79E7: ??? (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F25B08: Rf_eval (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F2779F: ??? (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F2791F: ??? (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F2958C: ??? (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F259AE: Rf_eval (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F61FA2: Rf_ReplIteration (in /usr/lib/R/lib/libR.so)
> ==3031==  Address 0xebd20c8 is 1,688 bytes inside a block of size 1,968
> free'd
> ==3031==at 0x4C2A82E: free (in
> /usr/lib/valgrind/vgpreload_memcheck-amd64-linux.so)
> ==3031==by 0x4F644AC: ??? (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F67965: ??? (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4F6933E: Rf_allocVector (in /usr/lib/R/lib/libR.so)
> ==3031==by 0x4E7468C: PutRNGstate (in /usr/lib/R/lib/libR.so)
> 

Re: [Rcpp-devel] Packaging - and even more packaging

[Douglas Bates]

On Sun, Feb 3, 2013 at 2:41 PM, Simon Zehnder  wrote:

> Dear Rcpp-Devels,
>
> maybe you have an answer for me as packaging is a little new for me:
>
> I took the NLopt library (http://ab-initio.mit.edu/wiki/index.php/NLopt)
> and installed it into a subfolder /nlopt-2.3 in my package-/src folder. The
> library depends on a header nlopt.hpp and a static library libnlopt.a.
>
> My Makervars in the /src-folder looks like this (using RcppArmadillo):
>
> PKG_CPPFLAGS = -Inlopt-2.3/include -g -O1
> PKG_LIBS = `$(R_HOME)/bin/Rscript -e "Rcpp:::LdFlags()"` $(LAPACK_LIBS)
> $(BLAS_LIBS) $(FLIBS) -Lnlopt-2.3/lib/libnlopt.a
>
> When I call R CMD INSTALL /MCpkgArmadillo I get
>
> Error in dyn.load(file, DLLpath = DLLpath, ...) :
>   unable to load shared object
> '/rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkgArmadillo/libs/MCpkgArmadillo.so':
>
> /rwthfs/rz/cluster/home/simo_109/opt/R/lib64/R/library/MCpkgArmadillo/libs/MCpkgArmadillo.so:
> undefined symbol: nlopt_set_munge
>
> So, it seems, that the linking to the static library has to be
> extended/modified in some way. How should I change it? I looked at
> Subsection 1.2.1.3 Compiling in subdirectories on
> http://cran.r-project.org/doc/manuals/R-exts.html but this shows an
> example, where there is code in subdirectories, that has to be archived
> into a static library.
>
> Is it even possible to put already a static library into the package and
> then compile? Or do I have to put in all source files of the NLopt library
> and create during compiling the package the static library of the NLopt
> library?
>

It is possible and we do so in the Matrix package where we have the .a
files from the various SuiteSparse C libraries but we modify the makefiles
in those libraries so that the .a file ends up in the src directory, not a
subdirectory of src.

It seems to work.


> How is it intended to be done?
>
>
> Best
>
> Simon
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] problem compiling with const in RcppEigen

[Douglas Bates]

Sorry but I am still struggling (even after reformatting the source code)
to understand what the problem is.  For one thing it is not clear to me
what the Cm_class is for.

The fundamental rule for me is that I always declare functions like

//[[Rcpp::Export]]
Eigen::VectorXd testfunc(const MapVecd ytrait, const MapMatd xmat) {
   ...
}

If you need to get a pointer to the contents of the matrix you could use a
const_cast but in that case I don't understand why you would use the Eigen
classes.  If all you want is the dimensions and a pointer to the contents
you should use NumericMatrix class.



On Thu, Jan 31, 2013 at 1:27 PM, M A  wrote:

> Right. The problem comes up because of the constructor with one of my
> classes. Here is a slightly extended minimal example:
>
> // [[Rcpp::depends(RcppEigen)]]
>
> #include 
> #include 
>
> typedef Eigen::Map MapVecd;
> typedef Eigen::Map MapMatd; // (A). This one works
> with (1) below
> //typedef Eigen::Map MapMatd; // (B). This one
> fails with (1) below
>
> class Cm_class {
> public:
> const MapMatd c_matrix;
>
> //Cm_class( double *p_pdmat, const int nr, const int nc ) : //
> Works with (A) or (B)
> Cm_class( const double *p_pdmat, const int nr, const int nc ) : //
> Works with (B) fails with (A)
> c_matrix(p_pdmat, nr, nc) {}
> };
>
> //Eigen::VectorXd testfunc(const MapVecd ytrait, MapMatd xmat) { //
> (1). This only works with (A)
> // [[Rcpp::export]]
> Eigen::VectorXd testfunc(const MapVecd ytrait) {
>
> return ytrait;
> }
>
> I really want testfunc to be as in (1) where xmat gets mapped into a
> double matrix, though I take it out above for simplicity. I have my
> class constructor taking a const double*, but that fails with the
> typedef (A), the way you recommend. To get that constructor
> initialization to work I have to use typedef (B), but that causes
> failure with (1), which is what I was originally writing about. I can
> certainly work around this by removing the const from the double* in
> the class constructor, etc, but I guess there's the larger question of
> whether the RcppEigen 'as' function should be able to form an
> Eigen::Map. Thanks for the info regarding the
> point of failure.
>
> Mark A
>
> On Thu, Jan 31, 2013 at 11:34 AM, Douglas Bates 
> wrote:
> > I should check the syntax before sending the message.
> >
> > On Thu, Jan 31, 2013 at 11:30 AM, Douglas Bates 
> wrote:
> >>
> >> The as function in RcppEigen can form an
> >>
> >> Eigen::Map
> >>
> >> but doesn't know how to form an
> >>
> >> Eigen::Map
> >>
> >> Generally I find that I confuse myself less when the const is outside
> the
> >> Eigen::Map.  In other words, I use something like
> >>
> >> typedef Eigen::Map mMat;
> >>
> >> const mMat M(mat_from_R);
> >
> >
> > That should be
> >
> > const mMat M = as(mat_from_R);
> >
> > Well, at least I think so.
> >
> > The bottom line is that I find it easier to use const in the declaration
> of
> > an instance of a class, not in the typedef.
> >>
> >>
> >> This tells me that M is a mapped MatrixXd object which is constructed
> from
> >> the SEXP mat_from_R and is read-only.
> >>
> >>
> >>
> >> On Thu, Jan 31, 2013 at 11:01 AM, M A  wrote:
> >>>
> >>> I am compiling some source code using sourceCpp and getting a
> >>> compiling failure when the source has the line with the const keyword
> >>> but no failure when it does not, in the below source.
> >>>
> >>> // file: minimaltest.cpp
> >>> // [[Rcpp::depends(RcppEigen)]]
> >>>
> >>> #include 
> >>> #include 
> >>>
> >>> typedef Eigen::Map MapVecd;
> >>> //typedef Eigen::Map MapMatd; // This one works
> >>> typedef Eigen::Map MapMatd; // This one fails
> >>>
> >>> // [[Rcpp::export]]
> >>> Eigen::VectorXd testfunc(const MapVecd ytrait, MapMatd xmat) {
> >>>
> >>> return ytrait;
> >>> }
> >>>
> >>> This is the error:
> >>>
> >>> > sourceCpp("minimaltest.cpp")
> >>> g++ -arch x86_64 -I/Library/Frameworks/R.framework/Resources/include
> >>> -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG
> >>> -I/usr/local/include
> >>>
> >>>
> -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include&q

Re: [Rcpp-devel] problem compiling with const in RcppEigen

[Douglas Bates]

I should check the syntax before sending the message.

On Thu, Jan 31, 2013 at 11:30 AM, Douglas Bates  wrote:

> The as function in RcppEigen can form an
>
> Eigen::Map
>
> but doesn't know how to form an
>
> Eigen::Map
>
> Generally I find that I confuse myself less when the const is outside the
> Eigen::Map.  In other words, I use something like
>
> typedef Eigen::Map mMat;
>
> const mMat M(mat_from_R);
>

That should be

const mMat M = as(mat_from_R);

Well, at least I think so.

The bottom line is that I find it easier to use const in the declaration of
an instance of a class, not in the typedef.

>
> This tells me that M is a mapped MatrixXd object which is constructed from
> the SEXP mat_from_R and is read-only.
>
>
>
> On Thu, Jan 31, 2013 at 11:01 AM, M A  wrote:
>
>> I am compiling some source code using sourceCpp and getting a
>> compiling failure when the source has the line with the const keyword
>> but no failure when it does not, in the below source.
>>
>> // file: minimaltest.cpp
>> // [[Rcpp::depends(RcppEigen)]]
>>
>> #include 
>> #include 
>>
>> typedef Eigen::Map MapVecd;
>> //typedef Eigen::Map MapMatd; // This one works
>> typedef Eigen::Map MapMatd; // This one fails
>>
>> // [[Rcpp::export]]
>> Eigen::VectorXd testfunc(const MapVecd ytrait, MapMatd xmat) {
>>
>> return ytrait;
>> }
>>
>> This is the error:
>>
>> > sourceCpp("minimaltest.cpp")
>> g++ -arch x86_64 -I/Library/Frameworks/R.framework/Resources/include
>> -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG
>> -I/usr/local/include
>>
>> -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include"
>>
>> -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include"
>>   -fPIC  -g -O2  -c minimaltest.cpp -o minimaltest.o
>> Error in sourceCpp("minimaltest.cpp") :
>>   Error 1 occurred building shared library.
>>
>> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include/Rcpp/internal/Exporter.h:
>> In constructor ‘Rcpp::traits::Exporter::Exporter(SEXPREC*) [with T
>> = Eigen::Map> -0x1, 0, -0x1, -0x1>,
>> 0, Eigen::Stride<0, 0> >]’:
>>
>> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include/Rcpp/as.h:67:
>>   instantiated from ‘T Rcpp::internal::as(SEXPREC*,
>> Rcpp::traits::r_type_generic_tag) [with T = Eigen::Map> Eigen::Matrix> -0x1, -0x1>, 0, Eigen::Stride<0, 0>
>> >]’
>>
>> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include/Rcpp/as.h:112:
>>   instantiated from ‘T Rcpp::as(SEXPREC*) [with T = Eigen::Map> Eigen::Matrix> -0x1, -0x1>, 0, Eigen::Stride<0, 0>
>> >]’
>> minimaltest.cpp:35:   instantiated from here
>>
>> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include/Rcpp/internal/Exporter.h:31:
>> error: no matching function for call to ‘Eigen::Map> Eigen::Matrix> -0x1, -0x1>, 0, Eigen::Stride<0, 0>
>> >::Map(SEXPREC*&)’
>>
>> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include/Eigen/src/Core/Map.h:164:
>> note: candidates are: Eigen::Map> StrideType>::Map(typename Eigen::MapBase> MapOptions, StrideType>,
>> (Eigen::internal::accessors_level> MapOptions, StrideType> >::has_write_access ?  WriteAccessors :
>> ReadOnlyAccessors)>::PointerType, typename
>> Eigen::internal::traits> StrideType> >::Index, typename
>> Eigen::internal::traits> StrideType> >::Index, const StrideType&) [with PlainObjectType = const
>> Eigen::Matrix> -0x1, -0x1>, int MapOptions = 0,
>> StrideType = Eigen::Stride<0, 0>]
>>
>> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include/Eigen/src/Core/Map.h:151:
>> note: Eigen::Map> StrideType>::Map(typename Eigen::MapBase> MapOptions, StrideType>,
>> (Eigen::internal::accessors_level> MapOptions, StrideType> >::has_write_access ?  WriteAccessors :
>> ReadOnlyAccessors)>::PointerType, typename
>> Eigen::internal::traits> StrideType> >::Index, const StrideType&) [with PlainObjectType = const
>> Eigen::Matrix> -0x1, -0x1>, int MapOptions = 0,
>> StrideType = Eigen::Stride<

Re: [Rcpp-devel] problem compiling with const in RcppEigen

[Douglas Bates]

The as function in RcppEigen can form an

Eigen::Map

but doesn't know how to form an

Eigen::Map

Generally I find that I confuse myself less when the const is outside the
Eigen::Map.  In other words, I use something like

typedef Eigen::Map mMat;

const mMat M(mat_from_R);

This tells me that M is a mapped MatrixXd object which is constructed from
the SEXP mat_from_R and is read-only.



On Thu, Jan 31, 2013 at 11:01 AM, M A  wrote:

> I am compiling some source code using sourceCpp and getting a
> compiling failure when the source has the line with the const keyword
> but no failure when it does not, in the below source.
>
> // file: minimaltest.cpp
> // [[Rcpp::depends(RcppEigen)]]
>
> #include 
> #include 
>
> typedef Eigen::Map MapVecd;
> //typedef Eigen::Map MapMatd; // This one works
> typedef Eigen::Map MapMatd; // This one fails
>
> // [[Rcpp::export]]
> Eigen::VectorXd testfunc(const MapVecd ytrait, MapMatd xmat) {
>
> return ytrait;
> }
>
> This is the error:
>
> > sourceCpp("minimaltest.cpp")
> g++ -arch x86_64 -I/Library/Frameworks/R.framework/Resources/include
> -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG
> -I/usr/local/include
>
> -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include"
>
> -I"/Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include"
>   -fPIC  -g -O2  -c minimaltest.cpp -o minimaltest.o
> Error in sourceCpp("minimaltest.cpp") :
>   Error 1 occurred building shared library.
>
> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include/Rcpp/internal/Exporter.h:
> In constructor ‘Rcpp::traits::Exporter::Exporter(SEXPREC*) [with T
> = Eigen::Map -0x1, 0, -0x1, -0x1>,
> 0, Eigen::Stride<0, 0> >]’:
>
> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include/Rcpp/as.h:67:
>   instantiated from ‘T Rcpp::internal::as(SEXPREC*,
> Rcpp::traits::r_type_generic_tag) [with T = Eigen::Map Eigen::Matrix -0x1, -0x1>, 0, Eigen::Stride<0, 0>
> >]’
>
> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include/Rcpp/as.h:112:
>   instantiated from ‘T Rcpp::as(SEXPREC*) [with T = Eigen::Map Eigen::Matrix -0x1, -0x1>, 0, Eigen::Stride<0, 0>
> >]’
> minimaltest.cpp:35:   instantiated from here
>
> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/Rcpp/include/Rcpp/internal/Exporter.h:31:
> error: no matching function for call to ‘Eigen::Map Eigen::Matrix -0x1, -0x1>, 0, Eigen::Stride<0, 0>
> >::Map(SEXPREC*&)’
>
> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include/Eigen/src/Core/Map.h:164:
> note: candidates are: Eigen::Map StrideType>::Map(typename Eigen::MapBase MapOptions, StrideType>,
> (Eigen::internal::accessors_level MapOptions, StrideType> >::has_write_access ?  WriteAccessors :
> ReadOnlyAccessors)>::PointerType, typename
> Eigen::internal::traits StrideType> >::Index, typename
> Eigen::internal::traits StrideType> >::Index, const StrideType&) [with PlainObjectType = const
> Eigen::Matrix -0x1, -0x1>, int MapOptions = 0,
> StrideType = Eigen::Stride<0, 0>]
>
> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include/Eigen/src/Core/Map.h:151:
> note: Eigen::Map StrideType>::Map(typename Eigen::MapBase MapOptions, StrideType>,
> (Eigen::internal::accessors_level MapOptions, StrideType> >::has_write_access ?  WriteAccessors :
> ReadOnlyAccessors)>::PointerType, typename
> Eigen::internal::traits StrideType> >::Index, const StrideType&) [with PlainObjectType = const
> Eigen::Matrix -0x1, -0x1>, int MapOptions = 0,
> StrideType = Eigen::Stride<0, 0>]
>
> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include/Eigen/src/Core/Map.h:139:
> note: Eigen::Map StrideType>::Map(typename Eigen::MapBase MapOptions, StrideType>,
> (Eigen::internal::accessors_level MapOptions, StrideType> >::has_write_access ?  WriteAccessors :
> ReadOnlyAccessors)>::PointerType, const StrideType&) [with
> PlainObjectType = const Eigen::Matrix -0x1, 0, -0x1, -0x1>,
> int MapOptions = 0, StrideType = Eigen::Stride<0, 0>]
>
> /Library/Frameworks/R.framework/Versions/2.15/Resources/library/RcppEigen/include/Eigen/src/Core/Map.h:106:
> note: Eigen::Map -0x1, -0x1, 0, -0x1,
> -0x1>, 0, Eigen::Stride<0, 0> >::Map(const
> Eigen::Map -0x1, 0, -0x1, -0x1>,
> 0, Eigen::Stride<0, 0> >&)
> make: *** [minimaltest.o] Error 1
>
> Here is my session info:
> > sessionInfo()
> R version 2.15.2 (2012-10-26)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_U

Re: [Rcpp-devel] Rcpp with C++ and Fortran

[Douglas Bates]

Could you clarify a bit, Rodney.  When you say you want to "link it with
Fortran code" are you using the word "link" in the sense of linking to
object code or do you mean that you will copy the Fortran source files into
your package's src directory and have them compiled with other source
files?  The second, which I would describe as "calling Fortran subroutines
from C++ code", is definitely easier.

To actually call a Fortran subroutine you need to know about appending
underscores on entry point names and the fact that all arguments to Fortran
functions are passed by reference.  That is, everything is a pointer.  It
might be easiest if you could post the calling sequence of a particular
subroutine so we could show the exact code.  I'm happy to carry on the
conversation off-list if this seems like too much noise for this list.
 Discussing calling Fortran code on a C++ list is a "keep 'em down on the
farm after they've seen Paris" situation.


On Wed, Jan 23, 2013 at 2:16 PM, Rodney Sparapani  wrote:

> Hi!
>
> Quick question:  has anyone used Rcpp with C++ and Fortran?
>
> Longer version:  We have developed some C++ code to do MCMC sampling
> with RcppEigen.  It works pretty well:  about a 70X speed
> improvement over R code.  We have been able to link this with some
> ARMS C code rather easily.  However, now we want to link it with
> Fortran code from the DPpackage R package.
>
> I have looked through as many Rcpp docs as I could get my hands, but I
> didn't see anything about Fortran (plus googling).  Anyways, what is
> the best way to call Fortran functions from our C++ code that will be
> used via Rcpp/inline?  I've played around with f2c to convert the
> Fortran code to C (but it makes me nervous ;o)  Should I be working
> from a package skeleton?  Basically, any advice is welcome.  Thanks!
>
> --
> Rodney Sparapani, PhD  Center for Patient Care and Outcomes Research
> Sr. Biostatistician   http://www.mcw.edu/pcor
> 4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
> WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA
> __**_
> 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
>
___
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

Re: [Rcpp-devel] RcppEigen: Getting and setting dimnames of sparse matrix

[Douglas Bates]

On Sat, Jan 19, 2013 at 10:38 AM, Douglas Bates  wrote:

>
On thinking about this a bit more, the dimnames should be cloned before
being assigned, otherwise you end up with two references to the same
storage.  And once you try to clone you get into the SlotProxy area where
you need to know what you are cloning.  One way of doing this is shown in
the enclosed.  An alternative is

Xout.slot("Dimnames") = clone(as((Xin.slot("Dimnames";

R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Matrix)
Loading required package: lattice
> library(RcppEigen)
Loading required package: Rcpp
> library(inline)
> 
> ugM <- new("dgCMatrix"
+, i = c(1L, 0L, 2L, 1L)
+, p = c(0L, 1L, 3L, 4L)
+, Dim = c(3L, 3L)
+, Dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
+, x = c(1, 1, 1, 1)
+, factors = list())
> ugM
3 x 3 sparse Matrix of class "dgCMatrix"
  a b c
a . 1 .
b 1 . 1
c . 1 .
> str(ugM)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i   : int [1:4] 1 0 2 1
  ..@ p   : int [1:4] 0 1 3 4
  ..@ Dim : int [1:2] 3 3
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:3] "a" "b" "c"
  .. ..$ : chr [1:3] "a" "b" "c"
  ..@ x   : num [1:4] 1 1 1 1
  ..@ factors : list()
> 
> src <- '
+ typedef Eigen::SparseMatrix SpMat;
+ S4Xin(XX_);
+ SpMat X(as(XX_));
+ S4Xout(wrap(X));
+ Xout.slot("Dimnames") = clone(List(Xin.slot("Dimnames")));
+ return(Xout);
+ '
> 
> names_ <- cxxfunction(signature(XX_="dsCMatrix"), plugin="RcppEigen",
+ body=src)
> 
> names_(ugM)
3 x 3 sparse Matrix of class "dgCMatrix"
  a b c
a . 1 .
b 1 . 1
c . 1 .
> 
> 
> proc.time()
   user  system elapsed 
 10.392   0.496  10.918 
___
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

Re: [Rcpp-devel] RcppEigen: Getting and setting dimnames of sparse matrix

[Douglas Bates]

The Eigen sparse matrix only contains the numeric information, it doesn't
contain information like Dimnames that are specific to the R structure.  If
you want to access the slots or set their values it is easiest to do that
with the Rcpp::S4 class object.

Of course it helps to use R's str function first to ensure that you know
what the structure actually is.

R version 2.15.2 (2012-10-26) -- "Trick or Treat"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Matrix)
Loading required package: lattice
> library(RcppEigen)
Loading required package: Rcpp
> library(inline)
> 
> ugM <- new("dgCMatrix"
+, i = c(1L, 0L, 2L, 1L)
+, p = c(0L, 1L, 3L, 4L)
+, Dim = c(3L, 3L)
+, Dimnames = list(c("a", "b", "c"), c("a", "b", "c"))
+, x = c(1, 1, 1, 1)
+, factors = list())
> ugM
3 x 3 sparse Matrix of class "dgCMatrix"
  a b c
a . 1 .
b 1 . 1
c . 1 .
> str(ugM)
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i   : int [1:4] 1 0 2 1
  ..@ p   : int [1:4] 0 1 3 4
  ..@ Dim : int [1:2] 3 3
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:3] "a" "b" "c"
  .. ..$ : chr [1:3] "a" "b" "c"
  ..@ x   : num [1:4] 1 1 1 1
  ..@ factors : list()
> 
> src <- '
+ typedef Eigen::SparseMatrix SpMat;
+ S4Xin(XX_);
+ SpMat X(as(XX_));
+ S4Xout(wrap(X));
+ Xout.slot("Dimnames") = Xin.slot("Dimnames");
+ return(Xout);
+ '
> 
> names_ <- cxxfunction(signature(XX_="dsCMatrix"), plugin="RcppEigen",
+ body=src)
> 
> names_(ugM)
3 x 3 sparse Matrix of class "dgCMatrix"
  a b c
a . 1 .
b 1 . 1
c . 1 .
> 
> 
> proc.time()
   user  system elapsed 
 10.356   0.436  10.859 
___
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

Re: [Rcpp-devel] Regarding question on a suitable IDE for C++ project using Rcpp in Window (7) OS

[Douglas Bates]

The inability to use Visual Studio for compiling code to use with R is not
peculiar to Rcpp.  The Windows version of R is compiled with mingw
compilers and requires compatible compilers for other code.  Remember that
R is an Open Source project and Visual Studio isn't.

As far as an IDE goes, RStudio (rstudio.org) provides integration of Rcpp
in a GUI environment.  Support for Rcpp in RStudio is likely to keep
improving as JJ Allaire is actively involved in the development of Rcpp.
 Many of us use ESS in emacs but that's because we predate GUI's. (I tell
you, sonny, back in my day we had terminals with 80 by 24 character screens
and 300 baud modems and we still created better code than these
whippersnappers can today! :-)



On Mon, Jan 7, 2013 at 10:28 AM, Lincoln Xu  wrote:

> Hi, Rcpp experts,
>
> While I am very impressed by the power of Rcpp described by a nice
> introduction in your website. It seems to me that it is not easy to set up
> an IDE for developing C++ using Rcpp in window. As stated in FAQ, there is
> NOT A CHANGE to use VS with Rcpp and R. I wonder if you could kindly
> point me to some online resource of workable IDE in window OS for the devel
> purpose ( I mean it can compile C++ code and debug the C++ code running
> from R.
>
> I am also very puzzled by the statement that there is NOT A CHANGE to use
> VS with Rcpp and R in FAQ. If I understand correctly, original Rcpp written
> by
> by Dominick Samperi was indeed working with VS. So what new technique
> deployed prevents Rcpp and R to works with VS. If all possible, could you
> please provide information for my education.
>
> Regards,
>
> Lincoln
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] receiving R matrices: question about RcppArmadillo example fastlm

[Douglas Bates]

The important question is whether the as copies the contents of
the matrix from the SEXP passed by R.  Creating an Rcpp::NumericMatrix does
not copy the contents and converting the NumericMatrix to an arma::mat with
that trailing 'false' also just copies the pointer to the contents, not the
contents itself.  I'm not sure what the "as" method does.

I prefer the approach in RcppEigen where you have two different types of
objects with "as" methods, the Eigen::MatrixXd object (copies the contents)
and the Eigen::Map object which only copies the pointer.
 I generally declare the objects created from SEXP's as

const Eigen::Map  X(Rcpp::as< ...>(Xs));

so the compiler can take responsibility for ensuring that I don't
unintentionally modify the contents.

There are several examples in the RcppEigen vignette that Dirk and I wrote.




On Wed, Jan 9, 2013 at 11:58 PM, Paul Johnson  wrote:

> Hello, Rcpp land.
>
> I'm still collecting idioms and examples for the usage of Rcpp. Today
> I'm looking at ways that people manage the importation of R matrices.
>
> Consider fasLm.r file in the RcppArmadillo package's examples folder.
> It uses this this three-step process to bring in a matrix from R. This
> leads up to the fLmTwoCasts function:
>
> src <- '
> Rcpp::NumericMatrix Xr(Xs);
> Rcpp::NumericVector yr(ys);
> int n = Xr.nrow(), k = Xr.ncol();
> arma::mat X(Xr.begin(), n, k, false);
> arma::colvec y(yr.begin(), yr.size(), false);
>...
>
> It seems to waste some effort. It allocates Xr, then an iterator, and
> then writes it into X. 3 objects to get one?
>
> If I were writing this, I would take the more direct Rcpp::as route,
> like the varSimulation.r example does
> in the same folder:
>
> src <- '
> arma::mat X = Rcpp::as(Xs);
> arma::colvec y = Rcpp::as(ys);
> int n = X.n_rows; int k = X.n_cols;
> ...
>
> I tested this, it gives the exact same answers.
>
> I was certain it would be quicker. It creates fewer objects, it is
> more direct.  But, as far as I can see in a simple test, the Rcpp::as
> way is not faster.  If anything, it is slower, by a minute fraction.
>
> I really thought it would be faster, though, and lighter on memory
> usage. What do you think? I'm looking at as.h in the Rcpp code to try
> to figure it out. But, wow! It that complicated, or what?
>
>
> Here's my benchmark example, which I include only as evidence that I
> really tried before asking here.
>
> ## built from RcppArmadillo examples/fastLm.r
> library(inline)
> library(rbenchmark)
>
> src <- '
> Rcpp::NumericMatrix Xr(Xs);
> Rcpp::NumericVector yr(ys);
> int n = Xr.nrow(), k = Xr.ncol();
> arma::mat X(Xr.begin(), n, k, false);
> arma::colvec y(yr.begin(), yr.size(), false);
> int df = n - k;
>
> // fit model y ~ X, extract residuals
> arma::colvec coef = arma::solve(X, y);
> arma::colvec res  = y - X*coef;
>
> double s2 = std::inner_product(res.begin(), res.end(),
>res.begin(), 0.0)/df;
> // std.errors of coefficients
> arma::colvec sderr = arma::sqrt(s2 *
>arma::diagvec(arma::pinv(arma::trans(X)*X)));
>
> return Rcpp::List::create(Rcpp::Named("coefficients")=coef,
>   Rcpp::Named("stderr")  =sderr,
>   Rcpp::Named("df")  =df);
> '
>
> fLmTwoCasts <- cxxfunction(signature(Xs="numeric", ys="numeric"),
>src, plugin="RcppArmadillo")
>
>
> src <- '
> arma::mat X = Rcpp::as(Xs);
> arma::colvec y = Rcpp::as(ys);
> int n = X.n_rows; int k = X.n_cols;
>
> int df = n - k;
>
> // fit model y ~ X, extract residuals
> arma::colvec coef = arma::solve(X, y);
> arma::colvec res  = y - X*coef;
>
> double s2 = std::inner_product(res.begin(), res.end(),
>res.begin(), 0.0)/df;
> // std.errors of coefficients
> arma::colvec sderr = arma::sqrt(s2 *
>arma::diagvec(arma::pinv(arma::trans(X)*X)));
>
> return Rcpp::List::create(Rcpp::Named("coefficients")=coef,
>   Rcpp::Named("stderr")  =sderr,
>   Rcpp::Named("df")  =df);
> '
>
> pjTwoCasts <- cxxfunction(signature(Xs="numeric", ys="numeric"),
>src, plugin="RcppArmadillo")
>
>
> N <- 1
>
> X <- cbind(1, rnorm(N), rpois(N, lambda=0.8), rnorm(N), rnorm(N),
> rnorm(N), rnorm(N), rnorm(N), rnorm(N), rnorm(N), rnorm(N), rnorm(N),
> rnorm(N), rnorm(N))
> y <- rnorm(N)
>
> res <- benchmark(
> fLmTwoCasts(X, y),
> pjTwoCasts(X, y),
> columns = c("test", "replications", "relative",
> "elapsed", "user.self", "sys.self"),
> order="relative",
> replications=5000)
>
> print(res)
>
> --
> Paul E. Johnson
> Professor, Political Science  Assoc. Director
> 1541 Lilac Lane, Room 504  Center for Research Methods
> University of Kansas University of Kansas
> http://pj.freef

Re: [Rcpp-devel] Wrapping a C struct in C++ for constructor/destructor

[Douglas Bates]

On Fri, Dec 14, 2012 at 12:07 PM, Richard Downe wrote:

> Since cholmod has specific functions for all the allocation/free tasks
> involved in the different data structs, if you're really ambitious, you
> could probably create a templated mgr c++ class, where the template args
> are the specific memory management functions for the wrapped struct.
>

Perhaps I can continue this discussion with you off-list.  I'm not sure
what you mean by an mgr C++ class but I am willing to learn.

One of the complications with the CHOLMOD memory-management model is the
assumption that all memory allocation and freeing is done by CHOLMOD
functions.  So, for example, the cholmod_free_sparse function frees all the
arrays in the struct then frees the struct.  We want to have the various
sparse matrices and Cholesky decompositions accessible from R.  One could
accomplish this by always using cholmod_alloc_sparse, etc. and copying the
R object into the Cholmod-allocated storage and vice-versa but, considering
that many of these matrices are read-only and for complex models applied to
large data sets the single largest object is the Cholesky factor, that is
wasteful both in memory usage and in time.

Also, essentially every cholmod_* function takes a pointer to a
cholmod_common struct which is assumed to be global or at least consistent
throughout the series of calls in a problem solution.  Maintaining a global
C struct like that through several .Call instances involves some dancing
around, especially if you want to modify members of the struct at the R
level.  We had to do some rather messy manipulation in the Matrix package
to achieve this.

Then you could automagically generate ctor/dtor/copyctor/oper= code for
> each struct, and access the internal elements through a templated _ptr*
> member or something similar, and provide as()/wrap() specializations for
> the specific structs you wish to export to R.


Well at least this discussion has taught me what ctor, dtor etc. mean.  I
have seen the names but hadn't looked them up to discover that ctor is an
abbreviation of constructor and dtor is an abbreviation of destructor.


>
>
>
>  On 12/14/2012 10:43 AM, Romain Francois wrote:
>
>> I would create a class that contains a pointer to the c struct.
>>
>> This obviously leaves the problem of freeing the memory, i.e does the c++
>> object own the pointer.
>>
>> Maybe a strategy using some reference counting smart pointer.
>>
>> In RcppEigen we have not resolved this, we simply exposed a free methid
>> for the class and the programmer has to call it. Not the best strategy
>> perhaps.
>>
>>
>> You could also prevent copy ctor and assignment op by making them private
>> without an implementation, and always pass them by reference. Then you can
>> consider that your c++ object owns the pointer.
>>
>>
>>
>>
>>
>> Le 14 déc. 2012 à 17:10, Douglas Bates  a écrit :
>>
>>  This is more a C++ question than an Rcpp question although, of course,
>>> the application will be through Rcpp.
>>>
>>> I want to use C functions defined in the CHOLMOD package using Rcpp.
>>>  The arguments to these functions and the returned values are usually
>>> pointers to C structs. There are several such structs defined in CHOLMOD
>>> and it gets tedious to use code that allocates storage for the struct,
>>> populates the struct, uses it and then frees it.  (If you look in the
>>> Matrix package sources you will see several inelegant utility functions and
>>> macros to perform these operations).
>>>
>>> For an Rcpp/C++ approach I would prefer to use constructors and
>>> destructors for the memory allocation and freeing.  Is the recommended
>>> approach to create a C++ class or struct that extends the C struct or to
>>> create a C++ class with the C struct as a member of the class?
>>>
>>> I should be able to create the code for myself if I just know where to
>>> look for examples.  I tried StackOverflow but I was unable to narrow down
>>> my search.  Pointers to sections of "C++ Annotations" or to StackOverflow
>>> questions would be appreciated.
>>> __**_
>>> 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<https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel>
>>>
>> __**_
>> Rcpp-devel mailing list
>> Rcpp-devel@lists.r-forge.r-**project.org
>> https://lists.r-forge.r-**project.org/cgi-bin/mailman/**
>

[Rcpp-devel] Wrapping a C struct in C++ for constructor/destructor

[Douglas Bates]

This is more a C++ question than an Rcpp question although, of course, the
application will be through Rcpp.

I want to use C functions defined in the CHOLMOD package using Rcpp.  The
arguments to these functions and the returned values are usually pointers
to C structs. There are several such structs defined in CHOLMOD and it gets
tedious to use code that allocates storage for the struct, populates the
struct, uses it and then frees it.  (If you look in the Matrix package
sources you will see several inelegant utility functions and macros to
perform these operations).

For an Rcpp/C++ approach I would prefer to use constructors and destructors
for the memory allocation and freeing.  Is the recommended approach to
create a C++ class or struct that extends the C struct or to create a C++
class with the C struct as a member of the class?

I should be able to create the code for myself if I just know where to look
for examples.  I tried StackOverflow but I was unable to narrow down my
search.  Pointers to sections of "C++ Annotations" or to StackOverflow
questions would be appreciated.
___
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

[Rcpp-devel] Wrapping a C struct in a C++ class/struct

[Douglas Bates]

This is more of a C++ question than an Rcpp question
___
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

Re: [Rcpp-devel] For a Matrix A, is A(i, j) or A[i, j] preferred as an accessor for an element

[Douglas Bates]

On Mon, Dec 10, 2012 at 10:38 AM, Romain Francois
wrote:

> Le 10/12/12 17:29, Douglas Bates a écrit :
>
>> or does it matter?
>>
>
> A[i,j] is wrong, not valid C or C++ code. so there is only one choice.
>

Indeed.  Thanks.  Somehow I managed to convince myself that I had used
A[i,j] at some point and the compiler accepted it.  Apparently my memory
isn't what it used to be.
___
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

[Rcpp-devel] For a Matrix A, is A(i, j) or A[i, j] preferred as an accessor for an element

[Douglas Bates]

or does it matter?
___
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

Re: [Rcpp-devel] Sparse matrices with RcppArmadillo

[Douglas Bates]

On Sat, Dec 8, 2012 at 10:35 AM, c s  wrote:

> Armadillo sparse matrices are stored in Compressed Sparse Column format:
>
>
> http://en.wikipedia.org/wiki/Sparse_matrix#Compressed_sparse_column_.28CSC_or_CCS.29
>
> This layout is used by a majority of external solvers.
>
> It would be far more efficient to take this layout into account when
> copying matrices, instead of blindly (and slowly) copying element by
> element.
>

dgCMatrix objects also use compressed sparse column format with zero-based
indices so it is likely that it would only be necessary to copy the
contents of the arrays of column pointers (the "p" slot), row indices (the
"i" slot) and values of the non-zeros (the "x" slot).  It may even be
possible to just map the pointers for a read-only sparse matrix.

Romain:
 In general it is very slow to insert non-zero elements into this format.
 In the worst case the entire structure must be copied and extended for
each insertion.  We have to keep telling people about this when they tell
us that sparse matrices in R are so slow to work with.

 You have seen what the layout of a dgCMatrix is so it is only necessary to
find the desired components of a arma::sp_mat object.  I'll take a look
once I get some things installed.  If it is desirable the as and wrap
methods for the Eigen::SparseMatrix and Eigen::MappedSparseMatrix
classes should be fairly easy to adapt for arma::sp_mat class.

>
> On Sunday, December 9, 2012, Romain Francois 
> wrote:
> > Le 08/12/12 09:45, Søren Højsgaard a écrit :
> >>
> >> Dear all,
> >>
> >> I want to use a matrix (of type "dgCMatrix" from the Matrix package) in
> RcppArmadillo, so I do:
> >>
> >> library(inline)
> >> src <- '
> >> using namespace arma;
> >> using namespace Rcpp;
> >> SpMat X = as >(XX_);
> >> '
> >> foo <- cxxfunction(signature(XX_=""), body=src, plugin="RcppArmadillo")
> >>
> >> - but this fails. It seems to me (browsing the web) that SpMat are
> supported, but I might be wrong here. I have no indication that dgCMatrix
> matrices can be "converted" to SpMat's. I know that I can work with
> dgCMatrix matrices with RcppEigen, but what I need is to extract
> submatrices and that is very easy to do with RcppArmadillo.
> >>
> >> Any thoughts? Thanks in advance.
> >> Regards
> >> Søren
> >
> > Doug might know better about the internals of Matrix types. This is just
> following the recipee from how these are handled in RcppEigen:
> >
> > #include 
> > // [[Rcpp::depends(RcppArmadillo)]]
> > using namespace Rcpp ;
> >
> > // [[Rcpp::export]]
> > void convert(S4 mat){
> > IntegerVector dims = mat.slot( "Dim" ) ;
> > IntegerVector i = mat.slot( "i" ) ;
> > IntegerVector p = mat.slot( "p" ) ;
> > NumericVector x = mat.slot( "x" ) ;
> >
> > int nrow = dims[0], ncol = dims[1] ;
> > arma::sp_mat res( nrow, ncol) ;
> > for(int j = 0; j < ncol; ++j) {
> > for (int k = p[j]; k < p[j + 1]; ++k) res( i[k], j ) = x[k];
> > }
> > std::cout << res << std::endl ;
> >
> > }
> >
> >
> > /*** R
> > require(Matrix)
> > i <- c(1,3:8); j <- c(2,9,6:10); x <- 7 * (1:7)
> > ( A <- sparseMatrix(i, j, x = x) )
> >
> > convert(A)
> > ***/
> >
> > I don't think there is a better way to fill multiple values ina SpMat,
> maybe Conrad has insights.
> >
> > Romain
> >
> > --
> > Romain Francois
> > Professional R Enthusiast
> > +33(0) 6 28 91 30 30
> >
> > R Graph Gallery: http://gallery.r-enthusiasts.com
> >
> > 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
> >
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] How to increase the coding efficiency

[Douglas Bates]

I don't know as much about Armadillo as I do about Eigen so I will cheat
and write using RcppEigen instead of RcppArmadillo.  Page 6 of the Eigen
tutorial at http://eigen.tuxfamily.org/dox/ discusses decompositions and
solving linear systems of equations.  One of the simplest ways of solving a
weighted least squares system, if you only want the solution itself, is
with a QR decomposition of which there are 3 in Eigen.  The simplest is
HouseholderQR (named after A.S. Householder who first described an
elementary reflector).  For a least squares solution (assuming the model
matrix has full column rank) you simply use the solve method. For a
weighted least squares solution you premultiply both the model matrix and
the response by the square roots of the weights.

It can be collapsed to a one-liner, as in the enclosed.  A test that it
works for unweighted least squares (i.e. using unit weights) is

> sourceCpp("/tmp/wtls.cpp")
> coef(lm(optden ~ carb, Formaldehyde))(Intercept)carb
0.005085714 0.876285714
> with(Formaldehyde, wtls(model.matrix(~carb), optden, rep.int(1,
length(optden
$betahat
[1] 0.005085714 0.876285714

The use of the other decompositions (Cholesky, ColPivHouseholderQR,
Symmetric Eigendecomposition, SVD) are described in the RcppEigen vignette.

By the way, the sourceCpp, evalCpp and cppFunction capabilities developed
by the Rcpp authors is very close to magic.  They are to be congratulated
on a remarkable accomplishment.


On Wed, Dec 5, 2012 at 1:16 PM, Tim Triche, Jr. wrote:

> this part will always make your code crawl along:
>
>
> On Wed, Dec 5, 2012 at 11:09 AM, Honglang Wang  > wrote:
>
>>   arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y;
>>
>
> first time I wrote a GLM engine, I wrote it the way statistics books
> illustrate it (i.e. actually invert things to do IWLS) and it crawled.  I
> took it to a physicist friend who went through alternate stages of disgust
> and laughter then told me never to invert anything I didn't have to.
>
> You should take Doug Bates' advice, it could save you a great deal of
> time.  Never bit fiddle when you can use a better algorithm.
>
>
> --
> *A model is a lie that helps you see the truth.*
> *
> *
> Howard Skipper
>
>
> ___
> 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
>
// [[Rcpp::depends(RcppEigen)]]
#include 

typedef Eigen::MatrixXd Mat;
typedef Eigen::MapMMat;
typedef Eigen::HouseholderQRQR;
typedef Eigen::VectorXd Vec;
typedef Eigen::MapMVec;

// [[Rcpp::export]]

Rcpp::List wtls(const MMat X, const MVec y, const MVec sqrtwts) {
return Rcpp::List::create(Rcpp::Named("betahat") =
			  QR(sqrtwts.asDiagonal()*X).solve(sqrtwts.asDiagonal()*y));
}
___
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

Re: [Rcpp-devel] How to increase the coding efficiency

[Douglas Bates]

On Tue, Dec 4, 2012 at 9:39 PM, Honglang Wang wrote:

> Yes, the main issue for my coding is the allocation of memory.  And I have
> fixed one of the biggest memory allocation issue: 4000 by 4000 diagonal
> matrix. And since I am not familiar with Rcpp and RcppArmadillo, I have no
> idea how to reuse the memory. I hope I can have some materials to learn
> this. Thanks.
>
> No, your main issue is not thinking about the computation.  As soon as you
write something like

arma::vec betahat = arma::inv(Inv)*arma::trans(D)*W*y;

you are in theory land which has very little relationship to practical
numerical linear algebra.  If you want to perform linear algebra
calculations like weighted least squares you should first take a bit of
time to learn about numerical linear algebra as opposed to theoretical
linear algebra.  They are very different disciplines.  In theoretical
linear algebra you write the solution to a system of linear equations as
above, using the inverse of the system matrix.  The first rule of numerical
linear algebra is that you never calculate the inverse of a matrix, unless
you only plan to do toy examples.  You mentioned sizes of 4000 by 4000
which means that the method you have chosen is doing thousands of times
more work than necessary (hint: how do you think that the inverse of a
matrix is calculated in practice? - ans: by solving n systems of equations,
which you are doing here when you could be solving only one).

Dirk and I wrote about 7 different methods of solving least squares
problems in our vignette on RcppEigen.  None of those methods involve
taking the inverse of an n by n matrix.

R and Rcpp and whatever other programming technologies come along will
never be a "special sauce" that takes the place of thinking about what you
are trying to do in a computation.

I could explain about using matrix decompositions, especially the Cholesky
and QR decompositions, to solve such problems but you have already ignored
all the other suggestions to think about the problem you are addressing and
decompose it into manageable chunks so I have no expectation that you would
pay attention to what I would write.


>
>> What exactly do these timings show?  A single call you your function?
>> How many calls?
>>
>> Here I called my function for 100 times.
>
>
>> Building on Romain's point: -- a portion of your function's runtime is
>> in memory allocation
>> (and you have a lot of allocations here).
>> If you're calling your function thousands or millions of times, then
>> it might pay to closely
>> examine your memory allocation strategies and figure out what's
>> temporary, for example.
>> It looks like you're already using  copy_aux_mem = false in a number
>> of places, but you're
>> allocating a lot of objects -- of approx what size?
>>
>> For example, wouldn't this work just as well with one less allocation?
>> arma::vec kk = t;
>> arma::uvec q1 = arma::find(arma::abs(tp)> kk.elem(q1) = ((1-arma::pow(tp.elem(q1)/h,2))/h)*0.75;
>> // done with q1.  let's reuse it.
>> q1 = arma::find(arma::abs(tp)>=h);
>> // was q2
>> kk.elem(q1).zeros();
>>
>> You could potentially allocate memory for temporary working space in
>> R, grab it with copy_aux_mem = false, write your temp results there,
>> and reuse these objects in subsequent function calls.  It doesn't make
>> sense to go to this trouble, though, if your core algorithm consumes
>> the bulk of runtime.
>>
>> Have you looked on the armadillo notes r.e. inv?  Matrix inversion has
>> O(>n^2).  You may be aided by pencil-and-paper math here.
>> http://arma.sourceforge.net/docs.html#inv
>>
>> Here my matrix for inverse is only 4 by 4, so I think it's ok.
>
>
>> best,
>> Christian
>>
>> > 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:
>>
>>
>>
>> --
>> A man, a plan, a cat, a ham, a yak, a yam, a hat, a canal – Panama!
>>
>
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] R vectorisation vs. C++ vectorisation

[Douglas Bates]

On Tue, Nov 20, 2012 at 9:54 AM, Hadley Wickham  wrote:

> > P.S. I don't think the sugar versions can be made any quicker, because
> > they have to allocate intermediate vectors, and do more memory copies.
>
> I don't think they _have_to - my understanding is that the expression
> template approach is general enough that you could compile away those
> intermediates.  It may be extremely difficult to do though!
>

That's the issue in delayed evaluation.  You do a lot of template magic to
discover how the result is going to be used before you evaluate it.  It has
its place but for cases like this it seems like what Ed Deming would have
called "burning the toast and then scraping it".


>
> Hadley
>
> --
> RStudio / Rice University
> http://had.co.nz/
> ___
> 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
>
___
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

Re: [Rcpp-devel] R vectorisation vs. C++ vectorisation

[Douglas Bates]

On Mon, Nov 19, 2012 at 11:36 AM, Hadley Wickham wrote:

> > To be fair with the R language, I would have compare with the use of real
> > primitives :
> >
> > vaccBubu <- function(age, female, ily) {
> > gender <- female * 1.25
> > gender[!female] <- 0.75
> > p <- (0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily) * gender
> > p[p < 0] <- 0
> > p[p > 1] <- 1
> > p
> > }
> >
>
> I think you could do it a bit more elegantly with
>
> gender <- 0.75 + female * 0.5
>
> but pmin and pmax are "real" primitives, so I think it's fine to use them.


In the folklore ifelse, pmin and pmax are know to be slow.  I know that
there was some work on pmin and pmax but I think they are still more
sluggish than desired.

In early work on generalized linear mixed models i did some profiling and
discovered that an inordinate amount of time was spent on evaluating the
inverse link, derivative and variance functions and a lot of that could be
traced to pmin and pmax.
___
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

Re: [Rcpp-devel] R vectorisation vs. C++ vectorisation

[Douglas Bates]

I'd reply to your questions that yes Sebastian Villemot and Elliot Saba are
working on the Debian packaging (last weekend there was a big commit that
reorganized the files to facilitate this) and Jeff has announced
comprehensions on distributed arrays, although I haven't tried using it yet.

However, the list administrator has told me not to discuss Julia on this
list so I won't :-)

On Mon, Nov 19, 2012 at 10:54 AM, Dirk Eddelbuettel  wrote:

>
> On 19 November 2012 at 10:47, Douglas Bates wrote:
> | Sigh.  Speaking as one of the "Julia guys" I should point out two things
> (not
> | that they will change Dirk's "cold, dead hands" attitude towards Julia
> :-)
>
> No, just "lazy apt-get-able awaiting" hands.
>
> | 1. Comprehensions provide what I feel is a clean syntax for sugar-like
> | operations in Julia
> |
> | 2. A problem with vectorization is the issue of multiple loops, hence the
> | number of attempts at implementing delayed evaluation in compiled code
> (Eigen)
> | and in add-on's to R.
> |
> | A translation of Hadley's vacc3 into Julia could be
> |
> | function vacc3a(age::Float64, female::Bool, ily::Float64){
> |   p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
> |   p *= female ? 1.25 : 0.75
> |   min(max(0., p), 1.)
> | }
> |
> | out = [vacc3a(age[i], female[i], ily[i]) for i in 1:length(age)]
> |
> | The comprehension collapses the
> |  1. Determine the length of the output vector
> |  2. Allocate the result
> |  3. Loop over indices populating the result
> |  4. Return the result
> |
> | to a single syntactic element that, in my opinion, is quite readable.
>
> Does that already give you parallel processing?  That is the one element
> missing here and I was thinking that in order completely push Hadley over
> the
> C++ cliff I should introduce him to OpenMP :)
>
> The next real trick will be to get these things done transparently on all
> cores.
>
> Dirk
>
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
>
___
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

Re: [Rcpp-devel] R vectorisation vs. C++ vectorisation

[Douglas Bates]

On Mon, Nov 19, 2012 at 9:56 AM, Dirk Eddelbuettel  wrote:

>
> On 19 November 2012 at 09:31, Hadley Wickham wrote:
> | Hi all,
> |
> | Inspired by "Rcpp is smoking fast for agent-based models in data
> | frames" (http://www.babelgraph.org/wp/?p=358), I've been doing some
>
> [ I liked that post, but we got flak afterwards as his example was not well
> chosen. The illustration of the language speed difference does of course
> hold. ]
>
> | exploration of vectorisation in R vs C++ at
> | https://gist.github.com/4111256
> |
> | I have five versions of the basic vaccinate function:
> |
> | * vacc1: vectorisation in R with a for loop
> | * vacc2: used vectorised R primitives
> | * vacc3: vectorised with loop in C++
> | * vacc4: vectorised with Rcpp sugar
> | * vacc5: vectorised with Rcpp sugar, explicitly labelled as containing
> | no missing values
> |
> | And the timings I get are as follows:
> |
> | Unit: microseconds
> | exprmin lq median uq max neval
> |  vacc1(age, female, ily) 6816.8 7139.4 7285.7 7823.9 10055.5   100
> |  vacc2(age, female, ily)  194.5  202.6  212.6  227.9   260.4   100
> |  vacc3(age, female, ily)   21.8   22.4   23.4   24.935.5   100
> |  vacc4(age, female, ily)   36.2   38.7   41.3   44.555.6   100
> |  vacc5(age, female, ily)   29.3   31.3   34.0   36.452.1   100
> |
> | Unsurprisingly the R loop (vacc1) is very slow, and proper
> | vectorisation speeds it up immensely.  Interestingly, however, the C++
> | loop still does considerably better (about 10x faster) - I'm not sure
> | exactly why this is the case, but I suspect it may be because it
> | avoids the many intermediate vectors that R requires.  The sugar
> | version is about half as fast, but this gets quite a bit faster with
> | explicit no missing flags.
> |
> | I'd love any feedback on my code 
> (https://gist.github.com
>


> /4111256 ) -
> | please let me know if I've missed anything obvious.
>
> I don't have a problem with sugar being a little slower that hand-rolling.
> The code is so much simpler and shorter. And we're still way faster than
> vectorised R.  I like that place.
>
> Somewhat off-topic/on-topic: I am still puzzled by how the Julia guys now
> revert back from vectorised code to hand-written loops because llvm does
> better on those.  Speed is good, but concise code with speed is better in
> my
> book.
>

Sigh.  Speaking as one of the "Julia guys" I should point out two things
(not that they will change Dirk's "cold, dead hands" attitude towards Julia
:-)

1. Comprehensions provide what I feel is a clean syntax for sugar-like
operations in Julia

2. A problem with vectorization is the issue of multiple loops, hence the
number of attempts at implementing delayed evaluation in compiled code
(Eigen) and in add-on's to R.

A translation of Hadley's vacc3 into Julia could be

function vacc3a(age::Float64, female::Bool, ily::Float64){
  p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
  p *= female ? 1.25 : 0.75
  min(max(0., p), 1.)
}

out = [vacc3a(age[i], female[i], ily[i]) for i in 1:length(age)]

The comprehension collapses the
 1. Determine the length of the output vector
 2. Allocate the result
 3. Loop over indices populating the result
 4. Return the result

to a single syntactic element that, in my opinion, is quite readable.


> Hence I would prefer to invoke the 80/20 rule as I think we have better
> targets to chase than to narrow that gap. But that's just my $0.02...
>
> If you can't sleep til both version have 20-some microsend medians then by
> all means go crazy ;-)
>
> Dirk
>
>
>
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> ___
> 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
>
___
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

Re: [Rcpp-devel] Matrix: column-major storage?

[Douglas Bates]

An Rcpp object uses the storage allocated in R hence it has the same
column-major ordering of the elements in a matrix.  The same is true of
matrix objects in RcppArmadillo.

Basically anything that uses Lapack and BLAS or calling sequences
compatible with Lapack or BLAS uses column-major ordering.

With RcppEigen the default is column-major but you can specify row-major.
 However, when wrapping the Eigen object to return it to R there will be a
copying operation to convert row major to column major.

You are probably safest staying with column major ordering if possible.


On Wed, Nov 7, 2012 at 1:25 PM, Ramon Diaz-Uriarte wrote:

>
> Dear All,
>
> I've recently started using Rcpp, so maybe this is a trivial question, but
> I have not been able to find an answer in previous posts or the docs.
>
>
> I thought that, for any of the Matrix types available, if we have
>
> x(i, j)
>
> then
>
> x(i, j + 1)
>
> would be contiguous in memory (whereas x(i + 1, j) could be extremely far
> away, depending on the size of the matrix).
>
> However, it seems that with Matrix objects it is the other way around, so
> that x(i, j) is next to x(i + 1, j), not to x(i, j + 1):
>
> #
>
> showThePointers <- cxxfunction(signature(dm = "integer"),
>  body = '
>
> int m = as(dm);
> IntegerMatrix xm(m, m);
> int *p;
> int i;
> int j;
>
> std::cout << std::endl << "* Row-major access *" << std::endl;
>
> for(i = 0; i < m; i++) {
>for(j = 0; j < m; j++) {
>  xm(i, j) = i * m + j;
>  p = &xm(i, j);
>  std::cout << "i = " << i << "; j = " <<
>  j << "; xm(i, j) = " << xm(i, j) << "; p " << p << std::endl;
>   }
> }
>
> std::cout << std::endl << " ** Column-major access *" << std::endl;
>
> for(j = 0; j < m; j++) {
>for(i = 0; i < m; i++) {
>  xm(i, j) = i * m + j;
>  p = &xm(i, j);
>  std::cout << "i = " << i << "; j = " <<
>  j << "; xm(i, j) = " << xm(i, j) << "; p " << p << std::endl;
>   }
> }
>
> ',
>  plugin = "Rcpp",
>  verbose = TRUE)
>
> showThePointers(3)
>
> #
>
>
>
>
> I noticed that because if we time functions that access and modify in
> row-major order, they are a lot slower than in column-major order (in my
> laptop, about four times slower in the examples below, about nine times if
> I then return the object to R).
>
>
> #
>
> byRow <- cxxfunction(signature(dm = "integer"),
>  body = '
>
> int m = as(dm);
> IntegerMatrix xm(m, m);
> int i;
> int j;
> for(i = 0; i < m; i++) {
>for(j = 0; j < m; j++) {
>  xm(i, j) = i * m + j;
>   }
> }
> ',
>  plugin = "Rcpp",
>  verbose = TRUE)
>
> byCol <- cxxfunction(signature(dm = "integer"),
>  body = '
>
> int m = as(dm);
> IntegerMatrix xm(m, m);
> int i;
> int j;
> for(j = 0; j < m; j++) {
>for(i = 0; i < m; i++) {
>  xm(i, j) = i * m + j;
>   }
> }
> ',
>  plugin = "Rcpp",
>  verbose = TRUE)
>
>
> nn <- 1
> benchmark(byRow(nn), byCol(nn),
>   columns=c("test", "replications", "elapsed", "relative",
>   "user.self", "sys.self"),
>   order="relative", replications=10)
>
>
> #
>
>
> I understand it is not possible to store things in row-major order with
> Rcpp Matrix types?
>
>
> Best,
>
>
> R.
>
>
> --
> Ramon Diaz-Uriarte
> Department of Biochemistry, Lab B-25
> Facultad de Medicina
> Universidad Autónoma de Madrid
> Arzobispo Morcillo, 4
> 28029 Madrid
> Spain
>
> Phone: +34-91-497-2412
>
> Email: rdia...@gmail.com
>ramon.d...@iib.uam.es
>
> http://ligarto.org/rdiaz
>
> ___
> 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
___
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

Re: [Rcpp-devel] Vectorized version of rbinom function via Rcpp

[Douglas Bates]

On Thu, Oct 25, 2012 at 9:29 AM, Dirk Eddelbuettel  wrote:
>
> Hi Carsten,
>
> On 25 October 2012 at 16:05, carsten Gerken wrote:
> | Hi All,
> |
> | I am trying to outsource the simulation of binomially distributed random
> | variables to C++ via Rcpp (e.g. to replicate the R command 
> rbinom(10,10,c(1:10)
> | /20)). This works well as long as the third argument of the function (i.e. 
> the
> | probability) is a single value and not a vector. I tried to compile the
> | following code:
> |
>
> [ Some editing here ]
>
>  rbinom_c <- rcpp( signature(n = "integer",m = "integer",p="numeric"), body="
>int N = as(n);
>int M = as(m);
>NumericVector P = NumericVector(p);
>NumericVector res = rbinom(N,M,P);
>return res;
>  ")
>
> | But the compiler expects a “double” variable in the third argument.
>
> Correct, as I recall.  Some (common) things are vectors, some are scalars.
> As a simple example in rnorm(n, m, s) the mean m and stdev s are expected to
> be scalars.
>
> | Is there a
> | way of fixing the code without using an explicit loop in the C++ code?
>
> I don't think so. The code currently has this hard-wired:
>
> inline NumericVector rbinom( int n, double nin, double pp ){
> return NumericVector( n, ::Rf_rbinom, nin, pp) ;
> }
>
> so until one of us has sufficient to extend this (and test it) will have to
> do the loop at your end.

And I suggest that before expending a lot of effort on this you check
whether you are likely to do better than the R function rbinom.  That
R function is vectorized in C code and you are unlikely to make things
much faster by using Rcpp, unless you want to use the random variates
in another C++ function.
___
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

Re: [Rcpp-devel] Slices of 3d arrays as matrices

[Douglas Bates]

On Wed, Oct 24, 2012 at 11:14 AM, Giovanni Petris  wrote:
> Hi Doug,
>
> Thank you for your answer and for your patience reading my confusing code. I 
> have read the jss paper on RcppEigen and the Eigen tutorial. I am very 
> excited about the opportunities that RcppEigen (and Rcpp) gives to the R 
> programmer.  Thank you so much for the very useful contribution!
>
> To get back quickly to the point in my post, the reason for accessing and 
> modifying an R object from within C++ code is that the final application will 
> involve large matrices and many iterations of a Gibbs sampler. Therefore, I 
> would like to have some workspace allocated once and for all at the beginning 
> of the calculations. In C I would probably declare a static double *workmat, 
> invisible from R, but permanent for the different calls to my C function.

I know the problem.  A lot of the effort in the RcppEigen-based C++
code for the lme4 package is designed to allow for in-place updating
of large structures.

One approach that is legitimate in R, but only because of special
semantics for environments, is to make a new environment and install
the objects that you want to modify in there.  Environments behave
differently from other R objects in that there is only ever one copy
of an environment so changes in the objects in an environment made in
one place are reflected in all other views of the environment.

> From: dmba...@gmail.com [dmba...@gmail.com] on behalf of Douglas Bates 
> [ba...@stat.wisc.edu]
> Sent: Wednesday, October 24, 2012 10:34 AM
> To: Giovanni Petris
> Cc: rcpp-devel@lists.r-forge.r-project.org
> Subject: Re: [Rcpp-devel] Slices of 3d arrays as matrices
>
> On Tue, Oct 23, 2012 at 5:52 PM, Giovanni Petris  wrote:
>> Thank you Doug and Dirk.
>>
>> I was afraid I had to move to Eigen... The problem now is that in the 
>> transition I am working with matrices both of type NumericMatrix and 
>> Map.
>>
>> The following snippet is a loop that is giving me troubles...
>>
>>wz1 <- lapply(1:K, function(k) matrix(0.0, N, M)) # this is defined in R, 
>> but used by the C++ function
>>
>> Rcpp::NumericMatrix cpp_u(u);
>> Rcpp::NumericMatrix cpp_mu(mu);
>> Rcpp::NumericVector cpp_Sigma(Sigma);
>> double *mm = cpp_Sigma.begin();
>> int K = cpp_mu.ncol();
>> int M = cpp_mu.nrow();
>> int N = cpp_u.nrow();
>> Mmat var1(NULL, M, M);
>> SelfAdjointEigenSolver var1_eigen(M);
>> Environment glob = Environment::global_env();
>> List cpp_wz1(glob.get("wz1"));
>> NumericMatrix tmpMatrix1;
>>
>>  for (int k=0; k < K; k++) {
>> tmpMatrix1 = as(cpp_wz1[k]);
>> for (int i=0; i < M; i++)
>> tmpMatrix1(_, i) = cpp_u(_, i) - cpp_mu(i, k);
>> new (&var1)Mmat(mm + k * M * M, M, M);
>> var1_eigen.compute(var1);
>> tmpMatrix1 = tmpMatrix1 * var1_eigen.operatorInverseSqrt(); // <== 
>> here is the problem
>>// more stuff here...
>> }
>
>> Do you have any suggestion for an analog of the line
>
>>   tmpMatrix1 = as(cpp_wz1[k]);
>
>> using an Eigen::Map object instead of the 
>> Rcpp::NumericMatrix object tmpMatrix1?
>
> This code is very confusing and breaks the fundamental rule of
> functional semantics in R.  It is very dangerous to reach out into the
> global environment, grab the value of a symbol and change it.  R does
> copy on write but this code changes the wz1 object without R being
> aware that it needs to copy it.
>
> In other words, you are going about things in an awkward way.  You
> have the ability to create the object to be returned within the C++
> code so this is a roundabout and dangerous way of creating the result.
>
> To stay within the Eigen framework, I would create the array to be
> returned as a matrix of dimension (M, M*K) then change the dimension
> in R to (M, M, K) upon return.
>
> There are several idioms for dealing with matrices passed from R as
> Eigen mapped matrices in C++ as Dirk and I explained in our vignette
> on RcppEigen.  It would be worthwhile reading that vignette if you
> want to use RcppEigen, and also read at least the tutorial in the
> Eigen3 documentation.
>
> I would start with something like the enclosed but that code doesn't
> compile and gives the usual page after page of reports when compiling
> Eigen-based C++ code with g++.  Regrettably I don't have time to debug
> it.
>
>> From: dmba...@gmail.com [dmba...@gmail.com] on behalf of Douglas Bates 
>> [ba...@stat.wisc.edu]
>> Sent: Tuesday, October 23, 2012 1:20 PM
>> To: Giovanni Pet

Re: [Rcpp-devel] R/Rcpp/RcppEigen Optimization WAS: NumericVector Double mismatch when indexing an array

[Douglas Bates]

On Tue, Oct 23, 2012 at 6:01 PM, Darren Cook  wrote:
>> What gives?  What I found that does speed up the code dramatically
>> is the -march switch.  I guess that can't be repo-ed because it is
>> CPU dependent, right?  Here's the important settings that I used to
>> compile R from source:
>>
>> CC="gcc"
>> CFLAGS="-g -O2 -march=amdfam10"
>> CXX="g++"
>> CXXFLAGS="-g -O2 -march=amdfam10"
>>
>> With these settings...
>>
>> MCMC CodeSwitchesRelative Time
>> Rnone or -O21.0
>> RcppEigennone or -O20.09
>> R-O2 -march=amdfam100.5
>> RcppEigen-O2 -march=amdfam100.013
>
> Wow, those are huge differences. Am I misreading, or does that say
> RcppEigen runs *seven* times quicker with the -march=amdfam10 option?
>
> Can that be explained? (e.g. does the AMDFAM10 processor have some
> feature that speeds it up 7 times in ideal conditions?)

Eigen does use SSE2 and SSE3 pipelined instructions when they are
determined to be available.
___
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

Re: [Rcpp-devel] Slices of 3d arrays as matrices

[Douglas Bates]

On Tue, Oct 23, 2012 at 5:52 PM, Giovanni Petris  wrote:
> Thank you Doug and Dirk.
>
> I was afraid I had to move to Eigen... The problem now is that in the 
> transition I am working with matrices both of type NumericMatrix and 
> Map.
>
> The following snippet is a loop that is giving me troubles...
>
>wz1 <- lapply(1:K, function(k) matrix(0.0, N, M)) # this is defined in R, 
> but used by the C++ function
>
> Rcpp::NumericMatrix cpp_u(u);
> Rcpp::NumericMatrix cpp_mu(mu);
> Rcpp::NumericVector cpp_Sigma(Sigma);
> double *mm = cpp_Sigma.begin();
> int K = cpp_mu.ncol();
> int M = cpp_mu.nrow();
> int N = cpp_u.nrow();
> Mmat var1(NULL, M, M);
> SelfAdjointEigenSolver var1_eigen(M);
> Environment glob = Environment::global_env();
> List cpp_wz1(glob.get("wz1"));
> NumericMatrix tmpMatrix1;
>
>  for (int k=0; k < K; k++) {
> tmpMatrix1 = as(cpp_wz1[k]);
> for (int i=0; i < M; i++)
> tmpMatrix1(_, i) = cpp_u(_, i) - cpp_mu(i, k);
> new (&var1)Mmat(mm + k * M * M, M, M);
> var1_eigen.compute(var1);
> tmpMatrix1 = tmpMatrix1 * var1_eigen.operatorInverseSqrt(); // <== 
> here is the problem
>// more stuff here...
> }

> Do you have any suggestion for an analog of the line

>   tmpMatrix1 = as(cpp_wz1[k]);

> using an Eigen::Map object instead of the Rcpp::NumericMatrix 
> object tmpMatrix1?

This code is very confusing and breaks the fundamental rule of
functional semantics in R.  It is very dangerous to reach out into the
global environment, grab the value of a symbol and change it.  R does
copy on write but this code changes the wz1 object without R being
aware that it needs to copy it.

In other words, you are going about things in an awkward way.  You
have the ability to create the object to be returned within the C++
code so this is a roundabout and dangerous way of creating the result.

To stay within the Eigen framework, I would create the array to be
returned as a matrix of dimension (M, M*K) then change the dimension
in R to (M, M, K) upon return.

There are several idioms for dealing with matrices passed from R as
Eigen mapped matrices in C++ as Dirk and I explained in our vignette
on RcppEigen.  It would be worthwhile reading that vignette if you
want to use RcppEigen, and also read at least the tutorial in the
Eigen3 documentation.

I would start with something like the enclosed but that code doesn't
compile and gives the usual page after page of reports when compiling
Eigen-based C++ code with g++.  Regrettably I don't have time to debug
it.

> From: dmba...@gmail.com [dmba...@gmail.com] on behalf of Douglas Bates 
> [ba...@stat.wisc.edu]
> Sent: Tuesday, October 23, 2012 1:20 PM
> To: Giovanni Petris
> Cc: rcpp-devel@lists.r-forge.r-project.org
> Subject: Re: [Rcpp-devel] Slices of 3d arrays as matrices
>
> On Tue, Oct 23, 2012 at 12:56 PM, Giovanni Petris  wrote:
>> Hello,
>>
>> I have a 3d array and I want to compute the eigendecomposition of each slice 
>> [,, k]. I am trying to make Rcpp see the relevant slice of the 3d array as a 
>> matrix. My experiments, along the lines illustrated below, have so far been 
>> unsuccessful. Does anybody have a suggestion? Thank you in advance!
>
>>NumericVector cpp_Sigma(Sigma); // 3d array, with dim(Sigma) = c(M,M,K), 
>> as a vector
>>NumericMatrix tmpMatrix;
>>List var_eigen;
>>Function eigen("eigen");
>>int k=2; // 3rd slice, for example
>>
>>tmpMatrix = cpp_Sigma[k * M * M];
>>tmpMatrix.attr("dim") = Dimension(M, M);
>>var_eigen = eigen(tmpMatrix);
>>return var_eigen;
>
> To begin with, it is not a good idea to call R's eigen function from
> within C++ code.  It would be much better to call code from the
> RcppEigen or RcppArmadillo packages.  The eigen function in R is just
> going to marshal its arguments and call a Lapack routine, which is
> done more cleanly in RcppArmadillo.
>
> Are the slices symmetric matrices?  If they are general matrices you
> may need to allow for complex eigenvalues and eigenvectors.
>
> Assuming that the slices are symmetric you could use something like
>
> typedef Eigen::Map Mmat;
> const NumericVector cpp_Sigma(Sigma);
> double *mm = cpp_Sigma.begin();
> int k = 2;
> const Mmat em(mm + k * M * M, M, M);
> const Eigen::SelfAdjointEigenSolver VLV(em);
> return wrap(VLV.eigenvalues());
>
> if you just want the eigenvalues.


foo.R
Description: Binary data
___
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

Re: [Rcpp-devel] Slices of 3d arrays as matrices

[Douglas Bates]

On Tue, Oct 23, 2012 at 12:56 PM, Giovanni Petris  wrote:
> Hello,
>
> I have a 3d array and I want to compute the eigendecomposition of each slice 
> [,, k]. I am trying to make Rcpp see the relevant slice of the 3d array as a 
> matrix. My experiments, along the lines illustrated below, have so far been 
> unsuccessful. Does anybody have a suggestion? Thank you in advance!

>NumericVector cpp_Sigma(Sigma); // 3d array, with dim(Sigma) = c(M,M,K), 
> as a vector
>NumericMatrix tmpMatrix;
>List var_eigen;
>Function eigen("eigen");
>int k=2; // 3rd slice, for example
>
>tmpMatrix = cpp_Sigma[k * M * M];
>tmpMatrix.attr("dim") = Dimension(M, M);
>var_eigen = eigen(tmpMatrix);
>return var_eigen;

To begin with, it is not a good idea to call R's eigen function from
within C++ code.  It would be much better to call code from the
RcppEigen or RcppArmadillo packages.  The eigen function in R is just
going to marshal its arguments and call a Lapack routine, which is
done more cleanly in RcppArmadillo.

Are the slices symmetric matrices?  If they are general matrices you
may need to allow for complex eigenvalues and eigenvectors.

Assuming that the slices are symmetric you could use something like

typedef Eigen::Map Mmat;
const NumericVector cpp_Sigma(Sigma);
double *mm = cpp_Sigma.begin();
int k = 2;
const Mmat em(mm + k * M * M, M, M);
const Eigen::SelfAdjointEigenSolver VLV(em);
return wrap(VLV.eigenvalues());

if you just want the eigenvalues.
___
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

Re: [Rcpp-devel] 'Nested' Rcpp functions using inline

[Douglas Bates]

Remember that the inline package is convenient for interactive
exploration during code development but usually the ultimate goal when
working with Rcpp is to produce a package including both R and C++
code.  When you start doing complicated things in the C++ code it is
probably time to look at creating a package.

On Wed, Oct 10, 2012 at 9:59 AM, Jeffrey Pollock  wrote:
> I think you need to use the `includes` argument of cxxfunction to include a
> pure c++ function ie;
>
> library(inline)
> library(Rcpp)
>
> inc <-'
> int fun1(double x, double y) {
> return (exp(x) - y) > 0 ? 1 : 0;
> }
> '
>
> body <-'
>
> NumericVector u   = as(u_r);
> doublemu  = as(mu_r),
> v   = log(mu);
> int   n   = u.size(),
> res = 0;
>
> for(int i=0; i res += fun1(u(i), v);
> }
>
> return wrap(res);
> '
>
> fun.test <- cxxfunction(signature(u_r = 'numeric', mu_r = 'numeric'), body,
> "Rcpp", inc)
>
>> identical(fun.test(1:4, 100), sum(exp(1:4) > log(100)))
> [1] TRUE
>
> On Wed, Oct 10, 2012 at 3:39 PM, Rubem Kaipper Ceratti
>  wrote:
>>
>> Hi,
>>
>> I'm currently using Rcpp and inline to speed up some R code I wrote, but
>> being a C++/Rcpp newbie I'm running into some difficulties. More
>> specifically, I'm having problems calling an Rcpp function inside another
>> one. As a toy example consider the code:
>> ##
>> library(Rcpp)
>> library(inline)
>>
>> code2.0 <-'
>>   NumericVector u   = as(u_r);
>>   doublemu  = as(mu_r),
>> v   = log(mu);
>>   int   n   = u.size(),
>> res = 0;
>>
>>   for(int i=0; i> res += (exp(u(i))-v) > 0 ? 1 : 0;
>>   }
>>
>>   return wrap(res);
>> '
>> sig2 <- signature(u_r = 'numeric', mu_r = 'numeric')
>> fun.test <- cxxfunction(sig2, code2.0, plugin = "Rcpp")
>>
>> fun.test(1:4, 100)
>> sum(exp(1:4) > log(100))  # ok!
>> ##
>>
>> It works just fine. Now, I'd like to break it down into two functions like
>> so:
>> ##
>> code1 <-'
>>   double x = as(x_r),
>>  y = as(y_r),
>>  ans;
>>
>>   ans = (exp(x)-y) > 0 ? 1 : 0;
>>
>>   return wrap(ans);
>> '
>>
>> code2.1 <-'
>>   NumericVector u   = as(u_r);
>>   doublemu  = as(mu_r),
>> v   = log(mu);
>>   int   n   = u.size(),
>> res = 0;
>>
>>   for(int i=0; i> res += as(fun1(u(i),v));
>>   }
>>
>>   return wrap(res);
>> '
>>
>> sig1 <- signature(x_r = 'numeric', y_r = 'numeric')
>> sig2 <- signature(u_r = 'numeric', mu_r = 'numeric')
>> fun.test <- cxxfunction(list(fun1 = sig1, fun2 = sig2),
>> list(code1, code2.1), plugin = "Rcpp")
>> ##
>>
>> But I get the error message:
>> # file5686a0c71e8.cpp: In function 'SEXPREC* fun2(SEXP, SEXP)':
>> #   file5686a0c71e8.cpp:52:33: error: invalid cast from type 'double' to
>> type 'SEXP'
>> # file5686a0c71e8.cpp:55:20: error: invalid cast from type
>> 'Rcpp::traits::storage_type<14>::type {aka double}' to type 'SEXP'
>>
>> I then tried to change the types of 'u(i)' and 'v' to SEXP, but I get the
>> same error:
>> ##
>> code2.2 <-'
>>   NumericVector u   = as(u_r);
>>   doublemu  = as(mu_r),
>> v   = log(mu);
>>   int   n   = u.size(),
>> res = 0;
>>   SEXP  us, vs = (SEXP) v;
>>
>>   for(int i=0; i> us = (SEXP) u(i);
>> res += as(fun1(us,vs));
>>   }
>>
>>   return wrap(res);
>> '
>>
>> fun.test <- cxxfunction(list(fun1 = sig1, fun2 = sig2),
>> list(code1, code2.2), plugin = "Rcpp")
>> ##
>> # file5686a0c71e8.cpp: In function 'SEXPREC* fun2(SEXP, SEXP)':
>> #   file5686a0c71e8.cpp:52:33: error: invalid cast from type 'double' to
>> type 'SEXP'
>> # file5686a0c71e8.cpp:55:20: error: invalid cast from type
>> 'Rcpp::traits::storage_type<14>::type {aka double}' to type 'SEXP'
>>
>> Is there any way to make it work?
>>
>> Thanks,
>> Rubem
>>
>>
>>
>> ___
>> 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
>
>
>
> ___
> 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
___
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

Re: [Rcpp-devel] Correlation in ArmadilloRcpp with missing values (nan)

[Douglas Bates]

On Tue, Oct 9, 2012 at 1:06 PM, mateusz.ka...@gmail.com
 wrote:
> Maybe I am not good in coding, but even though I managed to make uvec with
> ones and zeros, I cannot use it to select rows  in given column, compiler
> complains
> error: no matching function for call to
> ‘arma::Mat::submat(arma::uvec&, int&)’

Basically your problem comes down to looking at each pair of columns
and determining the joint missingness pattern.  You can try to force
the calculation into arma but I think you are better off just writing
the code in loops.  To be careful you should calculate the means of
each column on the pairwise non-missing values as in the enclosed
function.

> On 9 October 2012 19:55, Dirk Eddelbuettel  wrote:
>>
>>
>> On 9 October 2012 at 19:07, mateusz.ka...@gmail.com wrote:
>> | Can you provide an example how to convert Armadillo colvec to uvec
>> vector which
>> | I assume works as selector for rows?
>>
>> You may need to loop (or use STL iterators) to fill the uvec position by
>> position.  Then use the subset, and proceed with your correlation
>> calculation.  I don't think there is a shortcut.
>>
>> Dirk
>>
>> --
>> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
>
>


foo.R
Description: Binary data
___
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

Re: [Rcpp-devel] Correlation in ArmadilloRcpp with missing values (nan)

[Douglas Bates]

Just for the record, earlier Dirk mentioned functions R_IsNA, R_IsNaN
and R_IsFinite.  The last one should be R_finite.
___
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

Re: [Rcpp-devel] Correlation in ArmadilloRcpp with missing values (nan)

[Douglas Bates]

By the way, how are you calculating this correlation in R?  Are you
using the cor function?

I'm confused because the cor function in R does a bit of bookkeeping
then calls a C function using .Internal.  It seems unlikely to me that
one could make a C++/Rcpp function run much faster.

On Tue, Oct 9, 2012 at 12:07 PM, mateusz.ka...@gmail.com
 wrote:
> Can you provide an example how to convert Armadillo colvec to uvec vector
> which I assume works as selector for rows?
>
> Thanks
>
> On 9 October 2012 18:17, mateusz.ka...@gmail.com 
> wrote:
>>
>> Dear Dirk,
>>
>> I dont see how to do that in Armadillo, but I think I can create same size
>> NumericMatrix B = is.na(X);
>> This maybe I could use for indexing ? To save time, and not to introduce
>> extra looping.
>>
>> Also, I want to perform regression column X on Z, and column Y on Z.
>> Does arma::solve(...) handle nan values ? or I guess I have to check that
>> myself ?
>>
>> The function works nicely, and while R implementation takes at least 40min
>> (surely more because I terminated), Armadillo computes everything in less
>> then one minute. But I skipped columns with missing cells, which now I want
>> to include.
>>
>> Thanks,
>> Mateusz
>>
>>
>> On 9 October 2012 17:37, Dirk Eddelbuettel  wrote:
>>>
>>>
>>> On 9 October 2012 at 10:08, Douglas Bates wrote:
>>> | You may find it easier to use the Rcpp class NumericMatrix than to use
>>> | RcppArmadillo.  Detection of NA's is built in to R and Rcpp classes but
>>> not
>>> | RcppArmadillo. For each pair of columns, run a loop that checks for
>>> NA's at
>>> | each position in each column, skips the position if NA's are detected
>>> and
>>> | otherwise increments the squared sums, cross-product and number of
>>> elements.
>>>
>>> All true, but at the same time, once in RcppArmadillo or RcppEigen ... it
>>> is
>>> convenient to stay there.
>>>
>>> You should be able to set up a little "sweeper" function which applies
>>> one of
>>>
>>> R_IsNA  just NA
>>> R_IsNaN just NaN
>>> R_IsFinite  NA, NaN or Inf
>>>
>>> across a vector or matrix and returns you an index vector. Armadillo can
>>> use
>>> indexing vectors in ways that are similar in R.
>>>
>>> Dirk
>>>
>>>
>>> |
>>> | On Oct 9, 2012 9:53 AM, "mateusz.ka...@gmail.com"
>>> 
>>> | wrote:
>>> |
>>> | Hi,
>>> |
>>> | I have written small code in C++ using Armadillo and inline with
>>> | RcppArmadillo package.
>>> | The input is data.marix(X). Some cells might be NAs. Example in R:
>>> X =
>>> | matrix(sample(c(rnorm(10*9.9),NA)),ncol=10)
>>> |
>>> | I am calculating conditional correlation on columns of that matrix,
>>> just
>>> | picking vectors, so cor(X,Y).
>>> | The problem is that sometimes I might have empty cell in one or
>>> both
>>> | vectors, in that case I would like to skip that row, and procede
>>> with
>>> | calculating Pearson's correlation on remaining data. I know that
>>> there will
>>> | be difference in degrees of freedom, but I have over 100 rows, so
>>> skiping
>>> | few shouldnt matter that much.
>>> |
>>> | Basically my question boils down to solving the problem:
>>> | How to find which colvec cells are nan, and remove this index from
>>> both X
>>> | and Y colvec, before calculating correlation.
>>> |
>>> | I would be very grateful for help,
>>> |
>>> | Kind regards,
>>> | Mateusz Kaduk
>>> |
>>> | ___
>>> | 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
>>> |
>>> |
>>> | --
>>> | ___
>>> | 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
>>> --
>>> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
>>
>>
>
___
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

Re: [Rcpp-devel] Correlation in ArmadilloRcpp with missing values (nan)

[Douglas Bates]

You may find it easier to use the Rcpp class NumericMatrix than to use
RcppArmadillo.  Detection of NA's is built in to R and Rcpp classes but not
RcppArmadillo. For each pair of columns, run a loop that checks for NA's at
each position in each column, skips the position if NA's are detected and
otherwise increments the squared sums, cross-product and number of
elements.
On Oct 9, 2012 9:53 AM, "mateusz.ka...@gmail.com" 
wrote:

> Hi,
>
> I have written small code in C++ using Armadillo and inline with
> RcppArmadillo package.
> The input is data.marix(X). Some cells might be NAs. Example in R: X =
> matrix(sample(c(rnorm(10*9.9),NA)),ncol=10)
>
> I am calculating conditional correlation on columns of that matrix, just
> picking vectors, so cor(X,Y).
> The problem is that sometimes I might have empty cell in one or both
> vectors, in that case I would like to skip that row, and procede with
> calculating Pearson's correlation on remaining data. I know that there will
> be difference in degrees of freedom, but I have over 100 rows, so skiping
> few shouldnt matter that much.
>
> Basically my question boils down to solving the problem:
> How to find which colvec cells are nan, and remove this index from both X
> and Y colvec, before calculating correlation.
>
> I would be very grateful for help,
>
> Kind regards,
> Mateusz Kaduk
>
> ___
> 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
>
___
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

Re: [Rcpp-devel] How to use R function mvrnorm in c++?

[Douglas Bates]

On Mon, Oct 8, 2012 at 4:23 PM, Victor  wrote:
> Dear all,
>
> I want to use R function mvrnorm (package MASS) in my c++ scripts. Below is
> the R code:
>
> library(MASS)
>
> x<-matrix(c(1.5,0,0,0,1.5,0,0,0,1.5),3,3)
>
> mvrnorm(1,c(15,33,26),x)
>
>
> I will appreciate you very much if you give me some help and advice! Thank
> you so much!

Read the source code for the mvrnorm function.  You will find that it
is implemented by taking an eigenvalue-eigenvector decomposition of
the Sigma matrix.  If your Sigma matrix is fixed, then simply take the
decomposition in R and pass the components.  If it varies then use
RcppEigen or RcppArmadillo to obtain the decomposition.

This is not difficult but also not trivial.  You will need to do some
reading and experimenting.  Ask yourself if the amount of time you
will save in execution of the function warrants the time you will need
to spend to learn the necessary programming.
___
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

Re: [Rcpp-devel] Missing PROTECT()s in forward_exception_to_r()?

[Douglas Bates]

On Mon, Oct 8, 2012 at 6:42 AM, Dirk Eddelbuettel  wrote:
>
> Hi Ben,
>
> On 8 October 2012 at 10:24, North, Ben wrote:
> | Hi,
> |
> | I've been using Rcpp for a while now, and finding it very useful ---
> | thanks!
> |
> | Recently, though, I came across strange behaviour when running under
> | gctorture(TRUE) and C++ exceptions were being forwarded to R.  The error
> | message string at the R level was some essentially random pathname, and
> | the string in the C++ exception was nowhere to be seen.  The classname
> | at the R level was the logical TRUE value.
> |
> | The patch below fixed the problem for me.  I would not claim to fully
> | know what I'm doing, so perhaps somebody who does could review and then
> | make the proper fix :-)
>
> It's been a while since we had something like this, and Doug lead it the last
> time. From what I recall, Rf_install() is safe, and the reason we factor
> these out as below. The bad boy, as I call, is probably Rf_mkString so by 
> just doing
>
>  str_what  = Rf_install(exception_what.c_str());
>  str_class = Rf_install(exception_class.c_str());
>  str_ns= Rf_install("Rcpp"));
>  Rf_eval(Rf_lang3(cppExceptSym, str_what, str_class), str_ns);
>
> may be enough.  Could you test that?

Just glancing quickly at this, I think it may or may not work.
Rf_install looks up a symbol in the symbol table, installing it if
needed, and returns the symbol SEXP.  Rf_mkstring creates a string
SEXP (vector of R class "character" containing a single string).  The
string probably needs protection.  The symbol does not.

I imagine the R function cpp_exception expects strings.  I'm not sure
what happens if it gets symbols.

> |
> | Thanks,
> |
> | Ben.
> |
> | --- src/exceptions.cpp   2012-10-01 01:29:53.0 +0100
> | +++ ../Rcpp/src/exceptions.cpp 2012-10-08 11:14:39.435644100 +0100
> | @@ -161,14 +161,23 @@
> |exception_class = name ; /* just using the mangled name */
> |   }
> |  }
> | -SEXP cppExceptSym = Rf_install("cpp_exception"); // cannot cause a 
> gc() once in symbol table
> | -Rf_eval(
> | -Rf_lang3(
> | - cppExceptSym,
> | -  Rf_mkString(exception_what.c_str()),
> | -  Rf_mkString(exception_class.c_str())
> | - ), R_FindNamespace(Rf_mkString("Rcpp"))
> | -) ;
> | +
> | +SEXP cppExceptSym;
> | +SEXP str_what, str_class, str_ns, ns, expr;
> | +
> | +PROTECT(cppExceptSym = Rf_install("cpp_exception"));
> | +PROTECT(str_what = Rf_mkString(exception_what.c_str()));
> | +PROTECT(str_class = Rf_mkString(exception_class.c_str()));
> | +PROTECT(str_ns = Rf_mkString("Rcpp"));
> | +
> | +PROTECT(ns = R_FindNamespace(str_ns));
> | +PROTECT(expr = Rf_lang3(cppExceptSym, str_what, str_class));
> | +
> | +Rf_eval(expr, ns); /* Should not return. */
> | +
> | +// Do this in case somehow someone replaces the definition of
> | +// "cpp_exception" such that it does return:
> | +UNPROTECT(6);
> |  }
> |  #else
> |
> | 
> |
> | Susquehanna International Group Limited
> |
> | Susquehanna International Group Limited is a private company limited by 
> shares and registered in Ireland. Registration No.: 445356. Registered 
> Address: 4th Floor, Georges Dock House, IFSC, Dublin 1, Ireland.
> | Susquehanna Ireland Limited is a private company limited by shares and 
> registered in Ireland. Registration No.: 305632. Registered Address: 4th 
> Floor, Georges Dock House, IFSC, Dublin 1, Ireland.
> | Susquehanna International Securities Limited is a private company limited 
> by shares and registered in Ireland. Registration No.: 337946. Registered 
> Address: 4th Floor, Georges Dock House, IFSC, Dublin 1, Ireland.
> |
> | IMPORTANT: The information contained in this email and/or its attachments 
> is confidential. If you are not the intended recipient, please notify the 
> sender immediately by reply and immediately delete this message and all its 
> attachments. Any review, use, reproduction, disclosure or dissemination of 
> this message or any attachment by an unintended recipient is strictly 
> prohibited. Neither this message nor any attachment is intended as or should 
> be construed as an offer, solicitation or recommendation to buy or sell any 
> security or other financial instrument. Neither the sender, his or her 
> employer nor any of their respective affiliates makes any warranties as to 
> the completeness or accuracy of any of the information contained herein or 
> that this message or any of its attachments is free of viruses.
> | ___
> | 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
>
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> ___
> Rcpp-devel mailing list
> Rcpp-devel@lists.r-forge.r-project.org

Re: [Rcpp-devel] NumericVector Double mismatch when indexing an array

[Douglas Bates]

On Tue, Sep 25, 2012 at 1:53 PM, Rodney Sparapani  wrote:
> On 09/25/2012 01:37 PM, Goldfeld, Keith wrote:
>
>>>  code<- 'Rcpp::RNGScope scope;
>>
>>
>> +Rcpp::NumericVector rn(5);
>>
>> +for (int i=0; i<  5; i++) {
>>
>> +rn(i) = rnorm(1,0,1);
>>
>> +}
>>
>> +return Rcpp::wrap(rn);
>>
>> +  '
>>
>>>
>>
>>>  fun<- cxxfunction(body=code, plugin="Rcpp")
>
>
> You just need an explicit type conversion like so:
>
> funk2 <- cxxfunction(body='RNGScope scope;
>   NumericVector rn(5);
>   for (int i=0; i < 5; i++) rn[i] = as(rnorm(1,0,1));
>   return wrap(rn);',
>   includes="using namespace Rcpp;\n", plugin="Rcpp")
>
> funk2()

That works but you wouldn't want to use it as a general model because
you are creating vectors of length 1 then dereferencing the first
element in each of these vectors and copying it to another vector.
Obviously with n=5 the difference in execution time will be minimal
but with n=5 million it won't be.  Using the Rcpp sugar function rnorm
will be the easiest approach unless, like me, you get a little queasy
when using Rcpp sugar functions just because it is so hard to pick
them apart and see what is actually happening.

I would probably end up going back to the R API and calling norm_rand
or Rf_rnorm.  That latter returns a double from two double arguments,
mu and sigma.
___
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

Re: [Rcpp-devel] Getting and setting array dimensions - simpler ways?

[Douglas Bates]

On Mon, Sep 10, 2012 at 8:55 AM, Søren Højsgaard  wrote:
> Dear list,
>
> I want to do some operations on an array, here exemplified by taking log (the 
> operations I have in mind are made using RcppArmadillo):
>
>> pp  <- array(c(.1, .4, .4, .1), dim=c(2,2),dimnames=list(A=c("a1","a2"), 
>> B=c("b1","b2")))
>> log(pp)
> B
> Ab1 b2
>   a1 -2.3025851 -0.9162907
>   a2 -0.9162907 -2.3025851
>
> It is essential that dim and dimnames attributes are preserved. My take on 
> this task is:
>
> src <-'
>  using namespace arma;
>  using namespace Rcpp;
>  NumericVector p(pp_);
>  // Some computations using RcppArmadillo
>  vec pp = as(pp_);
>  vec aa = log(pp);
>  // Get the result back
>  NumericVector ans(p.length());
>  ans = aa;
>  // Set attributes
>  ans.attr("dim")  = p.attr("dim");
>  ans.attr("dimnames") = p.attr("dimnames");
>  return(wrap(ans));
>  '
>> library(inline)
>> loga <- cxxfunction(signature(pp_ = ""), src, plugin = "RcppArmadillo", 
>> verbose=F)
>> loga(pp)
> B
> Ab1 b2
>   a1 -2.3025851 -0.9162907
>   a2 -0.9162907 -2.3025851
>
> This works, but if anyone can spot a more elegant (fewer lines) way of doing 
> it, I would be happy to know...

Following what Dirk said:

I think an easier way is to create an Rcpp::NumericMatrix and clone
it, then use std::transform to take the exponential of the elements.

Something like (not tested)

Rcpp::NumericMatrix p(pp_);
Rcpp::NumericMatrix q = clone(p);
std::transform(p.begin(), p.end(), q.begin(), std::ptr_func(exp));
return q;

There is a method called something like input_transform for some Rcpp
objects that may make it even easier but a quick scan of the unit
tests didn't show a use.
___
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

Re: [Rcpp-devel] Dirk's benchmarking of code for creation of random binary matrices

[Douglas Bates]

Slightly different results if I clean up the eigenFloor version and
use unif_rand() instead of runif.  The sugar version is still faster,
however.

R version 2.15.1 (2012-06-22) -- "Roasted Marshmallows"
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(inline)
> library(compiler)
> library(rbenchmark)
> 
> n <- 500
> k <- 100
> scott <- function(N, K) {
+ mm <- matrix(0, N, K)
+ apply(mm, c(1, 2), function(x) sample(c(0, 1), 1))
+ }
> scottComp <- cmpfun(scott)
> ted <- function(N, K) {
+ matrix(rbinom(N * K, 1, 0.5), ncol = K, nrow = N)
+ }
> david <- function(m, n) {
+ matrix(sample(0:1, m * n, replace = TRUE), m, n)
+ }
> luis <- function(m, n) {
+  round(matrix(runif(m * n), m, n))
+ }
> armaFloor <- cxxfunction(signature(ns="integer", ks="integer"), plugin = 
> "RcppArmadillo", body='
+int n = Rcpp::as(ns);
+int k = Rcpp::as(ks);
+return wrap(arma::floor(arma::randu(n, k) + 0.5));
+ ')
> sugarFloor <- cxxfunction(signature(ns="integer", ks="integer"), plugin = 
> "Rcpp", body='
+int n = Rcpp::as(ns);
+int k = Rcpp::as(ks);
+Rcpp::RNGScope tmp;
+Rcpp::NumericVector draws = Rcpp::floor(Rcpp::runif(n*k)+0.5);
+return Rcpp::NumericMatrix(n, k, draws.begin());
+ ')
> eigenFloor <- cxxfunction(signature(ns="integer", ks="integer"), plugin = 
> "RcppEigen", body='
+int n(Rf_asInteger(ns)), k(Rf_asInteger(ks));
+Rcpp::RNGScope tmp;
+return Rcpp::wrap(Eigen::ArrayXXd(n, k).unaryExpr(std::ptr_fun(rbinary)));
+ ', includes='
+inline double rbinary(double x) {return std::floor(unif_rand() + 0.5);}
+ ')
Loading required package: lattice

Attaching package: ‘RcppEigen’

The following object(s) are masked from ‘package:RcppArmadillo’:

fastLm, fastLmPure

> res <- benchmark(scott(n, k), scottComp(n,k),
+  ted(n, k), david(n, k), luis(n, k),
+  armaFloor(n, k), sugarFloor(n, k),
+  eigenFloor(n, k),
+  order="relative", replications=100)
> print(res[,1:4])
  test replications elapsed   relative
6  armaFloor(n, k)  100   0.145   1.00
7 sugarFloor(n, k)  100   0.161   1.110345
8 eigenFloor(n, k)  100   0.189   1.303448
4  david(n, k)  100   0.219   1.510345
3ted(n, k)  100   0.571   3.937931
5   luis(n, k)  100   0.755   5.206897
2  scottComp(n, k)  100  45.467 313.565517
1  scott(n, k)  100  45.804 315.889655
> 
> proc.time()
   user  system elapsed 
109.610   1.048 110.839 
___
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

[Rcpp-devel] Dirk's benchmarking of code for creation of random binary matrices

[Douglas Bates]

Dirk Eddelbuettel recently blogged on "Faster creation of binomial
matrices" (http://dirk.eddelbuettel.com/blog/code/snippets/) comparing
various approaches in R and using C++ callable from R with the Rcpp
package. I used cut-and-paste to create a script from the code in his
posting  but omitting the arma and sugar functions that return the
wrong answer (uniform(0,1) random variates rather than random binary
(0/1) responses with equal probability).  I added another function
using RcppEigen with disappointing results.  It came back slower than
sugarFloor and slower than David Smith's pure R versilon.  The results
are enclosed.

The armaFloor version, which is fastest, uses a different random
number generator, and thus doesn't offer the same controls as the
other versions, e.g. using set.seed in R to create a reproducible
experiment.

I decided to compare these results with Julia (and I think I got Dirk
a bit irritated by doing so).  Again, there is the caveat that these
results are using a different uniform random number generator from R's
so you should not take comparative timings too seriously.  In Julia
(www.julialang.org) all functions are generics using multiple dispatch
(like S4 generics and methods in R) and a 500 by 100 matrix of uniform
(0,1) random variates can be created as rand(500, 100).  Applying
iround to this produces 0/1 responses stored as integers.
julia> iround(rand(500, 100))
500x100 Int64 Array:
 0  0  1  0  0  1  0  1  1  0  1  1  1  :  0  0  1  1  1  0  1  0  1  0  0  1
 0  1  0  0  0  1  0  1  0  1  0  1  0 0  1  0  1  1  1  0  0  0  0  1  1
 1  1  1  1  0  1  0  1  0  0  0  1  1 1  1  1  0  0  1  1  1  1  0  0  1
 0  1  1  1  1  1  1  1  1  1  0  1  0 0  0  1  0  0  1  0  1  0  1  0  0
 1  0  1  0  0  0  0  1  1  1  1  0  1 0  1  1  1  1  0  0  0  0  0  0  0
 1  1  1  1  1  0  0  0  1  0  0  0  0  :  1  1  0  1  1  0  0  0  0  1  1  1
 0  0  0  1  1  1  0  1  0  1  0  1  0 1  0  0  0  1  1  0  1  1  0  0  1
 1  1  0  1  1  1  1  0  0  1  0  1  1 1  0  1  0  0  1  0  1  0  0  0  0
 1  1  0  1  1  1  1  0  1  0  0  0  1 1  0  1  1  0  1  1  0  1  1  0  1
 0  1  0  0  1  0  0  0  1  0  1  0  1 0  1  0  1  1  1  0  0  1  1  1  1
 :  :  :
 1  0  0  0  0  1  1  0  0  1  1  1  0 1  0  0  1  1  1  1  0  1  0  0  1
 1  1  1  1  0  0  0  0  1  1  1  0  0 0  0  1  1  0  1  1  0  0  1  1  1
 1  1  0  0  1  0  0  0  0  1  1  1  0 1  0  0  0  1  1  1  0  1  0  1  1
 1  1  1  0  1  1  1  0  0  1  0  0  0 1  1  1  0  1  1  0  1  0  1  0  1
 1  1  0  1  0  1  1  0  1  0  0  1  0  :  1  1  1  0  0  1  1  0  1  1  0  0
 1  0  1  1  0  1  0  1  0  1  1  0  0 1  0  1  1  1  0  1  1  0  0  0  1
 1  0  1  1  0  0  1  0  1  1  0  0  0 0  0  0  0  0  1  1  1  1  0  1  0
 0  0  0  0  0  1  1  1  1  0  1  1  1 1  0  1  0  0  1  1  1  0  0  1  0
 1  1  1  1  0  1  0  0  0  0  1  1  0 1  1  0  1  0  0  1  1  1  0  0  0

Timing 100 replications of this and then repeating the timing 10 times
for consistency produces
julia> [@elapsed sum([@elapsed iround(rand(500, 100)) for k in 1:100])
for l in 1:10 ]
10-element Float64 Array:
 0.0764122
 0.0762029
 0.075958
 0.0767078
 0.07639
 0.0763202
 0.078331
 0.07691
 0.085947
 0.076138

julia> [@elapsed sum([@elapsed iround(rand(500, 100)) for k in 1:100])
for l in 1:10 ]
10-element Float64 Array:
 0.078619
 0.078573
 0.0764699
 0.079267
 0.0765851
 0.0762711
 0.077647
 0.0764351
 0.0777521
 0.078192

There is a slightly faster version using a comprehension, which is a
way of creating an array that essentially folds the loops inside an
expression.

julia> [iround(rand()) for i in 1:500, j in 1:100]
500x100 Int64 Array:
 0  0  0  0  0  1  0  0  1  0  0  1  1  :  1  1  0  1  1  1  0  1  1  1  1  0
 0  1  1  1  1  1  1  1  0  1  1  1  0 0  1  0  1  1  0  0  0  0  1  1  0
 0  1  1  1  0  0  0  0  0  0  1  0  0 1  0  1  1  0  0  1  0  1  1  1  1
 1  1  0  0  1  0  0  1  1  1  1  1  1 1  0  1  1  0  1  0  0  1  1  1  0
 0  0  1  1  0  1  0  1  0  0  0  0  1 0  0  1  0  1  1  1  1  1  0  0  1
 0  0  0  0  0  0  0  1  0  1  1  0  1  :  1  0  1  1  0  1  1  1  1  0  1  1
 1  0  1  1  0  1  0  0  1  0  1  0  1 0  0  1  1  1  1  0  0  0  0  1  0
 0  1  1  1  1  1  0  0  0  0  0  0  0 1  0  0  0  1  0  1  0  1  1  0  0
 0  0  0  0  1  0  0  0  0  0  0  0  1 0  0  1  0  1  1  1  0  0  0  0  1
 1  1  1  1  1  1  1  0  1  0  1  0  0 0  1  1  1  0  1  1  0  1  1  1  1
 :  :  :
 0  0  1  0  0  1  1  1  1  0  1  0  1 1  1  0  1  1  1  0  1  1  1  0  1
 0  1  0  0  0  0  0  0  1  1  0  1  1 0  0  0  1  0  0  1  1  0  1  0  0
 1  0  0  0  1  1  1  0  0  1  0  0  1 0  0  1  1  0  0  1  0  1  0  1  0
 1  1  0  1  0  1  1  1  0  1  1  0  0 0  0  1  1  1  1  1  1  1  1  0  1
 1  0  1  0  0  0  1  1  1  0  0  1  0  :  1  1  1  0  1  0  0  1  1  0  0  0
 0  1  0  1  0  0  1  1  1  0  1  0  0 1  0  1  1  1  1  1  0  0  0  0  1
 1  1  1 

Re: [Rcpp-devel] Advice on installing CHOLMOD/Suitesparse

[Douglas Bates]

On Fri, Aug 31, 2012 at 11:35 AM, Rodney Sparapani  wrote:
> On 08/31/2012 11:20 AM, Douglas Bates wrote:
>>
>> We seem to be talking at crossed purposes.
>>
>> What I meant to say is that for some mysterious reason the standard R
>> configuration defines
>>
>> CPPFLAGS=-I/usr/local/include
>>
>> or
>>
>> CPPFLAGS=-I/opt/local/include
>>
>
> Actually, it is not mysterious.  I just checked it with R 2.15.1
> and, if you build it with a blank CPPFLAGS in your environment, then
> within $RHOME/etc/Makeconf CPPFLAGS is blank as well.  So, as I say,
> the problem was self-inflicted since our standard set up had
> CPPFLAGS=-I/opt/local/include

No.  It's in the configure script for R

## We provide these defaults so that headers and libraries in
## '/usr/local' are found (by the native tools, mostly).
if test -f "/sw/etc/fink.conf"; then
  : ${CPPFLAGS="-I/sw/include -I/usr/local/include"}
  : ${LDFLAGS="-L/sw/lib -L/usr/local/lib"}
else
  : ${CPPFLAGS="-I/usr/local/include"}
  : ${LDFLAGS="-L/usr/local/${LIBnn}"}
fi

> But, that is a habit inherited from UNIX where /opt/local (or
> /usr/local) was quite necessary for all of the GNU/FOSS stuff.
> On GNU Linux, where "everything" is already included or added via
> the package manager, this is not needed and, as you suggest,
> dangerous.
>
>> whichever is available and that doesn't make sense to me because the
>> /usr/local/  or /opt/local directories are, well, local.  So any
>> package that looks for header files in the form
>>
>> #include
>>
>> where there is a foo directory in the $PKG/inst/include (source
>> package) or $PKG/include (installed package) directory will break if
>> there happens to be a/usr/local/include/foo/  directory.
>>
>> This makes being an author of packages like RcppArmadillo or RcppEigen
>> rather difficult because the author can't just deposit the header
>> files needed by those packages linking to the Rcpp* package in
>> $PKG/inst/include. He/she must anticipate what potential users may
>> have in their /usr/local/include or /opt/local/include directories.
>> One way around this, I suppose, is for me to rename
>> RcppEigen/inst/include/Eigen and RcppEigen/inst/include/unsupported to
>> other names and go in and modify every file in those directories to
>> use nonstandard names.  Of course, this makes maintenance awkward.
>>
>> Your solution of changing the name in include/Eigen/CholmodSupport to
>>   is dangerous because it can result in version
>> skew.  The CHOLMOD code will be pulled from the Matrix package and the
>> headers for that version of CHOLMOD are those in RcppEigenCholmod.h
>> If you use your own headers you may end up with a different version of
>> many of the structs that are defined in those headers, with
>> unpredictable consequences.
>>
>
> It's not ideal.  It was just a workaround to get RcppEigen working so
> that we can start working on a new project.
>
>> To me the cleanest way out of this is to move /opt/local/include/Eigen
>> and /opt/local/include/unsupported into a subdirectory like
>> /opt/local/include/eigen3.
>>
>
> That would definitely work.  Since we always have shiny new R releases,
> I'm just upgrading to 2.15.1 with CPPFLAGS= ;o)
>
>> Did you perhaps install suitesparse and Eigen just for the purposes of
>> using RcppEigen.  If so, you don't need to have them.  Between the
>> Matrix package and Rcpp and RcppEigen all the necessary code is
>> provided in R packages.
>
>
> No, I have been using CHOLMOD for a long time.  And, recently I have
> been using Eigen as well.  But, I see your point.  Everything that
> RcppEigen needs is provided (I should have checked the source
> immediately).  Thanks
>
> --
> Rodney Sparapani, PhD  Center for Patient Care and Outcomes Research
> Sr. Biostatistician   http://www.mcw.edu/pcor
> 4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
> WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA
___
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

Re: [Rcpp-devel] Advice on installing CHOLMOD/Suitesparse

[Douglas Bates]

We seem to be talking at crossed purposes.

What I meant to say is that for some mysterious reason the standard R
configuration defines

CPPFLAGS=-I/usr/liocal/include

or

CPPFLAGS=-I/opt/local/include

whichever is available and that doesn't make sense to me because the
/usr/local/ or /opt/local directories are, well, local.  So any
package that looks for header files in the form

#include 

where there is a foo directory in the $PKG/inst/include (source
package) or $PKG/include (installed package) directory will break if
there happens to be a /usr/local/include/foo/ directory.

This makes being an author of packages like RcppArmadillo or RcppEigen
rather difficult because the author can't just deposit the header
files needed by those packages linking to the Rcpp* package in
$PKG/inst/include. He/she must anticipate what potential users may
have in their /usr/local/include or /opt/local/include directories.
One way around this, I suppose, is for me to rename
RcppEigen/inst/include/Eigen and RcppEigen/inst/include/unsupported to
other names and go in and modify every file in those directories to
use nonstandard names.  Of course, this makes maintenance awkward.

Your solution of changing the name in include/Eigen/CholmodSupport to
 is dangerous because it can result in version
skew.  The CHOLMOD code will be pulled from the Matrix package and the
headers for that version of CHOLMOD are those in RcppEigenCholmod.h
If you use your own headers you may end up with a different version of
many of the structs that are defined in those headers, with
unpredictable consequences.

To me the cleanest way out of this is to move /opt/local/include/Eigen
and /opt/local/include/unsupported into a subdirectory like
/opt/local/include/eigen3.

Did you perhaps install suitesparse and Eigen just for the purposes of
using RcppEigen.  If so, you don't need to have them.  Between the
Matrix package and Rcpp and RcppEigen all the necessary code is
provided in R packages.

On Fri, Aug 31, 2012 at 8:51 AM, Rodney Sparapani  wrote:
> On 08/30/2012 05:53 PM, Douglas Bates wrote:
>
>>
>> You're right.  I didn't see that -I/opt/local/include was in the
>> compiler call.  If you look at $RHOME/etc/Makeconf you will probably
>> see that
>>
>> CPPFLAGS=-I/opt/local/include
>>
>
> Correct!
>
>
>> (in my case /usr/local/include) for reasons that are somewhat
>> mysterious to me.  It seems like the sort of thing that a prominent
>> member of R-core would yell and scream about if he had not put that
>> there himself.
>>
>
> I can't blame anyone else.  That's part of our standard setup which
> was working beautifully until 2 days ago.  But, I will re-think it.
>
>
>> It happens that the Ubuntu package for eigen3 installs the Eigen
>> header files in /usr/include/eigen3 so I don't get bitten by that.
>>
>> It seems that you options are 1) Change R's default CPPFLAGS, which
>> seems unlikely to happen or 2) change the location of the Eigen header
>> files.
>
>
> When the next release of R appears, #1 is probably what I should look
> into since it is probably self-inflicted and, otherwise, may lead to
> other unusual bug reports in the future.  For a quick and dirty fix,
> I tweaked my original hack (after playing around with the
> cpp flag -I- and giving up in frustration).  Replace line 9 of
> PREFIX/include/Eigen/CholmodSupport with
>
> #ifdef RcppEigen__RcppEigen__h
> #include 
> #else
> #include 
> #endif
>
> Thanks for all of your help!
>
> --
> Rodney Sparapani, PhD  Center for Patient Care and Outcomes Research
> Sr. Biostatistician   http://www.mcw.edu/pcor
> 4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
> WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA
___
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

Re: [Rcpp-devel] Advice on installing CHOLMOD/Suitesparse

[Douglas Bates]

On Thu, Aug 30, 2012 at 4:29 PM, Rodney Sparapani  wrote:
> On 08/30/2012 04:14 PM, Douglas Bates wrote:
>>
>> It is quite possible that you would have gotten a buggy version if you
>> checked out a copy of the SVN archive at some random time.  Making
>> such a change is something I would tend to forget until tests started
>> failing.
>
>
> I can confirm that it is not a buggy version.  What it is happening is
> that the pre-processor is seeing my Eigen 3.1.1 installation at
> /opt/local/include/Eigen
> instead of the RcppEigen 0.3.1 at
> /opt/local/lib64/R/library/RcppEigen/include/Eigen
>
> This must be happening because of the order of the -I switches...
>
>
> # R CMD INSTALL RcppEigen_0.3.1.tar.gz
> * installing to library '/opt/local/lib64/R/library'
> * installing *source* package 'RcppEigen' ...
> ** package 'RcppEigen' successfully unpacked and MD5 sums checked
> ** libs
> g++ -I/opt/local/lib64/R/include  -I/opt/local/include
> -I"/opt/local/lib64/R/library/Rcpp/include"  -I../inst/include ...
>
> Shouldn't -I../inst/include be first?  In src/Makevars,
> there is just PKG_CXXFLAGS=-I../inst/include.  So, I guess that
> should have been prepended to the CXXFLAGS rather than
> appended.  Am I making sense?

You're right.  I didn't see that -I/opt/local/include was in the
compiler call.  If you look at $RHOME/etc/Make.conf you will probably
see that

CPPFLAGS=-I/opt/local/include

(in my case /usr/local/include) for reasons that are somewhat
mysterious to me.  It seems like the sort of thing that a prominent
member of R-core would yell and scream about if he had not put that
there himself.

It happens that the Ubuntu package for eigen3 installs the Eigen
header files in /usr/include/eigen3 so I don't get bitten by that.

It seems that you options are 1) Change R's default CPPFLAGS, which
seems unlikely to happen or 2) change the location of the Eigen header
files.


> --
> Rodney Sparapani, PhD  Center for Patient Care and Outcomes Research
> Sr. Biostatistician   http://www.mcw.edu/pcor
> 4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
> WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA
___
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

Re: [Rcpp-devel] Advice on installing CHOLMOD/Suitesparse

[Douglas Bates]

Somehow you have obtained a buggy version of RcppEigen_0.3.1.tar.gz

I just downloaded that file from a CRAN mirror and this version has
the file RcppEigen/inst/include/Eigen/CholmodSupport starting with

#ifndef EIGEN_CHOLMODSUPPORT_MODULE_H
#define EIGEN_CHOLMODSUPPORT_MODULE_H

#include "SparseCore"

#include "src/Core/util/DisableStupidWarnings.h"

extern "C" {
  #include 
}

Check your copy of that file.  It looks like it has #include
 instead of #include 

It is quite possible that you would have gotten a buggy version if you
checked out a copy of the SVN archive at some random time.  Making
such a change is something I would tend to forget until tests started
failing.

On Thu, Aug 30, 2012 at 3:44 PM, Rodney Sparapani  wrote:
> On 08/30/2012 03:24 PM, Rodney Sparapani wrote:
>>
>> # R CMD INSTALL RcppEigen_0.3.1.tar.gz
>> * installing to library '/opt/local/lib64/R/library'
>> * installing *source* package 'RcppEigen' ...
>> ** package 'RcppEigen' successfully unpacked and MD5 sums checked
>> ** libs
>> g++ -I/opt/local/lib64/R/include  -I/opt/local/include
>> -I"/opt/local/lib64/R/library/Rcpp/include"  -I../inst/include -fpic  -g
>> -c RcppEigen.cpp -o RcppEigen.o
>> In file included from ../inst/include/RcppEigenForward.h:31,
>>   from ../inst/include/RcppEigen.h:25,
>>   from RcppEigen.cpp:22:
>> /opt/local/include/Eigen/CholmodSupport:9:23: error: cholmod.h: No such
>> file or directory
>>
>> Something is pulling in the Eigen header file
>> /opt/local/include/Eigen/CholmodSupport
>
>
> Actually it is getting pulled in from RcppEigenForward.h line 31
> just as the diagnostics suggest above.  Now I am confused again ;o)
>
>
> --
> Rodney Sparapani, PhD  Center for Patient Care and Outcomes Research
> Sr. Biostatistician   http://www.mcw.edu/pcor
> 4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
> WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA
___
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

Re: [Rcpp-devel] Advice on installing CHOLMOD/Suitesparse

[Douglas Bates]

On Thu, Aug 30, 2012 at 12:41 PM, Rodney Sparapani  wrote:
> On 08/29/2012 04:42 PM, Douglas Bates wrote:
>>
>> Perhaps you have an out-of-date version of the RcppEigen package (or
>> maybe the Mac version hasn't been built or ...).  In any case
>> RcppEigenForward should include RcppEigenCholmod.h, not cholmod.h and
>> the correct version would be at
>> /opt/local/lib64/R/library/RcppEigen/include/RcppEigenCholmod.h
>
>
> Thank you Doug!  I'm not sure what I did to deserve this.  I have 3
> systems with fairly similar installations, but just the one system
> is having this problem (of course, it's the one that I really need to
> use right now).  For others that may see this issue with v. 0.3.1...
>
> So, what Doug said is correct (as usual), RcppEigen does not require
> cholmod.h  However, Eigen 3.1.1 does require cholmod.h  Therefore, in
> PREFIX/include/Eigen/CholmodSupport on line 9, there is
> #include 

But the RcppEigen package should not access Eigen 3.1.1 headers other
than the ones that it provides.  Yes, in a sane world we would be able
to install Debian packages of libraries like Eigen 3.1.1, Suitesparse,
etc. and access the code from an R package.  However doing so in a
portable way is quite nontrivial.  The Rcpp and RcppEigen packages
provide their own header files in the inst/include subdirectory of the
source package which becomes the include  subdirectory of the
installed package.  Packages that use these headers should have
"LinkingTo: Rcpp, RcppEigen" in their DESCRIPTION file so that these
header files are found.

So RcppEigen should be completely separate from any Eigen or
Suitesparse packages that you have installed on your machine.  If your
code is accessing the system Eigen or Suitesparse libraries it will
all blow up, in the way you have described.  Your code should have the
line

#include 

and nothing else related to Eigen or Suitesparse.

This is all assuming that I haven't accidentally left a line referring
to cholmod.h in the RcppEigen include files.

> But, when you install RcppEigen, RcppEigenCholmod.h is included
> which conflicts with the previous definitions/declarations that
> were brought in above.  So, as a temporary fix, I have replaced
> line 9 above with #include 
>
> This allows you to install RcppEigen.  And the RcppEigen package
> will work.  But there may be some unforeseen consequences;
> particularly, for Eigen itself rather than RcppEigen.  I
> suspect a permanent solution would involve an edit of RcppEigenCholmod.h to
> resolve the conflicts.  After looking at
> the two files, it was not immediately clear what needs to be
> changed.  But, basically, we could detect whether CHOLMOD_H
> is defined and act accordingly.  If you guys think I'm on the
> right track, then I could probably figure this out (since I am the
> one that can replicate the issue that no one else seems to see).
> If not, then please advice.  Thanks
> --
> Rodney Sparapani, PhD  Center for Patient Care and Outcomes Research
> Sr. Biostatistician   http://www.mcw.edu/pcor
> 4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
> WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA
___
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

Re: [Rcpp-devel] Advice on installing CHOLMOD/Suitesparse

[Douglas Bates]

I just checked the results page on
http://cran.r-project.org/package=RcppEigen and the package doesn't
compile solaris-sparc and solaris-ix86 with the Sun compilers (which I
believe are very old, predating modern C++) and doesn't compile on
r-oldrel-macosx-ix86 but that is due to a compiler error in yet
another very old compiler.  It seems to have compiled and passed its
checks on r-release-macosx-ix86 and I think that is now using the
clang compiler.

On Wed, Aug 29, 2012 at 4:43 PM, Douglas Bates  wrote:
> On Wed, Aug 29, 2012 at 4:42 PM, Douglas Bates  wrote:
>> On Wed, Aug 29, 2012 at 4:32 PM, Rodney Sparapani  wrote:
>>> Hi Gang:
>>>
>>> I have been using CHOLMOD/Suitesparse, but I'm having an issue
>>> coaxing RcppEigen to see it.  I have Matrix, inline, Rcpp and
>>> RcppEigen installed as well as Suitesparse.  But, when I run
>
> You should not need to have Suitesparse installed.  All the
> Suitesparse code that is used in RcppEigen is included in the Matrix
> package.
>
>>> this simple program (that works with Rcpp alone)...
>>>
>>> require(inline)
>>> require(Rcpp)
>>> require(RcppEigen)
>>>
>>> N <- 100
>>> p <- 3
>>> q <- 10
>>>
>>> data <- list(X=matrix(NA, N, p), Z=matrix(NA, N, q),
>>>  t=double(N), y=double(N), v1=integer(N), v2=integer(N),
>>>  c=integer(N))
>>>
>>> par.info <- list(beta=list(init=0., prior=list(mean=0., prec=0.001)),
>>>   mu=list(init=c(0., 0.),
>>> prior=list(mean=c(0., 0.), prec=diag(0.001, 2))),
>>>   theta=list(init=c(rep(0., ncol(data$X)),
>>>rep(1., ncol(data$Z)),
>>>rep(0., ncol(data$X))),
>>> prior=list(mean=rep(0., ncol(data$Z)+2*ncol(data$X)),
>>>prec=diag(0.001, ncol(data$Z)+2*ncol(data$X
>>>   )
>>>
>>> src <- '
>>> List list(arg1), beta=list["beta"], mu=list["mu"],
>>> theta=list["theta"];
>>>
>>> return theta["init"];
>>> '
>>>
>>> check1 <- cxxfunction(signature(arg1="list"), src, #plugin="Rcpp")
>>>   plugin="RcppEigen")
>>>
>>> I get...
>>> In file included from
>>> /opt/local/lib64/R/library/RcppEigen/include/RcppEigenForward.h:31,
>>>  from
>>> /opt/local/lib64/R/library/RcppEigen/include/RcppEigen.h:25,
>>>  from file7289313ed735.cpp:3:
>>> /opt/local/include/Eigen/CholmodSupport:9:23: error: cholmod.h: No such 
>>> file or directory
>>
>> Perhaps you have an out-of-date version of the RcppEigen package (or
>> maybe the Mac version hasn't been built or ...).  In any case
>> RcppEigenForward should include RcppEigenCholmod.h, not cholmod.h and
>> the correct version would be at
>> /opt/local/lib64/R/library/RcppEigen/include/RcppEigenCholmod.h
>>
>>> make: *** [file7289313ed735.o] Error 1
>>>
>>> ERROR(s) during compilation: source code errors or compiler configuration
>>> errors!
>>>
>>> Program source:
>>>   1:
>>>   2: // includes from the plugin
>>>   3: #include 
>>>   4: #include 
>>>   5:
>>>   6:
>>>   7: #ifndef BEGIN_RCPP
>>>   8: #define BEGIN_RCPP
>>>   9: #endif
>>>  10:
>>>  11: #ifndef END_RCPP
>>>  12: #define END_RCPP
>>>  13: #endif
>>>  14:
>>>  15: using namespace Rcpp;
>>>  16:
>>>  17:
>>>  18: // user includes
>>>  19:
>>>  20:
>>>  21: // declarations
>>>  22: extern "C" {
>>>  23: SEXP file7289313ed735( SEXP arg1) ;
>>>  24: }
>>>  25:
>>>  26: // definition
>>>  27:
>>>  28: SEXP file7289313ed735( SEXP arg1 ){
>>>  29: BEGIN_RCPP
>>>  30:
>>>  31: List list(arg1), beta=list["beta"], mu=list["mu"],
>>>  32: theta=list["theta"];
>>>  33:
>>>  34: return theta["init"];
>>>  35:
>>>  36: END_RCPP
>>>  37: }
>>>  38:
>>>  39:
>>> Error in compileCode(f, code, language = language, verbose = verbose) :
>>>   Compilation ERROR, function(s)/method(s) not created! In file included
>>> from /opt/local/lib

Re: [Rcpp-devel] Advice on installing CHOLMOD/Suitesparse

[Douglas Bates]

On Wed, Aug 29, 2012 at 4:42 PM, Douglas Bates  wrote:
> On Wed, Aug 29, 2012 at 4:32 PM, Rodney Sparapani  wrote:
>> Hi Gang:
>>
>> I have been using CHOLMOD/Suitesparse, but I'm having an issue
>> coaxing RcppEigen to see it.  I have Matrix, inline, Rcpp and
>> RcppEigen installed as well as Suitesparse.  But, when I run

You should not need to have Suitesparse installed.  All the
Suitesparse code that is used in RcppEigen is included in the Matrix
package.

>> this simple program (that works with Rcpp alone)...
>>
>> require(inline)
>> require(Rcpp)
>> require(RcppEigen)
>>
>> N <- 100
>> p <- 3
>> q <- 10
>>
>> data <- list(X=matrix(NA, N, p), Z=matrix(NA, N, q),
>>  t=double(N), y=double(N), v1=integer(N), v2=integer(N),
>>  c=integer(N))
>>
>> par.info <- list(beta=list(init=0., prior=list(mean=0., prec=0.001)),
>>   mu=list(init=c(0., 0.),
>> prior=list(mean=c(0., 0.), prec=diag(0.001, 2))),
>>   theta=list(init=c(rep(0., ncol(data$X)),
>>rep(1., ncol(data$Z)),
>>rep(0., ncol(data$X))),
>> prior=list(mean=rep(0., ncol(data$Z)+2*ncol(data$X)),
>>prec=diag(0.001, ncol(data$Z)+2*ncol(data$X
>>   )
>>
>> src <- '
>> List list(arg1), beta=list["beta"], mu=list["mu"],
>> theta=list["theta"];
>>
>> return theta["init"];
>> '
>>
>> check1 <- cxxfunction(signature(arg1="list"), src, #plugin="Rcpp")
>>   plugin="RcppEigen")
>>
>> I get...
>> In file included from
>> /opt/local/lib64/R/library/RcppEigen/include/RcppEigenForward.h:31,
>>  from
>> /opt/local/lib64/R/library/RcppEigen/include/RcppEigen.h:25,
>>  from file7289313ed735.cpp:3:
>> /opt/local/include/Eigen/CholmodSupport:9:23: error: cholmod.h: No such file 
>> or directory
>
> Perhaps you have an out-of-date version of the RcppEigen package (or
> maybe the Mac version hasn't been built or ...).  In any case
> RcppEigenForward should include RcppEigenCholmod.h, not cholmod.h and
> the correct version would be at
> /opt/local/lib64/R/library/RcppEigen/include/RcppEigenCholmod.h
>
>> make: *** [file7289313ed735.o] Error 1
>>
>> ERROR(s) during compilation: source code errors or compiler configuration
>> errors!
>>
>> Program source:
>>   1:
>>   2: // includes from the plugin
>>   3: #include 
>>   4: #include 
>>   5:
>>   6:
>>   7: #ifndef BEGIN_RCPP
>>   8: #define BEGIN_RCPP
>>   9: #endif
>>  10:
>>  11: #ifndef END_RCPP
>>  12: #define END_RCPP
>>  13: #endif
>>  14:
>>  15: using namespace Rcpp;
>>  16:
>>  17:
>>  18: // user includes
>>  19:
>>  20:
>>  21: // declarations
>>  22: extern "C" {
>>  23: SEXP file7289313ed735( SEXP arg1) ;
>>  24: }
>>  25:
>>  26: // definition
>>  27:
>>  28: SEXP file7289313ed735( SEXP arg1 ){
>>  29: BEGIN_RCPP
>>  30:
>>  31: List list(arg1), beta=list["beta"], mu=list["mu"],
>>  32: theta=list["theta"];
>>  33:
>>  34: return theta["init"];
>>  35:
>>  36: END_RCPP
>>  37: }
>>  38:
>>  39:
>> Error in compileCode(f, code, language = language, verbose = verbose) :
>>   Compilation ERROR, function(s)/method(s) not created! In file included
>> from /opt/local/lib64/R/library/RcppEigen/include/RcppEigenForward.h:31,
>>  from
>> /opt/local/lib64/R/library/RcppEigen/include/RcppEigen.h:25,
>>  from file7289313ed735.cpp:3:
>> /opt/local/include/Eigen/CholmodSupport:9:23: error: cholmod.h: No such file
>> or directory
>> make: *** [file7289313ed735.o] Error 1
>> In addition: Warning message:
>> running command '/opt/local/lib64/R/bin/R CMD SHLIB file7289313ed735.cpp 2>
>> file7289313ed735.cpp.err.txt' had status 1
>>
>> So, I'm wondering...  Where does RcppEigen expect to see cholmod.h?
>> Perhaps, I should re-install RcppEigen from source rather than from
>> CRAN.  Or possibly install Suitesparse from source in either
>> /usr/local, or preferably, /opt/local?  What works for other?  Thanks
>>
>> Compiler Settings:
>> CC="gcc"
>> CFLAGS="-g&q

Re: [Rcpp-devel] Advice on installing CHOLMOD/Suitesparse

[Douglas Bates]

On Wed, Aug 29, 2012 at 4:32 PM, Rodney Sparapani  wrote:
> Hi Gang:
>
> I have been using CHOLMOD/Suitesparse, but I'm having an issue
> coaxing RcppEigen to see it.  I have Matrix, inline, Rcpp and
> RcppEigen installed as well as Suitesparse.  But, when I run
> this simple program (that works with Rcpp alone)...
>
> require(inline)
> require(Rcpp)
> require(RcppEigen)
>
> N <- 100
> p <- 3
> q <- 10
>
> data <- list(X=matrix(NA, N, p), Z=matrix(NA, N, q),
>  t=double(N), y=double(N), v1=integer(N), v2=integer(N),
>  c=integer(N))
>
> par.info <- list(beta=list(init=0., prior=list(mean=0., prec=0.001)),
>   mu=list(init=c(0., 0.),
> prior=list(mean=c(0., 0.), prec=diag(0.001, 2))),
>   theta=list(init=c(rep(0., ncol(data$X)),
>rep(1., ncol(data$Z)),
>rep(0., ncol(data$X))),
> prior=list(mean=rep(0., ncol(data$Z)+2*ncol(data$X)),
>prec=diag(0.001, ncol(data$Z)+2*ncol(data$X
>   )
>
> src <- '
> List list(arg1), beta=list["beta"], mu=list["mu"],
> theta=list["theta"];
>
> return theta["init"];
> '
>
> check1 <- cxxfunction(signature(arg1="list"), src, #plugin="Rcpp")
>   plugin="RcppEigen")
>
> I get...
> In file included from
> /opt/local/lib64/R/library/RcppEigen/include/RcppEigenForward.h:31,
>  from
> /opt/local/lib64/R/library/RcppEigen/include/RcppEigen.h:25,
>  from file7289313ed735.cpp:3:
> /opt/local/include/Eigen/CholmodSupport:9:23: error: cholmod.h: No such file 
> or directory

Perhaps you have an out-of-date version of the RcppEigen package (or
maybe the Mac version hasn't been built or ...).  In any case
RcppEigenForward should include RcppEigenCholmod.h, not cholmod.h and
the correct version would be at
/opt/local/lib64/R/library/RcppEigen/include/RcppEigenCholmod.h

> make: *** [file7289313ed735.o] Error 1
>
> ERROR(s) during compilation: source code errors or compiler configuration
> errors!
>
> Program source:
>   1:
>   2: // includes from the plugin
>   3: #include 
>   4: #include 
>   5:
>   6:
>   7: #ifndef BEGIN_RCPP
>   8: #define BEGIN_RCPP
>   9: #endif
>  10:
>  11: #ifndef END_RCPP
>  12: #define END_RCPP
>  13: #endif
>  14:
>  15: using namespace Rcpp;
>  16:
>  17:
>  18: // user includes
>  19:
>  20:
>  21: // declarations
>  22: extern "C" {
>  23: SEXP file7289313ed735( SEXP arg1) ;
>  24: }
>  25:
>  26: // definition
>  27:
>  28: SEXP file7289313ed735( SEXP arg1 ){
>  29: BEGIN_RCPP
>  30:
>  31: List list(arg1), beta=list["beta"], mu=list["mu"],
>  32: theta=list["theta"];
>  33:
>  34: return theta["init"];
>  35:
>  36: END_RCPP
>  37: }
>  38:
>  39:
> Error in compileCode(f, code, language = language, verbose = verbose) :
>   Compilation ERROR, function(s)/method(s) not created! In file included
> from /opt/local/lib64/R/library/RcppEigen/include/RcppEigenForward.h:31,
>  from
> /opt/local/lib64/R/library/RcppEigen/include/RcppEigen.h:25,
>  from file7289313ed735.cpp:3:
> /opt/local/include/Eigen/CholmodSupport:9:23: error: cholmod.h: No such file
> or directory
> make: *** [file7289313ed735.o] Error 1
> In addition: Warning message:
> running command '/opt/local/lib64/R/bin/R CMD SHLIB file7289313ed735.cpp 2>
> file7289313ed735.cpp.err.txt' had status 1
>
> So, I'm wondering...  Where does RcppEigen expect to see cholmod.h?
> Perhaps, I should re-install RcppEigen from source rather than from
> CRAN.  Or possibly install Suitesparse from source in either
> /usr/local, or preferably, /opt/local?  What works for other?  Thanks
>
> Compiler Settings:
> CC="gcc"
> CFLAGS="-g"
> CXX="g++"
> CXXFLAGS="-g"
> CPP=""
> CPPFLAGS="-I/opt/local/include"
> LDFLAGS="-L/opt/local/lib64 -Wl,-rpath -Wl,/opt/local/lib64"
> LIBS=""
> PKG_CONFIG_PATH="/opt/local/lib64/pkgconfig"
> F77="gfortran"
> FFLAGS="-g -ff2c"
> FC="gfortran"
> FCFLAGS=""
>
> sessionInfo()
> R version 2.14.2 (2012-02-29)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
>  [1] LC_CTYPE=en_US.ISO8859-1   LC_NUMERIC=C
>  [3] LC_TIME=en_US.ISO8859-1LC_COLLATE=en_US.ISO8859-1
>  [5] LC_MONETARY=en_US.ISO8859-1LC_MESSAGES=en_US.ISO8859-1
>  [7] LC_PAPER=C LC_NAME=C
>  [9] LC_ADDRESS=C   LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.ISO8859-1 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods   base
>
> other attached packages:
> [1] RcppEigen_0.2.0 Rcpp_0.9.10 inline_0.3.8
>
> loaded via a namespace (and not attached):
> [1] compiler_2.14.2 grid_2.14.2 lattice_0.20-6  Matrix_1.0-5
> [5] tools_2.14.2
>
> yum install suitesparse
> Loaded plugins: auto-update-debuginfo, refresh-packagekit, rhnplugin
> Setting up Install Process
> Package suitesparse-3.4.0-2.el6.x86_64 already installed and latest version
>

Re: [Rcpp-devel] Why inline function is much faster than .Call?

[Douglas Bates]

On Mon, Aug 27, 2012 at 3:17 PM, Peng Yu  wrote:
> Hi Dirk,
>
>> inlines uses .Call, so there is a slight logical problem here...  And yes,
>
> What is the logical problem?

Dirk is trying to point out that the two calls are using the same
mechanism.  All that the inline package does is wrap the process of
compiling, linking and dynamically loading the C++ function.
Executing a function compiled this way is exactly the same as
executing a function compiled separately.

Benchmarking a function like this only once is liable to be unreliable
because the current state of the garbage collection and heap
allocation can change running times.  It is better to run the function
a couple of times before benchmarking.

>> functions have overheads.
>
> If that is the case, .Call("test", xx) should be faster than test(xx).
> But it is not (see below). This looks wired to me.
>
>>
>> | library(microbenchmark)
>> | microbenchmark(test(xx))
>> | microbenchmark(test_inline(xx))
>> | microbenchmark(.Call('test', xx))
>>
>> The normal idiom is a _single_ call so that you can compare:
>>
>>   microbenchmark(test(xx), test_inline(xx), .Call('test', xx))
>
> Here is the output.
>
>> microbenchmark(test(xx), test_inline(xx), .Call('test', xx))
> Unit: microseconds
>expr  min   lq
>  median   uq
> max
> 1 .Call("test", xx) 94.8319993633537 96.0904991587174
> 97.5570002160050 99.7624988631316
> 128.830998885869
> 2   test_inline(xx)  4.440390799  5.1014999701572
> 6.393682121  6.615213163
> 30.338999863576
> 3  test(xx) 69.6830006934897 70.3875002842171
> 71.5935005911716 73.1535008185452
> 198.308999749889
>
> --
> Regards,
> Peng
> ___
> 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
___
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

Re: [Rcpp-devel] Convertion from SEXP to std::vector

[Douglas Bates]

On Wed, Aug 22, 2012 at 9:55 AM, Peng Yu  wrote:
> Hi,
>
> Rcpp::IntegerVector can not be converted to std::vector. Is
> there a class that can help convert SEXP to std::vector?
>
> The current walkaround can be something like the following. I'm
> wondering if there is a better solution?
>
> Rcpp::IntegerVector vec(vx);

> std::vector v;

You should at least specify the size there so that push_back doesn't
need to allocate/reallocate and possibly overallocate.

std::vector v(vec.size());

You may be able to allocate from iterators as in

std::vector v(vec.begin(), vec.end());

but I'm not sure if that will do the automatic conversion from signed
int to unsigned.  Look up the different forms of constructors for
std::vector and the STL algorithms such as std::transform.

> for (int i=0; i v.push_back(vec[i]);
> }
>
> --
> Regards,
> Peng
> ___
> 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
___
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

Re: [Rcpp-devel] How to manage testing code for Rcpp based project?

[Douglas Bates]

On Wed, Aug 22, 2012 at 12:16 AM, Davor Cubranic  wrote:
> On 12-08-21 09:54 PM, Darren Cook wrote:

 RUnit works, but the setup is not the easiest. If someone has ideas for
 better solutions, by all means implement those.  Unit testing has helped
 us
 many times, and I consider it to be a very important part of package
 development.
>>>
>>>
>>> How (or in what aspect) RUnit is better than devtools?
>>
>>
>> What is "devtools"? Is it a package for testing R software? (it sounds
>> more generic)
>
>
> See http://cran.r-project.org/web/packages/devtools/ and
> https://github.com/hadley/devtools
>
> It aims to make various phases of package development cycle easier, at the
> expense of requiring following certain conventions about naming and layout
> of source files and directories.
>
> I think what Peng meant as an alternative to Runit is not the devtools
> itself, but testthat. See http://cran.r-project.org/web/packages/testthat/
> and https://github.com/hadley/test_that
>
> Philosophically, testthat draws on Ruby test frameworks like rspec and other
> BDD tools. In practical terms, it means that the tests look like this:
>
> context("Basic tests")
>
> test_that("logical tests act as expected", {
>   expect_that(TRUE, is_true())
>   expect_that(FALSE, is_false())
> })
>
> test_that("equality holds", {
>   expect_that(5, equals(5))
>   expect_that(10, is_identical_to(10))
> })
>
> I haven't gotten very deep into using it, but in practical terms I find
> there is almost 1:1 correspondence in testing functions ("checkEquals" ->
> "expect_that(..., equals(...))", with testthat being a little more readable
> but requiring me to turn off ESS's smart-underscore. :-)

In the past it was difficult to have tests run through testthat use
the inline package because the two packages disagreed on the current
working directory.  Not sure if that has changed.
___
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

Re: [Rcpp-devel] RNGScope

[Douglas Bates]

I think you want to move the RNGScope declaration before the calls to
rcauchy and runif.

On Tue, Aug 14, 2012 at 9:17 AM, Rodney Sparapani  wrote:
> This is my first post so please be gentle ;o)  I have been through
> all of the docs that I can find (but I haven't looked at the source
> code directly).  Here is the behavior that I see:
>
>> require(inline)
> Loading required package: inline
>> source("helper.R")
>> source("corrMH.R")
>>
>> set.seed(66)
>>
>> for(i in 1:3) print(corrMH(0.2, 0.5, 100, 0))
> [1] -0.07071726
> [1] -8.085111
> [1] -4.004675
>>
>> corrMH.Rcpp <- cxxfunction(
> +  signature(a="numeric", b="numeric", N="numeric", x="numeric"),
> +  body=cat.file("corrMH.Rcpp.cpp"), plugin="Rcpp")
>>
>> set.seed(66)
>>
>> for(i in 1:3) print(corrMH.Rcpp(0.2, 0.5, 100, 0))
> [1] -0.07071726
> [1] -0.07071726
> [1] -0.07071726
>
> Notice that the Rcpp version of the function returns the same value
> each time (only 3 shown, but this is true no matter how many times
> that you do it).  I guess this is because we follow the RNGScope state
> of the RNG, but we don't update it.  However, updating is really needed
> in this case.  I hope I am making this clear.  Please let me know if
> you have any ideas.  Thanks
>
> Attachments seem to be allowed so I'm including my code...
>
> Red Hat Enterprise Linux Server release 6.2 (Santiago)
> g++ (GCC) 4.4.6 20110731 (Red Hat 4.4.6-3)
> R version 2.14.2 (2012-02-29)
>
>> installed.packages()["Rcpp", ]
>  Package  LibPath
>   "Rcpp" "/opt/local/lib64/R/library"
>  Version Priority
> "0.9.10"   NA
>  Depends  Imports
> "R (>= 2.12.0), methods"   NA
>LinkingTo Suggests
>   NA  "RUnit, inline, rbenchmark"
> Enhances  OS_type
>   NA   NA
>  LicenseBuilt
> "GPL (>= 2)" "2.14.2"
> --
> Rodney Sparapani, PhD  Center for Patient Care and Outcomes Research
> Sr. Biostatistician   http://www.mcw.edu/pcor
> 4 wheels good, 2 wheels better!   Medical College of Wisconsin (MCW)
> WWLD?:  What Would Lombardi Do?   Milwaukee, WI, USA
>
> ___
> 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
___
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

Re: [Rcpp-devel] clarification of interaction between RInside and RcppEigen

[Douglas Bates]

Probably the best place to start with RcppEigen is the vignette that
Dirk and I wrote for the RcppEigen package.  I can email you a copy of
the PDF file off-list if you wish.

With RcppEigen you can skip the creation of an Rcpp::NumericMatrix
object if you wish and go directly to a mapped Eigen Matrix.  There is
a specialization of the Rcpp::as templated function that allows for
creating an Eigen::MatrixXd or an Eigen::Map object
directly from an SEXP, which is the R's internal representation of
data.

On Tue, Aug 7, 2012 at 12:55 PM, Stephen J. Barr  wrote:
> Greetings,
>
> I am new to Rcpp, RInside and the entire family of related packages. I am
> experienced with R and also C++ using Eigen but I have never combined the
> two before. I would like to create a matrix using R and then play with it
> using Eigen. From what I have read so far, it seems like I can use RInside
> to create an Rcpp::NumericMatrix.
>
> My question is this: to transform a matrix from a Rcpp::NumericMatrix to an
> Eigen MatrixXd, do I use RcppEigen? Or is RcppEigen only for calling
> C++/Eigen from R and not the other way around?
>
> Thanks,
> Stephen
>
> ___
> 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
___
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

Re: [Rcpp-devel] How to manage testing code for Rcpp based project?

[Douglas Bates]

On Sat, Aug 4, 2012 at 11:48 AM, Peng Yu  wrote:
>> Seriously:  please think through what you are suggesting here. We are talking
>> about __compiled language__ and outside of a time machine I see very few
>> tools that would permit you to __test before installation__.
>
> You clearly exaggerated my point.
>
> For example, if you have 100 Rcpp based files in your package. But you
> just changed 1 file. To test this 1 file, you have to compile the
> whole package which take the time to compile all the 100 Rcpp files. A
> more efficient way is just to compile the file that was changed and
> possibly other files that depends on it.

But this is a matter of how you do the installation.  If you install
directly from the development directory and don't clean the object
files

R CMD INSTALL ~/Rpkgs/myPackage

then only the changed source files will be recompiled.  I prefer not
to do this because I don't want object files left lying around but
then I don't have 100 source files in any of my packages.
___
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

Re: [Rcpp-devel] How to manage testing code for Rcpp based project?

[Douglas Bates]

On Sat, Aug 4, 2012 at 11:22 AM, Peng Yu  wrote:
> Hi,
>
> I have trouble in finding a good strategy for organizing test cases to
> have complete coverage for Rcpp based code.
>
> The example in RcppArmadillo_0.3.2.4/inst/unitTests has some problem
> in the sense that, each file is too large --- if there is a little
> change the whole test file should be run, which takes more time than
> necessary. I think that having a more granular test file organization
> can be beneficial. Also, the test file is not in close proximity to
> the function being tested. This cause inefficient when navigating
> between the testing code and the function to be tested.
>
> There are many ways in solving these problems in other languages. For
> example, primitive test code can be embedded in the docstring in
> python. I'm wondering if anybody knows the best software engineering
> practice for working with Rcpp code.
>
> BTW, the test code in RcppArmadillo_0.3.2.4/inst/unitTests can be run
> only after RcppArmadillo is installed. For a developer, he/she does
> not always want to install the package before he/she can run the test
> case. "devtools" solves this problem for pure R package. I'm wondering
> if there is any equivalent thing for Rcpp based packages.

Because Rcpp-based packages are testing interfaces to compiled code,
they must either be installed or use the inline package to compile the
code on the fly.  Either approach is going to be slow.

As far as I know, it is still not possible to use the testthat
package, which is what devtools expects for a test harness, with
inline.  There are (or were) conflicting ideas of what the current
working directory is.
> Thanks!
>
> --
> Regards,
> Peng
> ___
> 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
___
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

Re: [Rcpp-devel] bounds checking disabled for RcppEigen?

[Douglas Bates]

I don't really think that bounds checking is slowing down this
calculation.  You can make the whole thing run faster by the simple
expedient of reordering the loops.  For me, this version

STLvectortest<-'
const NumericMatrix X = Rcpp::as(mat);
int p = X.nrow();
std::vector V(p*p);

// Copy data into STL vector and Eigen dynamic matrix
for (int i=0; ihttps://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel

Re: [Rcpp-devel] bounds checking disabled for RcppEigen?

[Douglas Bates]

On Wed, Aug 1, 2012 at 12:46 PM, Andrew Slaughter
 wrote:
> First of all, great job with Rcpp - it's an awesome tool, and finally gave
> me a reason to go back to my C++ books and get beyond "Hello World"!
>
> I did have a question about RcppEigen, though: specifically, when I use
> RcppEigen through inline, is bounds checking disabled? If not, how could I
> go about disabling it? I did try adding the eigen_no_debug flag to my
> includes, but it wouldn't compile. (I know the preferred way might be to
> build a package, but inline is so much faster for quickly testing changes to
> my code, I hate to give it up right now.)

Actually I think the debug checking is, by default, turned off in R
package compilation with the -DNDEBUG flag unless you work out a way
to override it.  I'm not sure if that applies in code compiled through
the inline package.

One of the things I would definitely do in your RcppEigen is to use a
Map not a MatrixXd for the input value, as this avoids a
copy operation.  Another thing to check is whether the colwise or
rowwise methods can be used to advantage.  I'll take a look and see if
I can come up with a faster formulation in RcppEigen.

> As a separate but related question: is there agreement about good dynamic
> (not fixed-size) C++ analogues to R's array? Especially anaologues that are
> Rcpp-friendly. :) I really want to be able to quickly loop over and access
> elements using something like "array(i,j,k,l,m)".
>
> BTW, sorry if these are simple questions, but I'm learning C/C++ for this
> project as I go along.
>
> Thanks,
> Andy
>
>
> Background for the curious:
>
> The vast majority of my code is looping over elements of multidimensional
> numeric arrays (not fixed at compile time) and performing simple
> multiplication, addition, and subtraction on those elements. Basically, the
> time required to access elements of a numeric array is pretty important. A
> much much smaller set of my code requires some matrix algebra and linear
> algebra routines (matrix multiplication, cholesky decomp, etc.) on small
> vectors and matrices. Right now I use BLAS/LAPACK, but it's harder for me to
> read and experiment with, so I thought I'd give Eigen a try.
>
> Since RcppArmadillo was much slower than using STL vectors last time I tried
> it, I decided to benchmark Eigen before rewriting any code (see code below
> comparing STL vector, NumericMatrix, and MatrixXd). It seems a good bit
> slower, which I suspect might be due to bounds checking when I access
> elements of the matrix.
>
> Ultimately, I guess I could mix and match if I have to (use STL vectors for
> the simple operations, and other classes for the matrix algebra stuff), but
> for some reason that "feels" wrong. but I'm no C++ person, so I'm
> probably wrong.
>
>
> # Here's a copy of the R syntax for the benchmark I mention above:
>
> library(RcppEigen)
> library(inline)
> library(rbenchmark)
>
> mat1<-matrix(rbinom(200*200, 1, 0.5), ncol=200)
>
> header<-'
> using namespace std;
> using namespace Eigen;
> '
> RcppMatrixtest<-'
> const NumericMatrix X = Rcpp::as(mat);
> int p = X.nrow();
> std::vector V(p*p);
> Eigen::MatrixXd EM(p,p);
>
> // Copy data into STL vector and Eigen dynamic matrix
> for (int i=0; i for (int j=0; j EM(i,j) = V[i + j*p] = X(i,j);
> }
> }
>
> double tran = 0.0;
>
> for (int i=0; i for (int j=0; j for (int k=0; k if ((j != i) && (k != i) && (j != k)) {
> tran += X(i,j)*X(j,k)*X(i,k);
> }
> }
> }
> }
>
> return Rcpp::wrap(tran);
> '
> rcpp<-cxxfunction(
> signature(
> mat="numeric"
> ),
> body=RcppMatrixtest,
> includes=header,
> plugin="RcppEigen",
> verbose=T
> )
>
> STLvectortest<-'
> const NumericMatrix X = Rcpp::as(mat);
> int p = X.nrow();
> std::vector V(p*p);
> Eigen::MatrixXd EM(p,p);
>
> // Copy data into STL vector and Eigen dynamic matrix
> for (int i=0; i for (int j=0; j EM(i,j) = V[i + j*p] = X(i,j);
> }
> }
>
> double tran = 0.0;
>
> for (int i=0; i for (int j=0; j for (int k=0; k if ((j != i) && (k != i) && (j != k)) {
> tran += V[i + j*p]*V[j + k*p]*V[i + k*p];
> }
> }
> }
> }
>
> return Rcpp::wrap(tran);
> '
> stl<-cxxfunction(
> signature(
> mat="numeric"
> ),
> body=STLvectortest,
> includes=header,
> plugin="RcppEigen",
> verbose=T
> )
>
> Reigentest<-'
> const NumericMatrix X = Rcpp::as(mat);
> int p = X.nrow();
> std::vector V(p*p);
> Eigen::MatrixXd EM(p,p);
>
> // Copy data into STL vector and Eigen dynamic matrix
> for (int i=0; i for (int j=0; j EM(i,j) = V[i + j*p] = X(i,j);
> }
> }
>
> double tran = 0.0;
>
> for (int i=0; i for (int j=0; j for (int k=0; k if ((j != i) && (k != i) && (j != k)) {
> tran += EM(i,j)*EM(j,k)*EM(i,k);
> }
> }
> }
> }
>
> return Rcpp::wrap(tran);
> '
>

Re: [Rcpp-devel] dyn.load error - symbol not found - expected in: flat namespace

[Douglas Bates]

One way forward is to use a program like c++filt to demangle the name
that is not found.  It is

re2::RE2::Arg::parse_float(char const*, int, void*)

On Wed, Jul 18, 2012 at 11:18 AM, Lescai, Francesco  wrote:
> I'm still having this problem.
> I progressively increased the complexity of my code, because my ultimate
> goal is to integrate a series of C++ classes to perform some functions
> external to R.
>
> My test examples worked well, and I just noticed in order to generate the
> shared object I have to compile together all the classes connected to each
> others, like this.
> PKG_CPPFLAGS=`Rscript -e 'Rcpp:::CxxFlags()'` PKG_LIBS=`Rscript -e
> 'Rcpp:::LdFlags()'` R CMD SHLIB ClassFive.cpp ClassOne.cpp ClassTwo.cpp
> ClassThree.cpp ClassFour.cpp
>
> When I finally compiled my real code, I didn't get any compilation error,
> but I still got from within R the error message:
> Error in dyn.load("ClassFive.so") :
>   unable to load shared object [...]
> Symbol not found: __ZN3re23RE23Arg11parse_floatEPKciPv
>
> I checked for the issues that generated the same problem before, i.e.
> constructor and deconstructors, and they are all ok.
> I couldn't find any explicit call to "parse_float" which seems to be
> reported in the error message.
> To start, I linked only one class to R and one method only of the very same
> class.
>
> I understand it might be quite difficult to identify the specific problem
> without going through lots of code (I have now 5 different classes, each
> with several methods).
> Is there a number of possible mistakes leading to such kind of messages,
> that could help me narrow down the cause?
>
> Thanks very much for any help,
> Francesco
>
>
> 
>> sessionInfo()
> R version 2.15.1 (2012-06-22)
> Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit)
>
> locale:
> [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
>
> attached base packages:
> [1] stats graphics  grDevices utils datasets  methods   base
>
> other attached packages:
> [1] inline_0.3.8 Rcpp_0.9.13
>
> loaded via a namespace (and not attached):
> [1] tcltk_2.15.1 tools_2.15.1
>
>
>
>
> On 6 Jul 2012, at 15:50, Lescai, Francesco wrote:
>
> Argh.. Apologies guys.. Found the error myself.
> Constructor and deconstructor must be specified with {} even if no code is
> foreseen for them.
>
> A novice error, hope at least highlighting it could be useful to other
> newbies like me :-)
>
> cheers,
> Francesco
>
>
> On 6 Jul 2012, at 15:37, Lescai, Francesco wrote:
>
> Hi there,
>
> I've seen other posts similar to this one, but I'm a complete novice in the
> use of Rcpp and couldn't really figure out how to solve the issue.
>
>
> I'm learning how to use Rcpp before connecting R to some C++ classes I'm
> developing.
>
> I started with a simple home made example, but in both cases compiling .cpp
> and header files or compiling inline code, I get the same outcome error
> "unable to load shared object" and then "Symbol not found" with some
> characters before and after my class name.
>
>
> I've seen Mac OS might have some issues, therefore I tested it also on an
> Ubuntu virtual machine, but the result is the same error message.
>
> Also, I'm using an R-devel version here but I'm having the same problem with
> R 14 as well.
>
> I'll copy below all the relevant information (bit lengthy, I'm sorry).
>
>
> I'd really appreciate some help here to point me in the right direction.
>
> thanks very much,
>
> Francesco
>
>
>
> --case 1 - external files ---
>
>
> PKG_CPPFLAGS=`Rscript -e 'Rcpp:::CxxFlags()'` PKG_LIBS=`Rscript -e
> 'Rcpp:::LdFlags()'` R CMD SHLIB example.cpp
>
> g++ -arch x86_64 -I/Library/Frameworks/R.framework/Resources/include
> -I/Library/Frameworks/R.framework/Resources/include/x86_64 -DNDEBUG
> -I/Library/Frameworks/R.framework/Versions/2.16/Resources/library/Rcpp/include
> -I/usr/local/include-fPIC  -g -O2  -c example.cpp -o example.o
>
> g++ -arch x86_64 -dynamiclib -Wl,-headerpad_max_install_names -undefined
> dynamic_lookup -single_module -multiply_defined suppress -L/usr/local/lib -o
> example.so example.o
> /Library/Frameworks/R.framework/Versions/2.16/Resources/library/Rcpp/lib/x86_64/libRcpp.a
> -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework
> -Wl,CoreFoundation
>
>
> library(Rcpp)
>
> library(inline)
>
> dyn.load("example.so")
>
> Error in dyn.load("example.so") :
>
> unable to load shared object
> '/Users/rehbfle/Documents/CPPexercise/RCPP/example.so':
>
> dlopen(/Users/rehbfle/Documents/CPPexercise/RCPP/example.so, 6): Symbol not
> found: __ZN7exampleD1Ev
>
> Referenced from: /Users/rehbfle/Documents/CPPexercise/RCPP/example.so
>
> Expected in: flat namespace
>
> in /Users/rehbfle/Documents/CPPexercise/RCPP/example.so
>
>
>
> -case 2 - inline code-
>
>
> fx<-cxxfunction(signature(), plugin="Rcpp", include=inc)
>
> Error in dyn.load(libLFile) :
>
> unable to load shared object
> '/var/folders/qj/p9

Re: [Rcpp-devel] Very Large Matrices in RcppArmadillo

[Douglas Bates]

On Tue, Jul 17, 2012 at 10:44 AM, French, Joshua
 wrote:
> Thank you all for the responses.
>
> Christian, I didn't know about the copy_aux_mem option.  I will have to
> take a look at that.
>
> Dirk, thanks for looking into the 64-bit matrix indices.
>
> Doug, the place in my code where I get the error is when I multiply
> matrices.  I might have matrices X and Y, where X is 30x500 and Y is
> 500x30 and I want Z = X * Y.

You will need to come up with another algorithm.  The size of Z is
3,000,000 by 3,000,000 and storing this as a dense matrix will require
about 67 terabytes of memory

> (300 * 300 * 8) / 2^30
[1] 67055.23

I don't think even Google has that much memory available.  You are
going to have to work with a decomposition.  For example a QR
decomposition can represent the matrix X as a (virtual) orthogonal
matrix Q of size 3,000,000 by 3,000,000 and R an upper triangular 500
by 500 matrix but stored in essentially the same amount of space as
the original matrix X.

> I could break up the matrices into
> smaller chunks and do the multiplication, but Z is later used in several
> other multi-step calculations (with addition and multiplication mostly) so
> I think that would be a last resort.  If I can get the 64-bit matrix
> indices working in RcppArmadillo, I think that will solve much of the
> problem, because I will only need to return very long vectors and not big
> matrices.
>
> Joshua
> --
> Joshua French, Ph.D.
> Assistant Professor
> Department of Mathematical and Statistical Sciences
> University of Colorado Denver
> joshua.fre...@ucdenver.edu
> http://math.ucdenver.edu/~jfrench/
> Ph:  303-556-6265  Fax:  303-556-8550
>
>
>
>
>
>
>
> On 7/17/12 8:56 AM, "Douglas Bates"  wrote:
>
>>On Tue, Jul 17, 2012 at 8:14 AM, Dirk Eddelbuettel  wrote:
>>>
>>> On 16 July 2012 at 23:30, French, Joshua wrote:
>>> | I am doing some linear algebra on large matrices in R and receiving
>>>the
>>> | following error:  "allocMatrix: too many elements specified".  From
>>>what I
>>> | understand, the error is caused by the fact that R uses 32-bit ints
>>>and not
>>> | 64-bit ints for matrix indices, so R doesn't have a way to represent
>>>all the
>>> | elements in the very large matrix.
>>> |
>>> | My two questions:
>>> |
>>> | 1.  Armadillo (and presumably RcppArmadillo) will not have this issue
>>>since
>>> | Armadillo provided support for 64-bit indices as of version 2.4.0.
>>>Is there a
>>> | way to easily utilize this functionality from within RcppArmadillo?
>>>
>>> I need to double check but this may have been a compile-time option you
>>>need
>>> to enable. In any event ... R indices are still limited so you may not
>>>be
>>> able to pass these back and forth.
>>>
>>> | 2.  I have found in the past that some of the speeds gains from
>>>RcppArmadillo
>>> | in comparison to pure R are lost when passing large matrices as
>>>arguments.
>>> |  There will always be overhead when passing arguments (especially
>>>large matrix
>>> | arguments) to pretty much any function.  Are there any tricks to
>>>minimize the
>>> | overhead when passing a non-sparse matrix argument of say 1,000,000
>>>by 500 from
>>> | R to Armadillo?
>>>
>>> I defer all question concerning sparse matrices to Doug and other users
>>>of
>>> sparse matrix code.  I live mostly in a small-to-medium size dense
>>>matrix world.
>>
>>Actually the question was about non-sparse matrices.  It looks as if
>>it is the number of rows in the matrix that will be problematic.  An
>>upper bound on the number of rows is the maximum integer value divided
>>by the number of columns.
>>
>>> .Machine$integer.max / 500
>>[1] 4294967
>>
>>I would try not to exceed about 1/2 to 1/4 of that bound.
>>
>>A simple way of handling data sets that are much larger than that is
>>to work with a sample of the rows.  If that is not feasible then I
>>would create a list of matrices each of size 1,000,000 by 500 or so
>>and vertically concatenate them in the  C++ code.  Of course, this
>>means a copying operation.  Also, when you are finished if you need to
>>pass results back to R then you face a similar problem getting a large
>>matrix in C++ back into R storage.
>>
>>You can create a read-only matrix in C++ using the storage from R as
>>described by Christian for RcppArmadillo or using the Eigen::Map class
>>in RcppEigen.
>>
>>What are you (Joshua) doing with these large matrices?  If the main
>>calculations involve X'X-type calculations you can carry out the
>>calculations on horizontal chunks and then assemble the results.
>
___
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

Re: [Rcpp-devel] Very Large Matrices in RcppArmillao

[Douglas Bates]

On Tue, Jul 17, 2012 at 8:14 AM, Dirk Eddelbuettel  wrote:
>
> On 16 July 2012 at 23:30, French, Joshua wrote:
> | I am doing some linear algebra on large matrices in R and receiving the
> | following error:  "allocMatrix: too many elements specified".  From what I
> | understand, the error is caused by the fact that R uses 32-bit ints and not
> | 64-bit ints for matrix indices, so R doesn't have a way to represent all the
> | elements in the very large matrix.
> |
> | My two questions:
> |
> | 1.  Armadillo (and presumably RcppArmadillo) will not have this issue since
> | Armadillo provided support for 64-bit indices as of version 2.4.0.  Is 
> there a
> | way to easily utilize this functionality from within RcppArmadillo?
>
> I need to double check but this may have been a compile-time option you need
> to enable. In any event ... R indices are still limited so you may not be
> able to pass these back and forth.
>
> | 2.  I have found in the past that some of the speeds gains from 
> RcppArmadillo
> | in comparison to pure R are lost when passing large matrices as arguments.
> |  There will always be overhead when passing arguments (especially large 
> matrix
> | arguments) to pretty much any function.  Are there any tricks to minimize 
> the
> | overhead when passing a non-sparse matrix argument of say 1,000,000 by 500 
> from
> | R to Armadillo?
>
> I defer all question concerning sparse matrices to Doug and other users of
> sparse matrix code.  I live mostly in a small-to-medium size dense matrix 
> world.

Actually the question was about non-sparse matrices.  It looks as if
it is the number of rows in the matrix that will be problematic.  An
upper bound on the number of rows is the maximum integer value divided
by the number of columns.

> .Machine$integer.max / 500
[1] 4294967

I would try not to exceed about 1/2 to 1/4 of that bound.

A simple way of handling data sets that are much larger than that is
to work with a sample of the rows.  If that is not feasible then I
would create a list of matrices each of size 1,000,000 by 500 or so
and vertically concatenate them in the  C++ code.  Of course, this
means a copying operation.  Also, when you are finished if you need to
pass results back to R then you face a similar problem getting a large
matrix in C++ back into R storage.

You can create a read-only matrix in C++ using the storage from R as
described by Christian for RcppArmadillo or using the Eigen::Map class
in RcppEigen.

What are you (Joshua) doing with these large matrices?  If the main
calculations involve X'X-type calculations you can carry out the
calculations on horizontal chunks and then assemble the results.
___
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

Re: [Rcpp-devel] replacing a splines ns() call with other Rcpp attached library function

[Douglas Bates]

On Fri, Jun 29, 2012 at 6:18 AM, Silkworth,David J.
 wrote:
> With respect to replies:
>  > a) fit a decent, modern, documented C++ implementation of what is needed
>
> Yes, this was the question.  Can anyone help?  I am failing right here.
>
>> and of course if a) fails, one can always 'rip out' what is underneath a
>> known R function and use that as a fallback.
>
> This is where I am.  I can easily (I think) 'rip out' kind Douglas Bate's two 
> C functions:
>
> /* Exports */
> SEXP spline_basis(SEXP knots, SEXP order, SEXP xvals, SEXP derivs);
> SEXP spline_value(SEXP knots, SEXP coeff, SEXP order, SEXP x, SEXP deriv);
>
> By initially just placing the entire contents of splines.c in the inline 
> header for testing.
>
> The challenge now is to send one or both of these the correct arguments, and 
> somehow come up with an order  (3 for the cubic spline) by length(xvals) 
> matrix from what appears to be individual vectors returning from calls here.
>
> Since the entire topic of the boor algorithm is new to me just this week this 
> remains a daunting challenge.  I plan to single step through the R code of 
> ns() and splineDesign(), which it calls, just to see what args are what and 
> what return objects are what.

His name is "de Boor", as in http://en.wikipedia.org/wiki/Carl_R._de_Boor
> We just can't stand for code in the middle of a double loop to go to R just 
> to get all the "boiler-plate" of the R code in ns().  We need to know how to 
> call spline_basis directly with confidence.
>
> I'm just looking to get one function to work for one problem.  No thought at 
> this time of actually attempting something like building an "Rf_ns" function 
> to bind into Rcpp.
>
> I suppose if one were to go that far, after understanding what I must, 
> Davor's black art method might be interesting.  But that is a long way off 
> from where I stand today.  We can rip code into one thesis project and be 
> done I hope.
>
>
>
>
>
___
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

Re: [Rcpp-devel] replacing a splines ns() call with other Rcpp attached library function

[Douglas Bates]

On Thu, Jun 28, 2012 at 8:08 AM, Silkworth,David J.
 wrote:
> Folks, I am bravely helping Gianluca Bonitta (a.k.a. bbo...@tin.it) develop
> an Rcpp solution for his thesis, he already has a functioning R prototype.

> He has been using ns() from the splines library (thank-you Doug Bates!), and
> we are interested in identifying an alternative call to one of the libraries
> that are already available to Rcpp. GSL seemed most likely to me, but this
> topic is a bit over my head, so I would appreciate some guidance.

Actually, I didn't write ns().  I only wrote the B-spline basis
functions and evaluations.  There is C source code in
$RSRC/library/splines/src/splines.c although it is not terribly easy
to understand - I was in a "terser is better" phase when that was
written.

You may find it easier to interface to some of the more recent C++
libraries (e.g. libnurbs.sourceforge.net) for the evaluation.
___
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

Re: [Rcpp-devel] check validity of object created from Rcpp module

[Douglas Bates]

On Tue, Jun 26, 2012 at 7:04 PM, Dirk Eddelbuettel  wrote:
>
> On 26 June 2012 at 17:21, Jiqiang Guo wrote:
> | Dear List,
> |
> | I am wondering if there is a way to check whether an object created from 
> Rcpp's
> | module is valid or not especially for objects loaded in another session.
> |
> | In detail, we know that we could not save and load objects created from 
> Rcpp's
> | module [section 5 of http://cran.r-project.org/web/packages/Rcpp/vignettes/
> | Rcpp-modules.pdf].  But that does not prevent users doing that. So if an
> | object created from Rcpp's module is loaded in another session, we would 
> end up
> | with errors or even segfault.  It would be nice we have a method to check if
> | it is valid or not so that when we are writing a package, we could do a 
> check
> | before calling those methods of that object.
>
> It's a very good idea. I don't have a good idea as to how we'd do this though.

There has to be one or more external pointers in the R object
somewhere (although I spent some time looking at the Rcpp-modules code
and came away very confused so I can't tell you exactly where.)  That
pointer will become a NULL pointer after a serialize/unserialize
sequence so you need to check for its being NULL.  This is why the
reference classes in the development version of lme4 have a field
named ptr and a method named Ptr.  The Ptr method first checks if ptr
is a null pointer and, if so, regenerates the object then returns ptr.
 I can say from experience that this all gets rather complicated.
However I haven't looked at John Chambers' new code for Rcpp classes
in R and he might have addressed the issue.

>
> Dirk
>
> | An illustrative example:
> |
> | require("Rcpp")
> | require("inline")
> |
> | inc <- '
> | class Foo {
> | private:
> |   int x_;
> | public:
> |   Foo(int x) : x_(x) { }
> |   int getx() const { return x_;}
> | };
> |
> | RCPP_MODULE(foo) {
> |   class_("Foo")
> |   .constructor()
> |   .method("getx", &Foo::getx)
> |   ;
> | }
> |
> | fx <- cxxfunction(signature(), "", include = inc, plugin = "Rcpp", verbose =
> | FALSE)
> | m <- Module("foo", getDynLib(fx))
> | cppo  = new(m$Foo, 5)
> |
> | print(cppo$getx())
> | save("cppo", file = 'cppo.RData')
> |
> | rm(list = ls())
> | load("cppo.RData")
> | print(cppo$getx())
> |
> | For the above code, we know the last line would not work.  If we have a
> | function to test the validity of cppo,then if not valid, we would not call
> | cppo's function and are able to tell users that cppo is not valid any more.
> |
> | Best,
> | Jiqiang
> |
> | --
> | ___
> | 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
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> ___
> 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
___
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

Re: [Rcpp-devel] Persistent C++ objects within R/Rcpp?

[Douglas Bates]

On Mon, Jun 25, 2012 at 3:38 PM, Andy Garcia  wrote:
> First, thank you to the Rcpp team for their work on Rcpp. It's an incredible
> resource to get C++ code executing from R.
>
> I've worked through setting up the correct development environment (Ian
> Fellow's post on Eclipse + Rcpp was great with a few tweaks for a Linux
> install) as well some of the examples. The next step I'd like to take is
> develop some prototype code for implementing an individual-based simulation.
>
> Let's say I have a function that creates the simulated population based on
> some number of input parameters:
>
> RcppExport SEXP CreatePopulation( /* some amount of input parameters*/ );
>
> I would then want other functions called in R via Rcpp to operate on the
> population or return specific results. The part I'm struggling with is how
> would the population state persist when created in CreatePopulation to be
> then access/modified in other functions. My initial though is to have some
> sort of population object, however I don't know how this would persist when
> the code returns from execution in C++. R can't store a C++ object to be
> passed back, correct?

It can, in a way.  Look up the Rcpp::ExternalPtr class.  You pass back
a pointer to the C++ object then pass it back down in calls to other
C++ code.  The thing you must be aware of is that
serialization/unserialization (i.e. save/load) on the R object will,
naturally, not be able to restore your C++ object so the external
pointer comes back as a NULL pointer.

> Could I just create a global object within C++?
>
> Any feedback would be appreciated,
> Andy
>
> ___
> 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
___
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

Re: [Rcpp-devel] Using sparse matrices from the Matrix package in Rcpp

[Douglas Bates]

On Sun, Jun 24, 2012 at 1:19 PM, Dirk Eddelbuettel  wrote:
>
> On 24 June 2012 at 19:56, Glenn Lawyer wrote:
> | I would suggest using the eigen library, Rccp has an interface to this.
> | Alternately, you could link to the boost graph library.
>
> Spot on.
>
> It was in fact the availability of (more) sparse matrix methods (even though
> they were (are ?) under development and not fully mature) that lead Doug from
> Armadillo to Eigen. I am sure he will follow-up here too.

Indeed.  There is good support for sparse matrices in Eigen including
the use of the Cholmod functions from the Matrix package, although the
latter requires a bit of effort.

With the release of Eigen 3.1 today, sparse matrices are now a
supported part of Eigen.  I have not yet incorporated Eigen 3.1 into
RcppEigen but will do so.

> Boost Graph is interesting as CRAN had the RBGL package to access it for such
> a long time.  It would be nice if someone looked more closely if we can in
> fact add appropriate as<>() and wrap() methods/functions to make interfacing
> this library more easy / more accessible from R.
>
> Dirk
>
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> ___
> 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
___
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

Re: [Rcpp-devel] Call R function

[Douglas Bates]

On Wed, Jun 20, 2012 at 4:30 PM, bbo...@tin.it  wrote:
>
> i must use it to evaluate  matropolis ratio (green ratio) in the following
> code in R :  i have  start just  now to study c and c++   (like you can see
> by my stupids questions (for you surely))  if  you have some suggestion im
> very happy to "heard" all suggestion ..thank for all  of the mailing
> list ...that teach me much
> tahnk you again and again
> OneTimeVariable<-function (logpost, start, n.iter, burn, thin, scale, ...)
> {
>     pb <- txtProgressBar(min = 0, max = n.iter, style = 3)
>     p = length(start)
>     vth = array(0, dim = c(n.iter/thin, p))
>     f0 = logpost(start, ...)
>     arate = array(0, dim = c(1, p))
>     th0 = start
>     for (i in (-burn):n.iter) {
>     setTxtProgressBar(pb, i)
>     for (j in 1:p) {
>     th1 = th0
>     th1[j] = th0[j] + rnorm(1) * scale[j]
>     f1 = logpost(th1, ...)
>     u = runif(1) <= min(exp(f1 - f0),1)
>     th0[j] = th1[j] * (u == 1) + th0[j] * (u == 0)
>     f0 = f1 * (u == 1) + f0 * (u == 0)
>     vth[floor(i/thin), j] = th0[j]
>     arate[j] = arate[j] + u
>     }
>     }
>     arate = arate/n.iter
>     stuff = list(par = mcmc(vth), accept = arate)
>     return(stuff)
> }

In the enclosed file I provide two alternative versions in R and one
using Rcpp but none of these are tested.  I suggest that you create a
test case and first find out if they produce the same distribution of
results (OneTimeV1 should be faster but the results will be different,
even after calling set.seed with the same seed).  Then use the
benchmark package to try see the relative speed.


foo.R
Description: Binary data
___
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

Re: [Rcpp-devel] Call R function (bbo...@tin.it)

[Douglas Bates]

On Wed, Jun 20, 2012 at 11:22 AM, Silkworth,David J.
 wrote:
> I am not going to ask why to do this, but here is a way how.
>
> I added your R function into the mypackage skeleton from RcppPackage
>
> R_user_F<-function (par) {
> y<-par[1]
> x<-par[2]
> return(dnorm(x)*dnorm(y))
> }
>
> I prefer to use the explicit return function for clarity.
>
> Now I install the package after a binary build.  It has become a part of R.
>
> The code you wish that calls this function is as follows:
>
> require(mypackage)
> require(inline)
>
> src <- '
> Rcpp::NumericVector par=(arg1);
> Rcpp::NumericVector par1(2);
> Rcpp::NumericVector par2(2);
> Rcpp::NumericVector f1(1);
> Rcpp::NumericVector f2(1);
> Rcpp::NumericVector outval(1);
> double a=2;
> double b=3;
> par1=par*b;
> par2=par*a;
> Environment mypackage("package:mypackage");
> Function R_userf = mypackage["R_user_F"];
> f1=R_userf(par1);
> f2=R_userf(par2);
> bool u = (Rf_runif(0.0,1.0) <= f2[0]-f1[0]  );
> outval[0]=u;
> return(outval);
> '
>
>  RcppGibbs  <- cxxfunction(signature(arg1 = "numeric"),
>  src, plugin = "Rcpp")
>
> I had to assume that you wished to return the Boolean result.  There is no 
> Rcpp::BooleanVector, so we have to use the knowledge that a Boolean in C++ is 
> just an integer 0 for false and 1 for true.

There is an Rcpp::LogicalVector (R uses the name logical instead of
boolean) or you can just call Rcpp::wrap on a bool value, as shown in
the code I posted.  The point is that, when passed back to R, the
elements in an Rcpp::LogicalVector are displayed as TRUE and FALSE.
An alternative is

return ::Rf_ScalarLogical(u);

which is a short-cut for creating an R logical vector of length 1
> You will have to run a quick test on the return value to get the Boolean in R.
>
> Notice that you have to dimension your Rcpp::NumericVector on declaration.
> You must address the value you wish when making calculations so f2-f1, 
> becomes f2[0]-f1[0] in this case.
>
> ___
> 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
___
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

Re: [Rcpp-devel] Call R function

[Douglas Bates]

I enclose a rewrite of your function

R version 2.15.0 (2012-03-30)
Copyright (C) 2012 The R Foundation for Statistical Computing
ISBN 3-900051-07-0
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> require(inline)
Loading required package: inline
> code <- '
+ Function  fun(ff);
+ NumericVector par(pp);
+ NumericVector par1(par), par2(par);
+ const double a(2.), b(3.);
+ 
+ par1[0] *= a;
+ par2[0] *= b;
+ NumericVector f1(fun(par1)), f2(fun(par2));
+ return wrap(Rf_runif(0., 1.) <= (f1[0] - f2[0]));
+ '
> RcppGibbs <- cxxfunction(signature(pp = "numeric", ff="function"),
+  code, plugin="Rcpp")
> ff <- function(par) dnorm(par[1]) * dnorm(par[2])
> 
> RcppGibbs(3., ff)
[1] FALSE
> 
> proc.time()
   user  system elapsed 
  4.888   0.388   5.390 


foo.R
Description: Binary data
___
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

Re: [Rcpp-devel] Call R function

[Douglas Bates]

On Wed, Jun 20, 2012 at 12:52 AM, bbo...@tin.it  wrote:
> Hello i would like to compute the following code ;
>
> Where R_user_F is a R function  (i could write this function directly inin
> line but is only for example because i want to call an R function more
> complex):
> so when compile appear this error  message:  error:
> "cannot convert 'Rcpp::sugar::Comparator_With_One_Value<14,
> Rcpp::sugar::greater_or_equal<14>, true,
> Rcpp::sugar::Minus_Vector_Vector<14, true, Rcpp::Vector<14>, true,
> Rcpp::Vector<14> > >' to 'bool' in initialization."
>
>  I can do for solve this problem? Thank You.
>
> The code:
>
> R_user_F<-function (par) {
> y<-par[1]
> x<-par[2]
> dnorm(x)*dnorm(y)
> }
>
>
>
> require(inline)
> code <- '
> NumericVector f1;
> NumericVector f2;
> NumericVector par1;
> NumericVector par2;
> NumericVector a=2;
> NumericVector b=3;
> NumericVector par =par;
> par1=par*b;
> par2=par*a;
> Function R_userf(fun);
> f1=R_userf(par1);
> f2=R_userf(par2);
> bool u = (Rf_runif(0.0,1.0) <= f2-f1  );
>
> '
> RcppGibbs <- cxxfunction(signature( par ="NumericVector", fun="function"),
> code,
> include='#include ',
> plugin="Rcpp")<-function (par) {
> y<-par[1]
> x<-par[2]
> dnorm(x)*dnorm(y)
> }

I don't understand the last call.  In cxxfunction the signature
argument uses the names of R classes, not Rcpp classes, so the
specification should be par="numeric".  Also, the second assignment
(RcppGibbs <- cxxfunction(...) <- function(par)) would be an error
because you can't assign to the cxxfunction call.

The error message relates to the expression Rf_runif(0.0, 1.0) <= f1 -
f2 because Rf_unif(0.0, 1.0) returns a double and f1 - f2 is a
NumericVector.  You should use (f1-f2)[0] to convert to a double.

With regard to your code, you are writing C++ in a C style where you
declare all the variables then assign them.  That isn't the best
approach in C++ because default constructors can be expensive. I
always find it easier to use constructors in declarations rather than
declarations followed by assignment in C++.  If you use a declaration
like

NumericVector f1;

then later assign

f1 = R_userf(par1)

you end up creating a NumericVector of length 0 then creating another
NumericVector of length > 0 to overwrite it.

I would begin by assigning the values of the arguments to Rcpp
objects.  I'll post a revised version shortly.
___
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

Re: [Rcpp-devel] Differences between RcppEigen and RcppArmadillo

[Douglas Bates]

On Sun, Jun 17, 2012 at 11:09 PM, c s  wrote:
> On Sun, Jun 17, 2012 at 4:21 AM, Douglas Bates wrote:
>> These comments may provoke a heated response from Conrad but,
>> if so, I don't plan to respond further.  Eigen and Armadillo are different
>> approaches, each with their own advantages and disadvantages.
>
> I'm actually in favour of heterogeneous solutions and support the use
> of C++ for performance critical implementations.
>
> Relying on only one solution (ie. homogeneity) is bad from a
> "survival" point of view, due to very similar reasons as in the
> biological world: a fault can affect the entire population, or a virus
> can propagate itself quite easily and do lots of nasty things.  For
> example, see the recent snafu with state-sponsored viruses attacking
> Windoze installations, ie, Stuxnet and Flame:
> http://www.h-online.com/security/features/FAQ-Flame-the-super-spy-1587063.html
>
> In my view, the C++ brethren Armadillo and Eigen are not really up
> against each other, but against Matlab.  Matlab is like crack or
> crystal meth.  Highly addictive, but in the end bad for you.
> Mathworks, the company that produces Matlab, hooks people in while
> they are university students.  It provides low prices to educational
> institutions, and charges through the roof for everyone else.
> Students become graduates, which are then employed by the companies
> that are in effect stuck with using Matlab.  This is a strategy known
> as vendor lock-in, exploited by monopolists.  Other descriptions
> consistent with the business model of Mathworks are "parasite" and
> "leech".

I completely and happily agree.

> One way of keeping Mathworks in check is to provide alternative
> solutions that have similar functionality (eg. easy visualisation of
> data and rapid prototyping).  One of them in GNU Octave, which really
> needs a JIT compiler to effectively compete with Matlab..  Another is
> R.  I believe R currently doesn't have a JIT compiler (I haven't
> checked lately), and hence the very useful Rcpp fills in the
> performance gap.

The recently-developed Julia language (JuliaLang.org) has a JIT and a
syntax that is familiar to Matlab and R users.  It is still very much
a work in progress but shows great promise.

> One also has to address the current reality that people still overly
> rely on Matlab, and have the problem of converting their Matlab code
> into production environments / products.  This is where C++ comes in,
> and libraries such as Armadillo (with its Matlab-like syntax) can be
> quite handy.
___
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

Re: [Rcpp-devel] Differences between RcppEigen and RcppArmadillo

[Douglas Bates]

On Thu, Jun 14, 2012 at 2:36 PM, Dirk Eddelbuettel  wrote:
>
> On 15 June 2012 at 02:56, c s wrote:
> | Simply installing ATLAS (which provides speed-ups for several Lapack
> | functions) on Debian/Ubuntu systems can already make a big difference.
> |  (Debian & Ubuntu use a trick to redirect Lapack and Blas calls to
> | ATLAS).  Under Mac OS X, the Accelerate framework provides fast
> | implementations of Lapack and Blas functions (eg. using
> | multi-threading).
>
> I found OpenBLAS to be faster than Atlas (with both coming from the
> distribution), and tend to install that now. It uses a multithreaded approach
> by default, but I have to double check one thing about cpu affinity which,
> when used from R, has been seen to interfere.  That is possible reason for
> the single-threaded performance here. I'll report back.
>
> | I've taken the modular approach to Armadillo (ie. using Lapack rather
> | than reimplementing decompositions), as it specifically allows other
> | specialist parties (such as Intel) to provide Lapack that is highly
> | optimised for particular architectures.  I myself would not be able to
> | keep up with the specific optimisations required for each CPU.  This
> | also "future-proofs" Armadillo for each new CPU generation.
> |
> | More importantly, numerically stable implementation of computational
> | decompositions/factorisations is notoriously difficult to get right.
> | The core algorithms in Lapack have been evolving for the past 20+
> | years, being exposed to a bazillion corner-cases.  Lapack itself is
> | related to Linpack and Eispack, which are even older.  I've been
> | exposed to software development long enough to know that in the end
> | only time can shake out all the bugs.  As such, using Lapack is far
> | less risky than reimplementing decompositions from scratch.  A
> | "home-made" matrix decomposition might be a bit faster on a particular
> | CPU, but you have far less knowledge as to when it's going to blow up
> | in your face.
> |
> | High-performance variants of Lapack, such as MKL, take an existing
> | proven implementation of a decomposition algorithm and recode parts of
> | it in assembly, and/or parallelise other parts.
>
> All excellent points which nobody disputes.  It just so happens that Eigen
> still looks better in some benchmarks but as you say we have to ensure we
> really compare apples to apples.

While I realize that Conrad has strong opinions on these issues and is
unlikely to be moved from those opinions, I would make a few points in
favor of Eigen.

- using BLAS/Lapack is effective because that is mature, well-tested
code.  It is ineffective because the reference implementations of
BLAS/Lapack are in Fortran 77 and suffer from all the idiosyncrasies
of that ancient language.  These include call-by-reference semantics,
even for scalars, lack of structures beyond pointers to arrays and the
inability to allocate storage.  Hence, the caller frequently must pass
an array 'work' and it length 'lwork' and, occasionally, 'rwork' and
'iwork'. There is a helpful convention for getting the length of the
work array correct which involves a preliminary call with lwork = -1
that does nothing but write the desired work array length into the
first element of the work array.  The caller must then allocate the
work array, change the value of lwork and call the routine a second
time.  Because the only structure available is a pointer to a
one-dimensional array, you must pass a general matrix as four
arguments M, N, A, LDA.  There are no machine-readable header files.
Jeff Bezanson recently joked about writing a parser for the upper-case
Fortran comments to create headers.  Does any of this sound like 21st
century software design?

- MKL blas is very good on Intel hardware.  Naturally Intel doesn't go
to great lengths to support AMD processors to the same extent.  MKL is
proprietary, closed-source software and plugging it into an
open-source software systems like R and Aramadillo feels like a
compromise.  There are open alternatives like Atlas and OpenBlas.
However, these all need to conform to the Fortran calling sequence
specifications (Lapacke and CBLAS do provide C calling sequences but
it is not a good idea to count on their availability at present).  And
notice that this approach contradicts the "Lapack code is complex but
widely tested so it should be left untouched" argument.  The way to
get good performance in OpenBlas and MKL is to rewrite the BLAS and
heavily-used parts of Lapack, like dgetrf, in a combination of C and
Assembler.  I would also point out that one of the reasons that Lapack
code is complex and difficult to understand is because it is written
in the archaic language Fortran 77.

- the approach in Eigen was to use templated C++ classes for dense and
sparse matrices and for the decompositions and solvers.  I don't think
readers of this group need to be reminded of the advantages of
templated object-oriented C++ code.  Eigen has a bit of a le

Re: [Rcpp-devel] Differences between RcppEigen and RcppArmadillo

[Douglas Bates]

On Wed, Jun 13, 2012 at 6:53 PM, Julian Smith  wrote:
> Doesn't svd in R by default compute D, U and V?

> http://stat.ethz.ch/R-manual/R-patched/library/base/html/svd.html

You're right but the default is the 'thin' U when X is n by p and n >=
p.  Does the svd in Armadillo return the full n by n matrix U?

Another thing to check is what underlying Lapack routine is called.
There are two: dgesvd, which is the older algorithm and the one with
the expected name if you know the Lapack conventions, and dgesdd which
is a newer and faster algorithm.  R uses dgesdd.

> On Wed, Jun 13, 2012 at 4:07 PM, Douglas Bates  wrote:
>>
>> On Wed, Jun 13, 2012 at 5:16 PM, Dirk Eddelbuettel  wrote:
>> >
>> > On 13 June 2012 at 15:05, Julian Smith wrote:
>> > | I agree that RcppEigen is a little bit faster, but ease of use is
>> > important to
>> > | me, so I feel like RcppArmadillo might win out in my application.
>> >
>> > Yup, that my personal view too.
>> >
>> > | | RcppArmadillo will use the very same LAPACK and BLAS libs your R
>> > session
>> > | | uses. So MKL, OpenBlas, ... are all options.  Eigen actually has its
>> > own
>> > | code
>> > | | outperforming LAPACK, so it doesn't  as much there.
>> > |
>> > | Why do you think R outperforms RcppArmadillo in this example below?
>> > Anyway to
>> > | speed this up?
>> >
>> > That is odd. "I guess it shouldn't." I shall take another look -- as I
>> > understand it both should go to the same underlying Lapack routine.  I
>> > may
>> > have to consult with Conrad on this.
>> >
>> > Thanks for posting a full and reproducible example!
>> >
>> > Dirk
>> >
>> > | require(RcppArmadillo)
>> > | require(inline)
>> > |
>> > | arma.code <- '
>> > |   using namespace arma;
>> > |   NumericMatrix Xr(Xs);
>> > |   int n = Xr.nrow(), k = Xr.ncol();
>> > |   mat X(Xr.begin(), n, k, false);
>> > |   mat U;
>> > |   vec s;
>> > |   mat V;
>> > |   svd(U, s, V, X);
>> > |   return wrap(s);
>> > | '
>>
>> Because the arma code is evaluating the singular vectors (U and V) as
>> well as the singular values (S) whereas the R code is only evaluating
>> the singular values.  There is considerably more effort required to
>> evaluate the singular vectors in addition to the singular values.
>>
>> > | rcppsvd <- cxxfunction(signature(Xs="numeric"),
>> > |     arma.code,
>> > |     plugin="RcppArmadillo")
>> > |
>> > | A<-matrix(rnorm(5000^2), 5000)
>> > |
>> > | > system.time(rcppsvd(A))
>> > |     user   system  elapsed
>> > | 1992.406    4.862 1988.737
>> > |
>> > | > system.time(svd(A))
>> > |    user  system elapsed
>> > | 652.496   2.641 652.614
>> > |
>> > | On Wed, Jun 13, 2012 at 11:43 AM, Dirk Eddelbuettel 
>> > wrote:
>> > |
>> > |
>> > |     On 13 June 2012 at 10:57, Julian Smith wrote:
>> > |     | I've been toying with both RcppArmadillo and RcppEigen the past
>> > few days
>> > |     and
>> > |     | don't know which library to continue using. RcppEigen seems
>> > really slick,
>> > |     but
>> > |     | appears to be lacking some of the decompositions I want and
>> > isn't nearly
>> > |     as
>> > |     | fast to code. RcppArmadillo seems about as fast, easier to code
>> > up etc.
>> > |     What
>> > |     | are some of the advantages/disadvantages of both?
>> > |
>> > |     That's pretty close.  I have been a fan of [Rcpp]Armadillo which I
>> > find
>> > |     easier to get my head around.  Doug, however, moved from
>> > [Rcpp]Armadillo
>> > |     to
>> > |     [Rcpp]Eigen as it has some things he needs.  Eigen should have a
>> > "larger"
>> > |     API
>> > |     than Armadillo, but I find the code and docs harder to navigate.
>> > |
>> > |     And you should find Eigen to be a little faster. Andreas Alfons
>> > went as far
>> > |     as building 'robustHD' using RcppArmadillo with a drop-in for
>> > RcppEigen (in
>> > |     package 'sparseLTSEigen'; both package names from memmory and

Re: [Rcpp-devel] Differences between RcppEigen and RcppArmadillo

[Douglas Bates]

On Wed, Jun 13, 2012 at 5:16 PM, Dirk Eddelbuettel  wrote:
>
> On 13 June 2012 at 15:05, Julian Smith wrote:
> | I agree that RcppEigen is a little bit faster, but ease of use is important 
> to
> | me, so I feel like RcppArmadillo might win out in my application.
>
> Yup, that my personal view too.
>
> | | RcppArmadillo will use the very same LAPACK and BLAS libs your R session
> | | uses. So MKL, OpenBlas, ... are all options.  Eigen actually has its own
> | code
> | | outperforming LAPACK, so it doesn't  as much there.
> |
> | Why do you think R outperforms RcppArmadillo in this example below? Anyway 
> to
> | speed this up?
>
> That is odd. "I guess it shouldn't." I shall take another look -- as I
> understand it both should go to the same underlying Lapack routine.  I may
> have to consult with Conrad on this.
>
> Thanks for posting a full and reproducible example!
>
> Dirk
>
> | require(RcppArmadillo)
> | require(inline)
> |
> | arma.code <- '
> |   using namespace arma;
> |   NumericMatrix Xr(Xs);
> |   int n = Xr.nrow(), k = Xr.ncol();
> |   mat X(Xr.begin(), n, k, false);
> |   mat U;
> |   vec s;
> |   mat V;
> |   svd(U, s, V, X);
> |   return wrap(s);
> | '

Because the arma code is evaluating the singular vectors (U and V) as
well as the singular values (S) whereas the R code is only evaluating
the singular values.  There is considerably more effort required to
evaluate the singular vectors in addition to the singular values.

> | rcppsvd <- cxxfunction(signature(Xs="numeric"),
> |     arma.code,
> |     plugin="RcppArmadillo")
> |
> | A<-matrix(rnorm(5000^2), 5000)
> |
> | > system.time(rcppsvd(A))
> |     user   system  elapsed
> | 1992.406    4.862 1988.737
> |
> | > system.time(svd(A))
> |    user  system elapsed
> | 652.496   2.641 652.614
> |
> | On Wed, Jun 13, 2012 at 11:43 AM, Dirk Eddelbuettel  wrote:
> |
> |
> |     On 13 June 2012 at 10:57, Julian Smith wrote:
> |     | I've been toying with both RcppArmadillo and RcppEigen the past few 
> days
> |     and
> |     | don't know which library to continue using. RcppEigen seems really 
> slick,
> |     but
> |     | appears to be lacking some of the decompositions I want and isn't 
> nearly
> |     as
> |     | fast to code. RcppArmadillo seems about as fast, easier to code up 
> etc.
> |     What
> |     | are some of the advantages/disadvantages of both?
> |
> |     That's pretty close.  I have been a fan of [Rcpp]Armadillo which I find
> |     easier to get my head around.  Doug, however, moved from [Rcpp]Armadillo
> |     to
> |     [Rcpp]Eigen as it has some things he needs.  Eigen should have a 
> "larger"
> |     API
> |     than Armadillo, but I find the code and docs harder to navigate.
> |
> |     And you should find Eigen to be a little faster. Andreas Alfons went as 
> far
> |     as building 'robustHD' using RcppArmadillo with a drop-in for RcppEigen 
> (in
> |     package 'sparseLTSEigen'; both package names from memmory and I may have
> |     mistyped).  He reported a performance gain of around 25% for his problem
> |     sets.  On the 'fastLm' benchmark, we find the fast Eigen-based
> |     decompositions
> |     to be much faster than Armadillo.
> |
> |     | Can you call LAPACK or BLAS from either? Is there a wrapper in 
> RcppEigen
> |     to
> |     | call LAPACK functions? Want some other decomposition methods, dont 
> like
> |     the
> |     | JacobiSVD method in Eigen.
> |
> |     You need to differentiate between the Eigen and Armadillo docs _for 
> their
> |     libraries_ and what happens when you access the Rcpp* variant from R.
> |
> |     RcppArmadillo will use the very same LAPACK and BLAS libs your R session
> |     uses. So MKL, OpenBlas, ... are all options.  Eigen actually has its own
> |     code
> |     outperforming LAPACK, so it doesn't  as much there.
> |
> |     Hope this helps,   Dirk (at useR!)
> |
> |     |
> |     | --
> |     | ___
> |     | 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
> |     --
> |     Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> |
> |
>
> --
> Dirk Eddelbuettel | e...@debian.org | http://dirk.eddelbuettel.com
> ___
> 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
___
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

[Rcpp-devel] Memory error when using Dimension

[Douglas Bates]

Some time ago I mentioned memory errors that seemed to be related the
the Rcpp::Dimension class but I couldn't pin them down.  I have a more
detailed, although not necessarily reproducible, example now.

In the merPredD::condVar method defined in the lme4/src/predModule.cpp
source file there is the line

Rcpp::NumericVector ansi(Rcpp::Dimension(ncti, ncti, nli));

When that method is eventually invoked from a version of lme4 compiled
with -ltcmalloc for R-2.15.0 (Ubuntu Version: 2.15.0-1precise0)
tcmalloc reports an attempt to free an invalid pointer.  The gdb
traceback is

> print(dotplot(ranef(fm1, postVar = TRUE), strip = FALSE)[[1]])
src/tcmalloc.cc:390] Attempt to free invalid pointer: 0x53d75d0

Program received signal SIGABRT, Aborted.
0x7726f445 in raise () from /lib/x86_64-linux-gnu/libc.so.6
(gdb) where
#0  0x7726f445 in raise () from /lib/x86_64-linux-gnu/libc.so.6
#1  0x77272bab in abort () from /lib/x86_64-linux-gnu/libc.so.6
#2  0x70e507fb in ?? () from /usr/lib/libtcmalloc.so.0
#3  0x70e50b04 in TCMalloc_CrashReporter::PrintfAndDie(char
const*, ...) () from /usr/lib/libtcmalloc.so.0
#4  0x70e45840 in ?? () from /usr/lib/libtcmalloc.so.0
#5  0x70e64191 in tc_delete () from /usr/lib/libtcmalloc.so.0
#6  0x707fcde6 in deallocate (__p=,
this=) at /usr/include/c++/4.6/ext/new_allocator.h:98
#7  _M_deallocate (__p=, this=,
__n=) at /usr/include/c++/4.6/bits/stl_vector.h:156
#8  ~_Vector_base (this=, __in_chrg=) at
/usr/include/c++/4.6/bits/stl_vector.h:142
#9  ~vector (this=0x7fff9620, __in_chrg=) at
/usr/include/c++/4.6/bits/stl_vector.h:351
#10 ~Dimension (this=0x7fff9620, __in_chrg=) at
/home/bates/R/x86_64-unknown-linux-gnu-library/2.15/Rcpp/include/Rcpp/Dimension.h:29
#11 lme4::merPredD::condVar (this=0x572, rho=...) at predModule.cpp:105
#12 0x707dc608 in merPredDcondVar (ptr=0x474c8c0,
rho=) at external.cpp:617
...
(gdb) up 10
#10 ~Dimension (this=0x7fff9620, __in_chrg=) at
/home/bates/R/x86_64-unknown-linux-gnu-library/2.15/Rcpp/include/Rcpp/Dimension.h:29
29  class Dimension {
(gdb) p dims
$1 = { >> = {_M_impl =
{> = {<__gnu_cxx::new_allocator> = {}, }, _M_start = 0x53d75d0,
  _M_finish = 0x53d75dc, _M_end_of_storage = 0x53d75dc}}, }

so the address in question is the _M_start of the dims
std::vector private member of the Rcpp::Dimension instance.  The
length of 12 bytes (difference between _M_finish and _M_start) is what
would be expected for a std::vector of length 3.

Those memory locations have the values we would expect
(gdb) up
#11 lme4::merPredD::condVar (this=0x572, rho=...) at predModule.cpp:105
105 Rcpp::NumericVector ansi(Rcpp::Dimension(ncti, ncti, nli));
(gdb) p ncti
$2 = 1
(gdb) p nli
$3 = 6
(gdb) down
#10 ~Dimension (this=0x7fff9620, __in_chrg=) at
/home/bates/R/x86_64-unknown-linux-gnu-library/2.15/Rcpp/include/Rcpp/Dimension.h:29
29  class Dimension {
(gdb) p {int}0x53d75d0
$4 = 1
(gdb) p {int}0x53d75d4
$5 = 1
(gdb) p {int}0x53d75d8
$6 = 6

To me it looks as if the pointer to the contents is being passed as
the pointer to a malloc'd chunk, which has a header describing the
size of the chunk.  This may be the common approach - I don't know
enough about the innards of malloc and free to tell.

Anyway, I can't see anything suspicious in the code in
Rcpp/src/Dimension.cpp.  For the time being I am going to fall back on
the standard malloc/free

Well, actually I did another thing which is to replace the line 105 by

Rcpp::NumericVector ansi(ncti * ncti * nli);
ansi.attr("dim") = Rcpp::IntegerVector::create(ncti, ncti, nli);

and now things are fine.  Still don't know what the problem was though.
___
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

Re: [Rcpp-devel] Including a compiled function with inline

[Douglas Bates]

On Mon, Jun 4, 2012 at 4:19 AM, Marie Auger-Methe
 wrote:
> Hi list,

> I am writing a Rcpp function (referred as f2). To be able to debug it, I am
> currently using inline to write f2. f2 uses a function that I've wrote and
> is now part of a package that I have made (referred as f1). f1 is an Rcpp
> function that is not called within R but only called from Rcpp functions. I
> was wondering how I could use this f1 function within f2 when using inline.
> From what I understand I should be able to use the argument includes from
> cxxfunction. I think it is related to this thread:
> http://thread.gmane.org/gmane.comp.lang.r.rcpp/2593/focus=2600
> But I can't quite figure out.

Unfortunately, it's complicated if you want to access the compiled
code from another package.   There is a mechanism for another
package's compiled C functions but it is somewhat unwieldy.  If you
want to access another package's compiled C++ classes and methods and
stand-alone functions you need to arrange for the compiled code to be
in a known location and have this location on the path for the loader.
 The Rcpp package does this but getting the details right on multiple
platforms is tricky.

Perhaps the easiest approach is to include not only the headers but
all the C++ code in the file listed in the includes argument to
cxxfunction.  It's a bit wasteful but at least it works. :-)

>
> Here is a small example:
>
> library(inline)
> library(fxwithinfx) # My package that contain f1
>
> src <- '
>   using namespace Rcpp;
>   double y = 9;
>   NumericVector z(1);
>
>   f1(y, z);
>   z[0] = z[0] + 1;
>   return z;
>   '
> f2 <- cxxfunction(signature(),  body = src, plugin = "Rcpp", includes =
> "#include ")
>
> # I get the following error:
>
> Error in compileCode(f, code, language = language, verbose = verbose) :
>   Compilation ERROR, function(s)/method(s) not created!
> file11741d6932c4.cpp:19:27: fatal error: fxwithinfx/f1.h: No such file or
> directory
> compilation terminated.
> # I tried to change the reference to f1.h by changing the path, but whatever
> path I give I still get an error. For example, if I use this instead:
> f2 <- cxxfunction(signature(),  body = src, plugin = "Rcpp",
>  includes = "#include
> ")
> # I get the following error:
>
> Error in compileCode(f, code, language = language, verbose = verbose) :
>   Compilation ERROR, function(s)/method(s) not created!
> file11743a6031d8.o:file11743a6031d8.cpp:(.text+0x104): undefined reference
> to `f1'
> collect2: ld returned 1 exit status
>
> # I'm actually not sure whether I should refer the actual .h location or to
> the location within my R libraries
>
> # Here are the header and  src codes for the f1
>
> f1.h:
>
> #ifndef _fxwithinfx_f1_H
> #define _fxwithinfx_f1_H
>
> #include 
>
> extern "C" int f1(double y, SEXP z) ;
>
> #endif
>
>
> f1.cpp:
>
> #include "f1.h"
>
> int f1(double y, SEXP z){
>   using namespace Rcpp;
>
>   NumericVector zz(z);
>   zz[0] = y + 1.5;
> }
>
>
>
> Any help will be appreciated!
>
> Marie
>
>
> ___
> 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
___
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