Start Rcpp implementation

Implemented calEmis in with Rcpp

.Rbuildignore

@@ -1,3 +1,5 @@
 ^make-Rda$
 ^schematic.png$
 ^LICENSE$
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

ALFAM2.Rproj

@@ -0,0 +1,17 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: knitr
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source

DESCRIPTION

@@ -1,13 +1,18 @@
 Package: ALFAM2
 Type: Package
 Title: Model and Data on Ammonia Emission from Field-Applied Manure
-Version: 2.1
-Date: 2022-03-07
-Author: Sasha D. Hafner, Christoph Haeni
-Maintainer: Sasha D. Hafner <sasha.hafner@eng.au.dk>
+Version: 2.1.1
+Date: 2022-06-03
+Authors@R: c(person("Sasha D.", "Hafner", role = c("aut", "cre"),
+                     email = "sasha.hafner@eng.au.dk"),
+              person("Christoph", "Haeni", role = "aut"),
+              person("Roland", "Fuss", role = "aut",
+                     email = "roland.fuss@thuenen.de"))
 Depends: R (>= 3.5.0)
+Imports: Rcpp (>= 1.0.8), stats, parallel
 Suggests: testthat, knitr
 Description: An implementation of the ALFAM2 model for ammonia emission from field-applied manure. Model can be used with default parameters or applied to specific cases (locations, application methods, etc.).
 License: GPL-3
 VignetteBuilder: knitr
 LazyData: true
+LinkingTo: Rcpp

NAMESPACE

@@ -1,3 +1,5 @@
 import(parallel)
+import(Rcpp)
 importFrom("stats", "na.omit", "setNames")
 export(alfam2, ALFAM2mod)
+useDynLib(ALFAM2)

R/RcppExports.R

@@ -0,0 +1,7 @@
+# Generated by using Rcpp::compileAttributes() -> do not edit by hand
+# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+rcpp_calcEmis <- function(ct, a0, u0, r1, r2, r3, f4) {
+    .Call('_ALFAM2_rcpp_calcEmis', PACKAGE = 'ALFAM2', ct, a0, u0, r1, r2, r3, f4)
+}
+

R/calcEmis.R

@@ -5,51 +5,15 @@ calcEmis <- function(ct, a0, u0, r1, r2, r3, f4, drop.rows) {
     stop('In calcEmis, ct is not sorted.')
   }
 
-  # Number of intervals
-  l <- length(ct)
-
-  # Empty result vectors
-  a <- u <- e <- numeric(l)
-
-  # Time step
-  ddt <- diff(c(0, ct))
-
-  # Initial values
-  ati0 <- a0
-  uti0 <- u0
-  eti <- 0
-
-  # Then calculate pools at end of ddt[1]
-  # ti = at transfer interval start
-  # Loop through the periods
-  for(i in 1:l) {
-
-    # Make incorporation transfer (at *start* of interval) (if none then f4 = 1 and ati = a[i])
-    ati <- f4[i] * ati0
-    uti <- (1 - f4[i]) * ati0 + uti0
-
-    # Calculate pools at *end* of ct[i]
-    a[i] <- ati*exp(-(r1[i] + r2[i])*ddt[i])
-    u[i] <- exp(-r3[i]*ddt[i]) * (r2[i] * ati * (exp((-r1[i]-r2[i]+r3[i]) * ddt[i]) - 1)/(-r1[i]-r2[i]+r3[i]) + uti)
-    e[i] <- eti + (uti - u[i]) +  (ati - a[i])
-
-    # save pools for next step
-    ati0 <- a[i]
-    uti0 <- u[i]
-    eti <- e[i]
-
-  }
-
-  # Now combine and drop rows that were added for incorporation in afMod()
-  out <- data.frame(ct = ct, dt = NA, 
-                    f0 = a0/(u0 + a0), r1 = r1, r2 = r2, r3 = r3, f4 = f4, 
-                    f = a, s = u, 
-                    j = NA, e = e, e.int = NA, er = e/(a0 + u0))
+  # call Rcpp function
+  out <- as.data.frame(rcpp_calcEmis(ct, a0, 
+                                     u0, r1, 
+                                     r2, r3, 
+                                     f4))
 
   out <- out[!drop.rows, ]
 
   # Recalculate ddt, e.int, and j.ave.int in case there were dropped rows
-  out$dt <- diff(c(0, out$ct))
   out$e.int <- diff(c(0, out$e))
   out$j <- out$e.int/out$dt
 

