Here is something faster. It avoids resizing/copying x
in the loop.
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericVector combine_points_Roland3(NumericVector x, double th) {
int i=0, j;
int n(x.size());
while(i < n-1)
{
if ((x[i+1] - x[i]) >= th)
{
i = i + 1;
}
else{
x[i] = (x[i+1] + x[i])/2;
n = n-1;
for (j=i+1; j<n; j++)
{
x[j]=x[j+1];
}
}
}
NumericVector y(n);
for (i = 0; i < n; i++) {
y[i] = x[i];
}
return y;
}
An R implementation of the same algorithm:
combine_points_Roland3R <- function(x, th) {
i <- 1
n <- length(x)
while(i < n) {
if ((x[i+1] - x[i]) >= th) {
i <- i + 1;
} else {
x[i] <- (x[i+1] + x[i])/2
n <- n-1
x[(i+1):n] <- x[(i+2):(n+1)]
}
}
x[1:n]
}
set.seed(41); d <- seq(1:4000) + runif(4000, 0, 1)
x2 <- combine_points_Roland2(d, 0.5)
x3 <- combine_points_Roland3(d, 0.5)
all.equal(x2, x3)
#TRUE
x4 <- combine_points_Roland3R(d, 0.5)
all.equal(x2, x4)
#TRUE
Benchmarks:
library(microbenchmark)
microbenchmark(combine_points_Roland2(d, 0.5),
combine_points_Roland3(d, 0.5),
combine_points_Roland3R(d, 0.5))
# Unit: microseconds
# expr min lq median uq max neval
# combine_points_Roland2(d, 0.5) 126458.64 131414.592 132355.4285 133422.2235 147306.728 100
# combine_points_Roland3(d, 0.5) 121.34 128.269 140.8955 143.3595 393.582 100
# combine_points_Roland3R(d, 0.5) 17564.24 18626.878 19155.6565 20910.2935 68707.888 100