[Rcpp-devel] Variable for loop

Amina Shahzadi aminashahzadi at gmail.com
Tue Sep 27 06:09:34 CEST 2016


Waoooo. Impressive----I have tried it. It is giving the same answer I want.
Now I learn it step by step how you did that. It is very different concept
from what I was using in rcpp.
Thank you so very much Avraham.

On Tue, Sep 27, 2016 at 4:52 PM, Avraham Adler <avraham.adler at gmail.com>
wrote:

> Hello, Amina.
>
> It is much, much, much easier to deal with small, simplified,
> examples, as Dirk and others have explained before. Changing your code
> to make n = 4 and the values seq(1, 12), we can see the pattern. For
> each slice, you are taking the product down the rows of "slice" number
> of entries. So the first slice is the 1-element product of the
> entries, or just a copy of the original matrix. Now, putting the index
> names back in familiar context (M rows and N columns), for the second
> slice, C(2, 1, 2) = A(1, 1) * A(2, 1). C(3, 1, 2) = A(2, 1) * A(3, 1),
> etc. For the third slice, we have 3-element products So C(1, , 3) and
> C(2, , 3) are 0. C(3, 1, 3) = A(1, 1) * A(2, 1) * A(3, 1), etc. In the
> last slice there is only one entry in the last row and it is the
> product of each column in its entirety.
>
> That being said, the way I would approach it is to use Armadillo's
> subsetting to get the needed vector and use Armadillo's product
> function to get their product.
>
> First, we'll turn your fixed R code into a function which we can reuse
> by creating any a you want and passing it to up_R
>
> up_R <- function(a){
>   m = dim(a)[[1]]
>   n = dim(a)[[2]]
>   b = array(as.double(0), dim = c(m, n, m))
>   for(i in 1:m){
>     for(j in 1:n){
>       for(q in 1:i){
>         if(q==1){
>           b[i, j, q] = a[i, j];
>         }
>         else{
>           b[i, j, q] = b[i-1, j, q-1] * a[i, j];
>         }
>       }
>     }
>   }
>   return(b)
> }
>
> The way I would approach the problem from first principles, knowing
> what you actually intended, is as follows:
>
> #include <RcppArmadillo.h>
> using namespace Rcpp;
> using namespace RcppArmadillo;
> using namespace arma;
> //[[Rcpp::depends(RcppArmadillo)]]
>
> //[[Rcpp::export]]
> cube up_C(mat a){
>   int m = a.n_rows;
>   int n = a.n_cols;
>   cube C = cube(m, n, m, fill::zeros);
>   C.slice(0) = a;
>   for (int slices = 1; slices < m; ++slices){
>     for (int j = 0; j < n; ++j) {
>       for (int i = slices; i < m; ++i) {
>         vec S = C.subcube(i - slices, j, 0, i, j, 0);
>         C(i, j, slices) = prod(S);
>       }
>     }
>   }
>   return(C);
> }
>
> This pulls the "slices" length subset of the jth column starting at
> slices and ending at i from the 0th slice (which is the original
> matrix) and puts the product in to the C(i, j, slices) entry.
>
> You can test the two are identical through identical(up_C(a), up_R(a))
> for any a, and the C version is, of course, much faster especially for
> large matrices.
>
> It often pays to start with something simpler, and see if there is a
> pattern you can identify which would help you conceptually, even if
> not completely optimized.
>
> Avi
>
>
>
>
> On Mon, Sep 26, 2016 at 10:31 PM, Amina Shahzadi
> <aminashahzadi at gmail.com> wrote:
> > Hi Dear
> > This is the corresponding R code which I want to have the same result
> from
> > rcpp code.
> >
> > set.seed(11)
> > a = matrix(runif(30), nrow=10, ncol=3)
> > n=10
> > m=3
> > b = array(as.double(0), dim = c(n, m, n))
> > for(i in 1:n){
> >   for(j in 1:m){
> >     for(q in 1:i){
> >       if(q==1){
> >         b[i, j, q] = a[i, j];
> >       }
> >       else{
> >         b[i, j, q] = b[i-1, j, q-1] * a[i, j];
> >       }
> >     }
> >   }
> > }
> >
> >
> > On Tue, Sep 27, 2016 at 3:22 PM, Avraham Adler <avraham.adler at gmail.com>
> > wrote:
> >>
> >> Please give a simple example of an input and its expected output.
> >> Unfortunately "run of a variable for loop" is too general for me to
> >> understand.
> >>
> >> Avi
> >>
> >> On Mon, Sep 26, 2016 at 10:07 PM, Amina Shahzadi
> >> <aminashahzadi at gmail.com> wrote:
> >> > Hi Dear
> >> >
> >> > My purpose is to make run of a variable for loop. Here I have assumed
> >> > else
> >> > statement to be zero. Otherwise it could be anything for example in
> the
> >> > following code.
> >> > In this code, my main aim to run the slices loop according to rows.
> >> > does it now look ok?
> >> >
> >> >
> >> > Thank you
> >> >
> >> > #include <RcppArmadillo.h>
> >> > using namespace Rcpp;
> >> > using namespace RcppArmadillo;
> >> > //[[Rcpp::depends(RcppArmadillo)]]
> >> > //[[Rcpp::export]]
> >> > arma::cube up(arma::mat a){
> >> >   int m = a.n_cols;
> >> >   int n = a.n_rows;
> >> >   int p = a.n_rows;
> >> >   arma::cube b(n, m, p);
> >> >   for(int i=0; i<n; i++){
> >> >       for(int j=0; j<m; j++){
> >> >         for(int q=0; q<i; q++){
> >> >          if(q==0){
> >> >            b(i, j, q) = a(i, j);
> >> >          }
> >> >          else{
> >> >            b(i, j, q) = b(i-1, j, q-1) * a(i, j);
> >> >          }
> >> >       }
> >> >     }
> >> >   }
> >> >   return b;
> >> > }
> >> >
> >> > On Tue, Sep 27, 2016 at 2:51 PM, Avraham Adler <
> avraham.adler at gmail.com>
> >> > wrote:
> >> >>
> >> >> Hello.
> >> >>
> >> >> Once again, it is very unclear what you want to do. Can you please
> >> >> explain, in English not code what your procedure intends to do, the
> >> >> input you expect, and the output you expect?
> >> >>
> >> >> What it LOOKS like you want to do is to create an N x M x N cube
> where
> >> >> the first slice is your matrix and the remaining slices are all 0. If
> >> >> that is the case, There are much, much simpler ways to do it than to
> >> >> traverse all N²M cells. The following should work.
> >> >>
> >> >> #include <RcppArmadillo.h>
> >> >> using namespace Rcpp;
> >> >> using namespace RcppArmadillo;
> >> >> //[[Rcpp::depends(RcppArmadillo)]]
> >> >>
> >> >> //[[Rcpp::export]]
> >> >> arma::cube fillup(arma::mat a){
> >> >>   int m = a.n_cols;
> >> >>   int n = a.n_rows;
> >> >>   arma::cube C = arma::cube(n, m, n, arma::fill::zeros);
> >> >>
> >> >>   C.slice(0) = a;
> >> >>
> >> >>   return(C);
> >> >> }
> >> >>
> >> >> Avi
> >> >>
> >> >> On Mon, Sep 26, 2016 at 5:59 PM, Amina Shahzadi
> >> >> <aminashahzadi at gmail.com>
> >> >> wrote:
> >> >> > Hi Dear
> >> >> >
> >> >> > I have a problem in using a variable for loop in using
> RcppArmadillo
> >> >> > library.
> >> >> > I have pasting here my code. It is executing but not giving the
> same
> >> >> > results
> >> >> > as its R code version gives. The results produced by it are really
> >> >> > weird. I
> >> >> > have checked it step by step. It is because of the for (int q=0;
> q<i;
> >> >> > q++).
> >> >> > I request tp please help how to handle it in cpp.
> >> >> >
> >> >> > The another question is I want to multiply the cube b(i, ,) by a
> >> >> > scalar.
> >> >> > How
> >> >> > can we consider the entire columns and slices of a cube for each of
> >> >> > the
> >> >> > rows.  "b(span(i), span(), span())"  is not working for me.
> >> >> >
> >> >> > Thank you
> >> >> >
> >> >> > #include <RcppArmadillo.h>
> >> >> > using namespace Rcpp;
> >> >> > using namespace RcppArmadillo;
> >> >> > //[[Rcpp::depends(RcppArmadillo)]]
> >> >> > //[[Rcpp::export]]
> >> >> > arma::cube up(arma::mat a){
> >> >> >   int m = a.n_cols;
> >> >> >   int n = a.n_rows;
> >> >> >   int p = a.n_rows;
> >> >> >   arma::cube b(n, m, p);
> >> >> >   for(int i=0; i<n; i++){
> >> >> >       for(int j=0; j<m; j++){
> >> >> >         for(int q=0; q<i; q++){
> >> >> >          if(q==0){
> >> >> >            b(i, j, q) = a(i, j);
> >> >> >          }
> >> >> >          else{
> >> >> >            b(i, j, q) = 0.0;
> >> >> >          }
> >> >> >       }
> >> >> >     }
> >> >> >   }
> >> >> >   return b;
> >> >> > }
> >> >> >
> >> >> > --
> >> >> > Amina Shahzadi
> >> >
> >> >
> >> >
> >> >
> >> > --
> >> > Amina Shahzadi
> >
> >
> >
> >
> > --
> > Amina Shahzadi
>



-- 
*Amina Shahzadi*
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20160927/3c724a51/attachment-0001.html>


More information about the Rcpp-devel mailing list