Самый быстрый способ получить неотрицательный компонент

Question

Что такое более быстрый способ получить неотрицательную составляющую двойного вектора? То есть,

pmax(x, 0) 

Моя попытка использования Rcpp:

//' @title Parallel maximum 
//' @description A faster \code{pmax()}. 
//' 
//' @name pmaxC 
//' @param x A numeric vector. 
//' @param a A single numeric value. 
//' @return The parallel maximum of the input values. 
//' @note This function will always be faster than \code{pmax(x, a)} when \code{a} is a single value, but can be slower than \code{pmax.int(x, a)} when \code{x} is short. Use this function when comparing a numeric vector with a single value. 
//' @export pmaxC 

#include <Rcpp.h> 
using namespace Rcpp; 

// [[Rcpp::export]] 
NumericVector pmaxC(NumericVector x, double a) { 
    int n = x.length(); 
    NumericVector out(n); 

    for (int i = 0; i < n; ++i) { 
    double xi = x[i]; 
    if (xi < a) { 
     out[i] = a; 
    } else { 
     out[i] = xi; 
    } 
    } 

    return out; 
} 

Это скромное улучшение:

set.seed(5) 
x <- rnorm(1e6) 

microbenchmark(pmax(x, 0), pmaxC(x, 0)) 
Unit: milliseconds 
     expr  min  lq  mean median  uq  max neval cld 
    pmax(x, 0) 8.500419 8.621341 11.09672 10.132045 10.791020 58.44972 100 a 
pmaxC(x, 0) 5.624480 5.709262 8.83968 7.598093 7.907853 53.91339 100 a 

Ни неприемлемо медленно, но, учитывая, что общий сценарий, я было интересно, развил ли пакет более быстрый подход.

Answer 1

Операция, которую вы выполняете, довольно проста, поэтому я не уверен, что есть много возможностей для улучшения в отношении вашего алгоритма выше. Однако, если вы действительно нужно выжать дополнительную производительность, это кажется хорошим кандидатом для распараллеливания. Вот возможный вариант осуществления с использованием RcppParallel:

// [[Rcpp::depends(RcppParallel)]] 
#include <RcppParallel.h> 
#include <Rcpp.h> 

struct Pmax : public RcppParallel::Worker { 
    struct Apply { 
     double mx; 
     Apply(double mx_) 
      : mx(mx_) 
     {} 

     double operator()(const double x) const 
     { 
      return x > mx ? x : mx; 
     } 
    }; 

    const RcppParallel::RVector<double> input; 
    RcppParallel::RVector<double> output; 

    Apply f; 

    Pmax(const Rcpp::NumericVector input_, 
     Rcpp::NumericVector output_, 
     double mx_) 
    : input(input_), output(output_), f(mx_) 
    {} 

    void operator()(std::size_t begin, std::size_t end) 
    { 
     std::transform(
      input.begin() + begin, 
      input.begin() + end, 
      output.begin() + begin, 
      f 
     ); 
    } 
}; 

// [[Rcpp::export]] 
Rcpp::NumericVector par_pmax(Rcpp::NumericVector x, double y) 
{ 
    Rcpp::NumericVector res = Rcpp::no_init_vector(x.size()); 
    Pmax p(x, res, y); 
    RcppParallel::parallelFor(0, x.size(), p); 

    return res; 
} 

---

Тестирование это с примера данных, я получаю разумное улучшение:

set.seed(5) 
x <- rnorm(1e6) 

all.equal(pmax(x, 0), par_pmax(x, 0)) 
#[1] TRUE 

microbenchmark::microbenchmark(
    pmax(x, 0), 
    pmaxC(x, 0), 
    par_pmax(x, 0), 
    times = 500L 
) 
# Unit: milliseconds 
#   expr  min  lq  mean median  uq  max neval 
#  pmax(x, 0) 11.843528 12.193126 14.972588 13.030448 16.799250 102.09895 500 
#  pmaxC(x, 0) 7.804883 8.036879 10.462070 8.772635 12.407587 69.08290 500 
# par_pmax(x, 0) 2.244691 2.443971 4.552169 2.624008 6.359027 65.99233 500