src/RcppExports.cpp

@@ -0,0 +1,39 @@
+// Generated by using Rcpp::compileAttributes() -> do not edit by hand
+// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
+
+#include <Rcpp.h>
+
+using namespace Rcpp;
+
+#ifdef RCPP_USE_GLOBAL_ROSTREAM
+Rcpp::Rostream<true>&  Rcpp::Rcout = Rcpp::Rcpp_cout_get();
+Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
+#endif
+
+// rcpp_calcEmis
+List rcpp_calcEmis(const NumericVector ct, const double a0, const double u0, const NumericVector r1, const NumericVector r2, const NumericVector r3, const NumericVector f4);
+RcppExport SEXP _ALFAM2_rcpp_calcEmis(SEXP ctSEXP, SEXP a0SEXP, SEXP u0SEXP, SEXP r1SEXP, SEXP r2SEXP, SEXP r3SEXP, SEXP f4SEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< const NumericVector >::type ct(ctSEXP);
+    Rcpp::traits::input_parameter< const double >::type a0(a0SEXP);
+    Rcpp::traits::input_parameter< const double >::type u0(u0SEXP);
+    Rcpp::traits::input_parameter< const NumericVector >::type r1(r1SEXP);
+    Rcpp::traits::input_parameter< const NumericVector >::type r2(r2SEXP);
+    Rcpp::traits::input_parameter< const NumericVector >::type r3(r3SEXP);
+    Rcpp::traits::input_parameter< const NumericVector >::type f4(f4SEXP);
+    rcpp_result_gen = Rcpp::wrap(rcpp_calcEmis(ct, a0, u0, r1, r2, r3, f4));
+    return rcpp_result_gen;
+END_RCPP
+}
+
+static const R_CallMethodDef CallEntries[] = {
+    {"_ALFAM2_rcpp_calcEmis", (DL_FUNC) &_ALFAM2_rcpp_calcEmis, 7},
+    {NULL, NULL, 0}
+};
+
+RcppExport void R_init_ALFAM2(DllInfo *dll) {
+    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
+    R_useDynamicSymbols(dll, FALSE);
+}

src/rcpp_calcEmis.cpp

@@ -0,0 +1,54 @@
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// [[Rcpp::export]]
+List rcpp_calcEmis(const NumericVector ct, const double a0, 
+                   const double u0, const NumericVector r1, 
+                   const NumericVector r2, const NumericVector r3, 
+                   const NumericVector f4) {
+    //Number of intervals
+    R_xlen_t l = ct.length();
+    NumericVector a(l), u(l), e(l), ddt(l);
+
+    //time steps
+    ddt[0] = ct[0] - 0;
+    for (R_xlen_t i = 1; i < l; ++i) {
+      ddt[i] = ct[i] - ct[i-1];
+    }
+
+    //initial values
+    double ati0 = a0;
+    double uti0 = u0;
+    double eti = 0;
+    double ati, uti;
+
+    for (R_xlen_t i = 0; i < l; ++i) {
+      // Make incorporation transfer (at *start* of interval) (if none then f4 = 1 and ati = a[i])
+      ati = f4[i] * ati0;
+      uti = (1 - f4[i]) * ati0 + uti0;
+
+      //Calculate pools at *end* of ct[i]
+      a[i] = ati * exp(-(r1[i] + r2[i]) * ddt[i]);
+      u[i] = exp(-r3[i] * ddt[i]) * (r2[i] * ati * (exp((-r1[i] - r2[i] + r3[i]) * ddt[i]) - 1.0)/(-r1[i] - r2[i] + r3[i]) + uti);
+      e[i] = eti + (uti - u[i]) +  (ati - a[i]);
+
+      //save pools for next step
+        ati0 = a[i];
+        uti0 = u[i];
+        eti = e[i];
+    }
+
+    return List::create(_["ct"] = ct, 
+                        _["dt"] = ddt,
+                        _["f0"] = a0/(u0 + a0),
+                        _["r1"] = r1, 
+                        _["r2"] = r2,
+                        _["r3"] = r3, 
+                        _["f4"] = f4, 
+                        _["f"] = a, 
+                        _["s"] = u, 
+                        _["j"] = NA_REAL, 
+                        _["e"] = e, 
+                        _["e.int"] = NA_REAL, 
+                        _["er"] = e/(a0 + u0));
+}