[Rcpp-devel] Variable for loop

[Avraham Adler]

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