C++ from http://adv-r.had.co.nz/Rcpp.html

agentBasedModel.cpp

// http://adv-r.had.co.nz/Rcpp.html#rcpp-classes
 
#include <Rcpp.h>
using namespace Rcpp;
 
double vacc3a(double age, bool female, bool ily){
  double p = 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily;
  p = p * (female ? 1.25 : 0.75);
  p = std::max(p, 0.0);
  p = std::min(p, 1.0);
  return p;
}
 
// [[Rcpp::export]]
NumericVector vacc3(NumericVector age, LogicalVector female, 
                    LogicalVector ily) {
  int n = age.size();
  NumericVector out(n);
 
  for(int i = 0; i < n; ++i) {
    out[i] = vacc3a(age[i], female[i], ily[i]);
  }
 
  return out;
}
 
/*** R
vacc1a <- function(age, female, ily) {
  p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
  p <- p * if (female) 1.25 else 0.75
  p <- max(0, p)
  p <- min(1, p)
  p
}
vacc1 <- function(age, female, ily) {
  n <- length(age)
  out <- numeric(n)
  for (i in seq_len(n)) {
    out[i] <- vacc1a(age[i], female[i], ily[i])
  }
  out
}
vacc2 <- function(age, female, ily) {
  p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
  p <- p * ifelse(female, 1.25, 0.75)
  p <- pmax(0, p)
  p <- pmin(1, p)
  p
}
vaccMe <- function(age, female, ily) {
  p <- 0.25 + 0.3 * 1 / (1 - exp(0.04 * age)) + 0.1 * ily
  p <- p + 0.75 + 0.5 * female // ifelse(female, 1.25, 0.75)
  p <- p[p < 0] <- 0 // pmax(0, p)
  p <- p[p > 1] <- 1 // pmin(1, p)
  p
}

n <- 1000
age <- rnorm(n, mean = 50, sd = 10)
female <- sample(c(T, F), n, rep = TRUE)
ily <- sample(c(T, F), n, prob = c(0.8, 0.2), rep = TRUE)

stopifnot(
#  all.equal(vacc1(age, female, ily), vacc1a(age, female, ily)),
  all.equal(vacc1(age, female, ily), vacc2(age, female, ily))
# ,all.equal(vacc1(age, female, ily), vaccMe(age, female, ily))
  ,all.equal(vacc1(age, female, ily), vacc3(age, female, ily))
)

microbenchmark(
  vacc1a = vacc1a(age, female, ily),
  vacc1 = vacc1(age, female, ily),
  vacc2 = vacc2(age, female, ily),
  vacc3 = vacc3(age, female, ily),
  vaccMe = vaccMe(age, female, ily)
)
*/

gibbs.cpp

// http://adv-r.had.co.nz/Rcpp.html#rcpp-classes
 
#include <Rcpp.h>
using namespace Rcpp;
 
// [[Rcpp::export]]
NumericMatrix gibbs_cpp(int N, int thin) {
  NumericMatrix mat(N, 2);
  double x = 0, y = 0;
 
  for(int i = 0; i < N; i++) {
    for(int j = 0; j < thin; j++) {
      x = rgamma(1, 3, 1 / (y * y + 4))[0];
      y = rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))[0];
    }
    mat(i, 0) = x;
    mat(i, 1) = y;
  }
 
  return(mat);
}
 
/*** R
gibbs_r <- function(N, thin) {
  mat <- matrix(nrow = N, ncol = 2)
  x <- y <- 0

  for (i in 1:N) {
    for (j in 1:thin) {
      x <- rgamma(1, 3, y * y + 4)
      y <- rnorm(1, 1 / (x + 1), 1 / sqrt(2 * (x + 1)))
    }
    mat[i, ] <- c(x, y)
  }
  mat
}

microbenchmark(
  gibbs_r(100, 10),
  gibbs_cpp(100, 10)
)
*/

index.html

...

stats.cpp

#include <Rcpp.h>
using namespace Rcpp;
 
