Migrate to Rcpp

Author: hafen

DESCRIPTION

@@ -14,9 +14,11 @@ NeedsCompilation: no
 Imports:
     lattice,
     yaImpute,
-    stats
+    stats,
+    Rcpp
 Suggests:
     testthat
+LinkingTo: Rcpp
 License: BSD_3_clause + file LICENSE
 URL: https://github.com/stlplus
 Note: This is experimental software distributed free of charge and

NAMESPACE

@@ -26,6 +26,7 @@ export(remainder)
 export(seasonal)
 export(stlplus)
 export(trend)
+importFrom(Rcpp,sourceCpp)
 importFrom(lattice,panel.abline)
 importFrom(lattice,panel.loess)
 importFrom(lattice,panel.segments)

NEWS.md

@@ -2,3 +2,4 @@ Version 0.5.0
 ----------------------------------------------------------------------
 
 - Major overhaul of code to match more modern R package development standards
+- Update old C interface to Rcpp

R/RcppExports.R

@@ -0,0 +1,15 @@
+# This file was generated by Rcpp::compileAttributes
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+c_interp <- function(m, fits, slopes, at) {
+    .Call('stlplus_c_interp', PACKAGE = 'stlplus', m, fits, slopes, at)
+}
+
+c_loess <- function(xx, yy, degree, span, ww, m, l_idx, max_dist) {
+    .Call('stlplus_c_loess', PACKAGE = 'stlplus', xx, yy, degree, span, ww, m, l_idx, max_dist)
+}
+
+c_ma <- function(x, n_p) {
+    .Call('stlplus_c_ma', PACKAGE = 'stlplus', x, n_p)
+}
+

R/cpp_wrappers.R

@@ -0,0 +1,8 @@
+.interp <- function(m, fits, slopes, at) {
+  if(any(is.nan(fits))) {
+    ind <- !is.nan(fits)
+    c_interp(m[ind], fits[ind], slopes[ind], at)
+  } else {
+    c_interp(m, fits, slopes, at)
+  }
+}

R/loess_stl.R

@@ -1,4 +1,6 @@
 #' @importFrom yaImpute ann
+#' @importFrom Rcpp sourceCpp
+#' @useDynLib stlplus
 .loess_stlplus <- function(x = NULL, y, span, degree, weights = NULL,
   m = c(1:length(y)), y_idx = !is.na(y), noNA = all(y_idx), blend = 0,
   jump = ceiling(span / 10), at = c(1:length(y))) {
@@ -53,16 +55,13 @@
     # max_dist <- max_dist * (span / n)
     max_dist <- max_dist + (span - n) / 2
 
-  out <- .C("loess_stlplus", as.double(x[y_idx]), as.double(y[y_idx]),
-    as.integer(n), as.integer(degree), as.integer(span),
-    as.double(weights[y_idx]), as.integer(m), as.integer(n_m),
-    as.integer(l_idx - 1), as.integer(r_idx - 1), as.double(max_dist),
-    result = double(n_m), slopes = double(n_m), PACKAGE = "stlplus")
+  out <- c_loess(x[y_idx], y[y_idx], degree, span, weights[y_idx],
+    m, l_idx - 1, as.double(max_dist))
 
   res1 <- out$result
   # do interpolation
   if(jump > 1)
-    res1 <- .stlplus_interp(m, out$result, out$slope, at)
+    res1 <- .interp(m, out$result, out$slope, at)
     # res1 <- approx(x = m, y = out$result, xout = at)$y
 
   if(blend > 0 && blend <= 1 && degree >= 1) {
@@ -98,17 +97,14 @@
     # right now, a lot of unnecessary calculation is done at the interior
     # where blending doesn't matter
 
-    tmp <- .C("loess_stlplus", as.double(x[y_idx]), as.double(y[y_idx]),
-      as.integer(n), as.integer(0), as.integer(sp0), as.double(weights[y_idx]),
-      as.integer(m2), as.integer(n_m2), as.integer(l_idx2-1),
-      as.integer(r_idx2-1),  as.double(max_dist2),
-      result = double(n_m2), slopes = double(n_m2), PACKAGE = "stlplus")
+    tmp <- c_loess(x[y_idx], y[y_idx], 0, sp0, weights[y_idx],
+      m2, l_idx2-1, max_dist2)
 
     if(jump > 1) {
-      res2_left <- .stlplus_interp(left,
+      res2_left <- .interp(left,
         head(tmp$result, length(left)),
         head(tmp$slope, length(left)), left_interp)
-      res2_right <- .stlplus_interp(right,
+      res2_right <- .interp(right,
         tail(tmp$result, length(right)),
         tail(tmp$slope, length(right)), right_interp)
     } else {

R/ma.R

@@ -43,7 +43,6 @@
 #' @example man-roxygen/ex-stlplus.R
 #' @export
 #' @rdname stlplus
-#' @useDynLib stlplus
 stlplus <- function(x, t = NULL, n.p, s.window, s.degree = 1,
   t.window = NULL, t.degree = 1,
   fc.window = NULL, fc.degree = NULL, fc.name = NULL,
@@ -301,14 +300,14 @@ stlplus.default <- function(x, t = NULL, n.p, s.window, s.degree = 1,
           tmps <- .loess_stlplus(y = cycleSub, span = s.window, degree = s.degree,
             m = cs.ev, weights = w[cycleSubIndices == i], blend = s.blend,
             jump = s.jump,  at = c(0:(cycleSub.length + 1)))
-          C[c(cs1, cycleSubIndices, cs2)==i] <- tmps
+          C[c(cs1, cycleSubIndices, cs2) == i] <- tmps
           # approx(x = cs.ev, y = tmps, xout = c(0:(cycleSub.length + 1)))$y
         }
       }
 
       # Step 3: Low-pass filtering of collection of all the cycle-subseries
       # moving averages
-      ma3 <- .ma(C, n.p)
+      ma3 <- c_ma(C, n.p)
 
       l.ev <- seq(1, n, by = l.jump)
       if(tail(l.ev, 1) != n) l.ev <- c(l.ev, n)