[Rcpp-devel] Working with indicators in Armadillo

[Conrad S]

On Sun, Feb 24, 2013 at 6:33 AM, Simon Zehnder <szehnder at uni-bonn.de> wrote:
>   1 arma::mat repS = arma::repmat(S, 1, K);
>   2 arma::mat compM = arma::ones(S.n_rows, K);
>   3 arma::rowvec par_post(K);
>   4
>   5 for(unsigned int k = 0; k < K; ++k) {
>   6               compM.col(k) = compM.col(k) * (k + 1);
>   7 }
>   8
>   9 arma::umat ind = (repS == compM);
> 10 par_post = sum(ind);

There is a performance issue on line 6.  I'll get back to this.

First, the source of the compilation error is a type mismatch.  The
result of sum(ind) is a "umat" (with one row) while par_post has the
type "rowvec".

"umat" is a shorthand (typedef) for Mat<uword>, where uword is an
unsigned integer type.
"rowvec" is a typedef for Row<double>, which in effect is Mat<double>
(albeit with one row).

So the underlying element type for the two Mat objects is different:
one is uword, and the other is double. As such, you have a type
mismatch.

You can use the conv_to() function to convert the result of sum() to
be compatible with rowvec, eg.
par_post = conv_to<rowvec>::from(sum(ind)).
Alternatively, you can multiply sum(ind) by 1.0 which should then
convert it to a compatible time.

About the performance issue on line 6.  Your code will run much
quicker if you change it from:
   compM.col(k) = compM.col(k) * (k + 1);
to:
   compM.col(k) *= (k + 1);

In the former, Armadillo aliasing detector will indicate that one of
the inputs is an alias of the output, and that one of the inputs is a
submatrix.  When using submatrices, detecting safe/non-destructive
aliasing is hard to do (difficult and tedious to code) so Armadillo
will currently make a copy of compM.col(k) before assigning the result
to compM.col(k).  You can avoid this by using the *= operator.