// [[Rcpp::export]]
double f1(NumericVector x) {
  int n = x.size();
  double y = 0;
 
  for(int i = 0; i < n; ++i) {
    y += x[i] / n;
  }
  return y;
}
 
// [[Rcpp::export]]
NumericVector f2(NumericVector x) {
  int n = x.size();
  NumericVector out(n);
 
  out[0] = x[0];
  for(int i = 1; i < n; ++i) {
    out[i] = out[i - 1] + x[i];
  }
  return out;
}
 
// [[Rcpp::export]]
bool f3(LogicalVector x) {
  int n = x.size();
 
  for(int i = 0; i < n; ++i) {
    if (x[i]) return true;
  }
  return false;
}
 
// [[Rcpp::export]]
int f4(Function pred, List x) {
  int n = x.size();
 
  for(int i = 0; i < n; ++i) {
    LogicalVector res = pred(x[i]);
    if (res[0]) return i + 1;
  }
  return 0;
}
 
// [[Rcpp::export]]
NumericVector f5(NumericVector x, NumericVector y) {
  int n = std::max(x.size(), y.size());
  NumericVector x1 = rep_len(x, n);
  NumericVector y1 = rep_len(y, n);
 
  NumericVector out(n);
 
  for (int i = 0; i < n; ++i) {
    out[i] = std::min(x1[i], y1[i]);
  }
 
  return out;
}

sugar.cpp

// http://adv-r.had.co.nz/Rcpp.html#rcpp-classes
 
#include <Rcpp.h>
using namespace Rcpp;
 
// [[Rcpp::plugins(cpp11)]]
#include <Rcpp.h>
#include <unordered_set>
using namespace Rcpp;
 
#include <Rcpp.h>
using namespace Rcpp;
 
// [[Rcpp::export]]
std::map<double, int> tableC(NumericVector x) {
  std::map<double, int> counts;
 
  int n = x.size();
  for (int i = 0; i < n; i++) {
    counts[x[i]]++;
  }
 
  return counts;
}
/*** R
x0 <- runif(1e5)
microbenchmark(tableC(x0))
str(tableC(x0))
*/
 
// [[Rcpp::export]]
LogicalVector duplicatedC(IntegerVector x) {
  std::unordered_set<int> seen;
  int n = x.size();
  LogicalVector out(n);
 
  for (int i = 0; i < n; ++i) {
    out[i] = !seen.insert(x[i]).second;
  }
 
  return out;
}
/*** R#
x0 <- runif(1e5)
microbenchmark(duplicatedC(x0))
str(duplicatedC(x0))
*/
 
// [[Rcpp::export]]
List rleC(NumericVector x) {
  std::vector<int> lengths;
  std::vector<double> values;
 
  // Initialise first value
  int i = 0;
  double prev = x[0];
  values.push_back(prev);
  lengths.push_back(1);
 
  NumericVector::iterator it;
  for(it = x.begin() + 1; it != x.end(); ++it) {
    if (prev == *it) {
      lengths[i]++;
    } else {
      values.push_back(*it);
      lengths.push_back(1);
 
      i++;
      prev = *it;
    }
  }
 
  return List::create(
    _["lengths"] = lengths, 
    _["values"] = values
  );
}
 
 
// [[Rcpp::export]]
bool any_naC(NumericVector x) {
  return is_true(any(is_na(x)));
}
 
/*** 
x0 <- runif(1e5)

microbenchmark(rleC(x0))
str(rleC(x0))
*/
/*** 
x1 <- c(x0, NA)
x2 <- c(NA, x0)
any_naR <- function(x) any(is.na(x))

microbenchmark(
  any_naR(x0), any_naC(x0),
  any_naR(x1), any_naC(x1),
  any_naR(x2), any_naC(x2)
)
*/
 
// [[Rcpp::export]]
double sum3(NumericVector x) {
  double total = 0;
  
  NumericVector::iterator it;
  for(it = x.begin(); it != x.end(); ++it) {
    total += *it;
  }
  return total;
}

z C++ and R examples.md

Summary goes here.