[Rcpp-devel] variable affectation in Rcpp Rcpp Armadillo

Nicolas Heslot nh269 at cornell.edu
Fri Jun 17 01:05:56 CEST 2011


Hi dear List,

Another question related to the first one is about the compilation of that
code.
If the source of error is removed, the code compile correctly using
cxxfunction.
However, I get a warning message when I try to produce  a package out of it
using
package.skeleton( "fun2", fun ):
In dump(item, file = file.path(code_dir, sprintf("%s.R", list0[item]))) :
the deparse code  of an S4 object is not always sourceable.

and the compilation on the same machine using R CMD fail

** R
Error in parse(outFile) : D:/fun2/R/fun.R:5:21: '<' unexpected
4:     all_rd_)
5: .Primitive(".Call")(<
                       ^
ERROR: unable to collate files for package 'fun2'

I don't have that problem with the example provided in the documentation but
I don't understand the warning message or how it can be adressed.


Thank you very much for your help

Nicolas

On Thu, Jun 16, 2011 at 1:54 PM, Nicolas  wrote:

> I ran into the same kind of problem again, but this time I just can't find
> what the correct syntax is:
> I need to update a matrix stored in a list, by summing with an arma:: mat
> object
> the syntax MyList[j]  += ArmaObject does not work.
> On the other hand I can access the matrix object stored in those list and
> use them as arma::mat.
> I tried various syntax using wrap or as but I can't get it working.
> The error message is:
> error: no match for 'operator+=' in 'Rcpp::Vector<RTYPE>::operator[](int)
> [with int RTYPE = 19, Rcpp::Vector<RTYPE>::Proxy =
> Rcpp::internal::generic_proxy<19>](k) += arma::operator-(const
> arma::Base<typename T1::elem_type, T1>&, const arma::Base<typename
> T1::elem_type, T2>&) [with T1 = arma::Glue<arma::eOp<arma::subview<double>,
> arma::eop_sc
>
> And the full inline code is:
>
> ## loop code
> src <- '
> // n markers
> int P = as<int>(P_);
> // n individuals
> int N = as<int>(N_);
> // marker polynome degree
> int degM = as<int>(degM_);
> // variance
> //int sigma2 = as<int>(sigma2_);
> //Rcpp::NumericVector sigma2(sigma2_);
> double sigma2 = as<double>(sigma2_);
> // design matrix
> arma::mat gen_a = Rcpp::as<arma::mat>(gen_a_);
> // covariate X
> arma::mat X = Rcpp::as<arma::mat>(X_);
> // alpha
> arma::mat alpha = Rcpp::as<arma::mat>(alpha_);
> // mu
> arma::mat mu = Rcpp::as<arma::mat>(mu_);
> // a_old
> arma::mat a_old = Rcpp::as<arma::mat>(a_old_);
> // tau2
> arma::mat tau2 = Rcpp::as<arma::mat>(tau2_);
> // list
> List all_ui = as<List>(all_ui_);
> List all_ui2 = as<List>(all_ui2_);
> List all_ui3 = as<List>(all_ui3_);
> List all_corMat = as<List>(all_corMat_);
> List all_rd = as<List>(all_rd_);
> arma::mat tmp;
> arma::mat tmp2 = arma::zeros<arma::mat>((degM+1),(degM+1));
> arma::mat tmp3 = arma::zeros<arma::mat>(1,(degM+1));
>
> arma::mat ui;
> arma::mat ui2;
> arma::mat ui3;
> arma::mat sig;
> arma::mat sigma;
> arma::mat sigma3;
> arma::mat rdi;
> arma::mat all_corMati;
> arma::mat diag;
> arma::mat diag2 = arma::zeros<arma::mat>((degM+1),(degM+1));
> arma::mat sol;
> arma::mat aVar_j;
> arma::mat aMu_j;
> int n = 1;
> int p = 0;
> arma::mat U;
> arma::vec s;
> arma::mat V;
> arma::mat rdn;
> arma::mat ans;
> arma::mat vsqrt;
> arma::mat a = arma::zeros<arma::mat>(P,(degM+1));
> List test(N);
>         for (int j = 0; j < P; j++) {
>
>         tmp2.fill(0);
>         tmp3.fill(0);
>               for (int i = 0; i < N; i++) {
>               ui = Rcpp::as<arma::mat>(all_ui[i]);
>               ui2 = Rcpp::as<arma::mat>(all_ui2[i]);
>               ui3 = Rcpp::as<arma::mat>(all_ui3[i]);
>               rdi = Rcpp::as<arma::mat>(all_rd[i]);
>               all_corMati = Rcpp::as<arma::mat>(all_corMat[i]);
>               sigma = sigma2*all_corMati;
>               // calculate mean and variance
>               sigma3 = sigma;
>               diag=diagmat(sigma3.fill(1));
>               tmp = (solve(sigma,diag)* gen_a(i,j))*ui3;
>               tmp2 += (gen_a(i,j)*trans(ui))*tmp;
>               tmp3 += (rdi - mu*trans(ui3)-(X(i)*alpha)*trans(ui2) +
> gen_a(i,j)*a.row(j)*trans(ui))*tmp;
>               };
>         sol = 1/sigma2 *diagmat(diag2.fill(1/tau2(j)))+tmp2;
>         aVar_j= solve(sol,diagmat(diag2.fill(1)) );
>         aMu_j = aVar_j*trans(tmp3);
>         ///////
>         p = aVar_j.n_cols;
>         svd(U,s,V,aVar_j);
>         vsqrt = trans(V*(trans(U)%repmat(trans(sqrt(trans(s))),1,p)));
>         rdn.randn(n,p);
>         ans =rdn*vsqrt + trans(repmat(aMu_j,1,n));
>         a.row(j) = ans;
>               for (int k = 0; k < N; k++) {
>               int k=1;
>               ui = Rcpp::as<arma::mat>(all_ui[k]);
> //////////////////////// Here is the issue!
>               all_rd[k] += gen_a(k,j)*a_old.row(j)*trans(ui)-
> gen_a(k,j)*a.row(j)*trans(ui);
>               };
>         };
> return Rcpp::List::create(wrap(a),wrap(all_rd));'
> ####
> fun <- cxxfunction(signature(P_ = "integer",N_ = "integer",degM_ =
> "integer", sigma2_ = "double", gen_a_="matrix", X_="matrix",alpha_="matrix",
>  mu_="matrix", a_old_="matrix",tau2_="matrix", all_ui_="List",
> all_ui2_="List", all_ui3_="List", all_corMat_="List",
> all_rd_="List"),body=src,plugin="RcppArmadillo")
>
> Thank you very much for your help
>
> Nicolas
>
>
>
>
> I finally figured out, it was a syntax error, as the new assigned object
>> was not an arma:: mat but an element from a list
>> I took care of the problem by using ui = Rcpp::as<arma::mat>(all_ui[i]) in
>> the reaffectation
>>
>> Thank you again all for your help
>>
>>
>> On Wed, Jun 15, 2011 at 11:00 PM, Dirk Eddelbuettel <edd at debian.org>wrote:
>>
>>>
>>> On 15 June 2011 at 20:17, Nicolas Heslot wrote:
>>> | Any idea about this variable reaffectation thing?
>>>
>>> Well ... that still make little sense. Variables are not protected; there
>>> is
>>> nothing stopping you from re-assigning values, be it element-wise, row or
>>> column-wise, or as a whole matrix. So I suspect you have a simple syntax
>>> error.
>>>
>>> But as you still have not shown an actualy error it remains idle
>>> speculation
>>> on our side.
>>>
>>> | It seems to be specific to arma classes like arma::mat
>>>
>>> Not really. We all use these and read values from / assign values into at
>>> will.  There are numerous examples in the archives of this mailing list,
>>> in
>>> the unit tests, and in some of the slides from past presentations.
>>>
>>> Dirk
>>>
>>> --
>>> Gauss once played himself in a zero-sum game and won $50.
>>>                      -- #11 at http://www.gaussfacts.com
>>>
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20110616/398f38e1/attachment-0001.htm>


More information about the Rcpp-devel mailing list