[Rcpp-devel] strange code problem

[Steve Lianoglou]

Hi David,

Interesting problem.

You'll find that if you change how you initialize your `res` vector,
you will sidestep this problem, eg. from this:

    NumericVector res(x);

to:

    NumericVector res(A.size());

Why? Because in your original instantiation, as you "save" your
results in `res`, you are trampling over your input.

Here's a short example:

R> s2 = "
  NumericVector x(x_);
  x[0] = 100;
  return x;
"

R> f <- cxxfunction(signature(x_ = "numeric"), body = s2, plugin="Rcpp")
R> x <- c(0.5, 1.0, 2.0)
R> f(x)
## [1] 100   1   2

## But look, we also changed x!
R> x
## [1] 100   1   2

## Now try again with this `x`, and ask yourself why the outcome is different
R> x <- 1:5
R> f(x)
## [1] 100   2   3   4   5

R> x
[1] 1 2 3 4 5

Welcome to the wild west.

HTH,
-steve



On Wed, Feb 15, 2012 at 1:09 PM, Carslaw, David <david.carslaw at kcl.ac.uk> 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



-- 
Steve Lianoglou
Graduate Student: Computational Systems Biology
 | Memorial Sloan-Kettering Cancer Center
 | Weill Medical College of Cornell University
Contact Info: http://cbio.mskcc.org/~lianos/contact