[Rcpp-devel] strange code problem

[Dirk Eddelbuettel]

Hi David,

I didn't have time to study this in detail but it smells like a tricky issue
one which has been discussed here.  It has to do with whether

  -- your argument vector is integer (x <- 1:20 or even x <- 1L:20L) or
     numeric (seq(1.0, 20.0., by=1.0)

  -- because you use a NumericVector A(x)

  -- AND because you init. results to be the same vector   res(x)   as
     opposed to res(lnr)   where lenr is a scalar --> new vector rather
     proxy object to x

  -- I think you overwrite too much in the final loop from 0 to lenr-1

If we run your example as 

  x <- seq(1.0, 20.0, by=1.0)
  print(fun(x, 8, 75))
  print(fun(x, 8, 75))

  x <- x/2
  print(fun(x, 8, 75))
  print(fun(x, 8, 75))

it turns to NA much sooner:

edd at max:~$ r /tmp/david.r 
Loading required package: methods
 [1]       NA       NA       NA       NA       NA       NA       NA 4.500000
 [9] 5.062500 5.570312 6.016602 6.393677 6.692886 6.904497 7.017559 7.457254
[17] 7.881598 8.295509 8.705373 9.119335
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
 [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
edd at max:~$ 

I didn't change your program, other than by indenting.  It is included below.

Dirk


library(inline)

src <- '
   Rcpp::NumericVector A(x);
   double capr = as<double>(cap); // data capture %
   int lenr = as<int>(len); // window size %
   NumericVector res(x); // for results
   LogicalVector NA(x);
   NumericVector missing(1);
   missing[0] = NA_REAL;
   int n = A.size();
   double sum = 0.0;
   double sumNA = 0; // number of missings
   NA = is_na(A) ; // logical vector of missings
   for (int i = 0; i <= (n - lenr); i++) {
     sum = 0; // initialise
     sumNA = 0;
     for (int j = i; j < i + lenr; j++) {
       //sum += A(j); // increment sum.
       if (NA[j]) {
         sumNA += 1;
       } else {
         sum += A[j];
       }
     }
     // Calculate moving average and display.
     if (1 - sumNA / lenr < capr / 100) {
       res(i + lenr - 1) = missing[0];
     } else {
       res(i + lenr - 1) = sum / (lenr - sumNA);
     }
   }
   // pad out missing data
   for (int i = 0; i < lenr - 1; i++) {
     res(i) = missing[0];
   }
   return res;
'

fun <- cxxfunction(signature(x = "numeric", len = "int", cap = "double"), body=src, plugin="Rcpp")

## x <- 1:20
## print(fun(x, 8, 75))
## print(fun(x, 8, 75))

## x <- x/2
## print(fun(x, 8, 75))
## print(fun(x, 8, 75))

#x <- seq(1.0, 20.0, by=1.0)
x <- 1L:20L
print(fun(x, 8, 75))
print(fun(x, 8, 75))

x <- x/2
print(fun(x, 8, 75))
print(fun(x, 8, 75))




On 15 February 2012 at 18:09, Carslaw, David wrote:
| Hi all
| 
|  
| 
| I?m new to Rcpp ? and C++, but I am keen to use it to speed up simple parts of
| my R package and to learn more.  I thought I had made a good start, but I have
| encountered a strange problem.  Now, this may well be due to my poor knowledge
| of C++, but I can?t find where the problem is and I?d appreciate some help.
| 
|  
| 
| In the code below it calculates rolling means. With a simple test data set it
| works as I want e.g.
| 
|  
| 
| x <- 1:20
| 
| fun(x, 8, 75)
| 
| [1]   NA   NA   NA   NA   NA   NA   NA  4.5  5.5  6.5  7.5  8.5  9.5 10.5 11.5
| 12.5 13.5 14.5 15.5 16.5
| 
|  
| 
| And if I run it again I get the same result.
| 
|  
| 
| Now if I divide x by 2 and run it again I get:
| 
|  
| 
| x = x/2
| 
| > fun(x, 8, 75)
| 
|  [1]       NA       NA       NA       NA       NA       NA       NA 2.250000 2.531250 2.785156 3.008301
| 
| [12] 3.196838 3.346443 3.452249 3.508780 3.728627 3.940799 4.147755 4.352686 4.559667
| 
|  
| 
| Which is also correct.  BUT if I run the fun(x, 8, 75) again I get all NAs:
| 
|  
| 
| fun(x, 8, 75)
| 
|  [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
| 
|  
| 
| Now I appreciate that I maybe I?m doing something silly ? but can you tell me
| what?!  This is on Windows 7, Rcpp_0.9.9.  Any help, much appreciated.
| 
|  
| 
| Thanks
| 
|  
| 
| David
| 
|  
| 
|  
| 
| library(inline)
| 
| src <- '
| 
| Rcpp::
| 
| NumericVector A(x);
| 
| double capr = as<double>(cap); // data capture %
| 
| int lenr = as<int>(len); // window size %
| 
| NumericVector res(x); // for results
| 
| LogicalVector NA(x);
| 
| NumericVector missing(1);
| 
| missing[0] = NA_REAL;
| 
| int n = A.size();
| 
| double sum = 0.0;
| 
| double sumNA = 0; // number of missings
| 
| NA = is_na(A) ; // logical vector of missings
| 
| for (int i = 0; i <= (n - lenr); i++) {
| 
|   sum = 0; // initialise
| 
|   sumNA = 0;
| 
|   for (int j = i; j < i + lenr; j++) {
| 
|     //sum += A(j); // increment sum.
| 
|     if (NA[j]) {
| 
|       sumNA += 1;
| 
|     }
| 
|     else
| 
|       {
| 
|         sum += A[j];
| 
|       }
| 
|   }
| 
|  
| 
|   // Calculate moving average and display.
| 
|  
| 
|   if (1 - sumNA / lenr < capr / 100) {
| 
|     res(i + lenr - 1) = missing[0];
| 
|   }
| 
|   else
| 
|     {
| 
|       res(i + lenr - 1) = sum / (lenr - sumNA);
| 
|     }
| 
|  
| 
| }
| 
|  
| 
| // pad out missing data
| 
| for (int i = 0; i < lenr - 1; i++) {
| 
|   res(i) = missing[0];
| 
| }
| 
| return res;
| 
| '
| 
| fun <- cxxfunction(signature(x = "numeric", len = "int", cap = "double"), body
| = src, plugin="Rcpp")
| 
|  
| 
|  
| 
|  
| 
| David Carslaw
| 
|  
| 
| Science Policy Group
| 
| Environmental Research Group
| 
| MRC-HPA Centre for Environment and Health
| 
| King's College London
| 
| Room 4.129 Franklin Wilkins Building
| 
| Stamford Street
| 
| London SE1 9NH
| 
| david.carslaw at kcl.ac.uk
| 
| T 020 7848 3844
| 
|  
| 
| 
| ----------------------------------------------------------------------
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
-- 
"Outside of a dog, a book is a man's best friend. Inside of a dog, it is too
dark to read." -- Groucho Marx