[Rcppdevel] strange code problem
Dirk Eddelbuettel
edd at debian.org
Wed Feb 15 19:50:49 CET 2012
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 lenr1
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

 MRCHPA 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




 
 _______________________________________________
 Rcppdevel mailing list
 Rcppdevel at lists.rforge.rproject.org
 https://lists.rforge.rproject.org/cgibin/mailman/listinfo/rcppdevel

"Outside of a dog, a book is a man's best friend. Inside of a dog, it is too
dark to read."  Groucho Marx
More information about the Rcppdevel
mailing list