<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"><head><META HTTP-EQUIV="Content-Type" CONTENT="text/html; charset=us-ascii"><meta name=Generator content="Microsoft Word 14 (filtered medium)"><style><!--
/* Font Definitions */
@font-face
{font-family:"Cambria Math";
panose-1:2 4 5 3 5 4 6 3 2 4;}
@font-face
{font-family:Calibri;
panose-1:2 15 5 2 2 2 4 3 2 4;}
@font-face
{font-family:Consolas;
panose-1:2 11 6 9 2 2 4 3 2 4;}
/* Style Definitions */
p.MsoNormal, li.MsoNormal, div.MsoNormal
{margin:0cm;
margin-bottom:.0001pt;
font-size:11.0pt;
font-family:"Calibri","sans-serif";
mso-fareast-language:EN-US;}
a:link, span.MsoHyperlink
{mso-style-priority:99;
color:blue;
text-decoration:underline;}
a:visited, span.MsoHyperlinkFollowed
{mso-style-priority:99;
color:purple;
text-decoration:underline;}
pre
{mso-style-priority:99;
mso-style-link:"HTML Preformatted Char";
margin:0cm;
margin-bottom:.0001pt;
font-size:10.0pt;
font-family:"Courier New";}
span.EmailStyle17
{mso-style-type:personal-compose;
font-family:"Calibri","sans-serif";
color:windowtext;}
span.HTMLPreformattedChar
{mso-style-name:"HTML Preformatted Char";
mso-style-priority:99;
mso-style-link:"HTML Preformatted";
font-family:"Courier New";
mso-fareast-language:EN-GB;}
span.gjwpqfqdc4
{mso-style-name:gjwpqfqdc4;}
span.gjwpqfqdo4
{mso-style-name:gjwpqfqdo4;}
.MsoChpDefault
{mso-style-type:export-only;
mso-fareast-language:EN-US;}
@page WordSection1
{size:612.0pt 792.0pt;
margin:72.0pt 72.0pt 72.0pt 72.0pt;}
div.WordSection1
{page:WordSection1;}
--></style><!--[if gte mso 9]><xml>
<o:shapedefaults v:ext="edit" spidmax="1026" />
</xml><![endif]--><!--[if gte mso 9]><xml>
<o:shapelayout v:ext="edit">
<o:idmap v:ext="edit" data="1" />
</o:shapelayout></xml><![endif]--></head><body lang=EN-GB link=blue vlink=purple><div class=WordSection1><p class=MsoNormal><span style='font-size:12.0pt'>Hi all<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>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.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>In the code below it calculates rolling means. With a simple test data set it works as I want e.g.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:Consolas;color:blue;mso-fareast-language:EN-GB'>x <- 1:20<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:Consolas;color:blue;mso-fareast-language:EN-GB'>fun(x, 8, 75)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:10.0pt;font-family:Consolas;color:black;mso-fareast-language:EN-GB'> [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<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>And if I run it again I get the same result.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>Now if I divide x by 2 and run it again I get:<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><pre><span class=gjwpqfqdc4><span style='font-family:Consolas;color:blue'>x = x/2<o:p></o:p></span></span></pre><pre><span class=gjwpqfqdo4><span style='font-family:Consolas;color:blue'>> </span></span><span class=gjwpqfqdc4><span style='font-family:Consolas;color:blue'>fun(x, 8, 75)<o:p></o:p></span></span></pre><pre><span style='font-family:Consolas;color:black'> [1] NA NA NA NA NA NA NA 2.250000 2.531250 2.785156 3.008301<o:p></o:p></span></pre><pre><span style='font-family:Consolas;color:black'>[12] 3.196838 3.346443 3.452249 3.508780 3.728627 3.940799 4.147755 4.352686 4.559667<o:p></o:p></span></pre><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>Which is also correct. BUT if I run the fun(x, 8, 75) again I get all NAs:<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><pre><span class=gjwpqfqdc4><span style='font-family:Consolas;color:blue'>fun(x, 8, 75)<o:p></o:p></span></span></pre><pre><span style='font-family:Consolas;color:black'> [1] NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA<o:p></o:p></span></pre><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>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.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>Thanks<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>David<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>library(inline)<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>src <- '<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>Rcpp::<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>NumericVector A(x);<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>double capr = as<double>(cap); // data capture %<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>int lenr = as<int>(len); // window size %<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>NumericVector res(x); // for results<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>LogicalVector NA(x);<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>NumericVector missing(1);<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>missing[0] = NA_REAL;<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>int n = A.size();<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>double sum = 0.0;<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>double sumNA = 0; // number of missings<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>NA = is_na(A) ; // logical vector of missings<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>for (int i = 0; i <= (n - lenr); i++) {<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> sum = 0; // initialise<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> sumNA = 0;<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> for (int j = i; j < i + lenr; j++) {<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> //sum += A(j); // increment sum.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> if (NA[j]) {<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> sumNA += 1;<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> }<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> else<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> {<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> sum += A[j];<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> }<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> }<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> // Calculate moving average and display.<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> if (1 - sumNA / lenr < capr / 100) {<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> res(i + lenr - 1) = missing[0];<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> }<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> else<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> {<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> res(i + lenr - 1) = sum / (lenr - sumNA);<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> }<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> }<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>// pad out missing data<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>for (int i = 0; i < lenr - 1; i++) {<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> res(i) = missing[0];<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> }<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>return res;<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'> '<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'>fun <- cxxfunction(signature(x = "numeric", len = "int", cap = "double"), body = src, plugin="Rcpp")<o:p></o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='font-size:12.0pt'><o:p> </o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'><o:p> </o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>David Carslaw<o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'><o:p> </o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>Science Policy Group<o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>Environmental Research Group<o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>MRC-HPA Centre for Environment and Health <o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>King's College London <o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>Room 4.129 Franklin Wilkins Building <o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>Stamford Street <o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>London SE1 9NH <o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'><a href="mailto:david.carslaw@kcl.ac.uk"><span style='color:blue'>david.carslaw@kcl.ac.uk</span></a><o:p></o:p></span></p><p class=MsoNormal><span style='mso-fareast-language:EN-GB'>T 020 7848 3844 <o:p></o:p></span></p><p class=MsoNormal><o:p> </o:p></p></div></body></html>