Скользящее среднее функции в C++ для использования С Р


Я пытаюсь улучшить код и ускорить в C++ (Rcpp) а (по центру) взвешенное скользящее среднее функции я закодированный.

Пример того, что функция roll_mean делает. Обратите внимание, что функция работает независимо от того, какого размера х и адаптируется к обоим хвосты мои данные

w=c(1/2,1,1/2)
x=c(4,2,6,12)
res=c(2,5,7,3) 
res=c((x[1:2]*w[2:3])/sum(w[2:3]),x[1:3]*w[1:3]/sum(w[1:3]),x[2:4]*w[1:3]/sum(w[1:3]),x[3:4]*w[1:2]/sum(w[1:2]))

Файл PartialMA.cpp

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector roll_mean(const NumericVector& x,
                        const NumericVector& w) {

  int n = x.size();
  int w_size = w.size();
  int size = (w_size - 1) / 2;

  NumericVector res(n);
  int i, ind_x, ind_w;

  double tmp_wsum, tmp_xwsum, tmp_w;

  for (i = 0; i < n; i++) {
    tmp_xwsum = 0;
    tmp_wsum = 0;
    for (ind_x = i - size, ind_w = 0; ind_x < i + size; ind_x++, ind_w++) {
      if((ind_x >= 0) & (ind_x < n)){
      tmp_w = w(ind_w);
      tmp_xwsum += x(ind_x) *  tmp_w;
      tmp_wsum += tmp_w;
      }
    }
    res[i] = tmp_xwsum / tmp_wsum;
  }

  return res;
}

Я пытался заменить петли +, если заявление с этой целью минимизации количества итераций:

for (ind_x = std::max(0, i - size), ind_w = std::max(0, size-1); ind_x < std::min(n, i + size); ind_x++, ind_w++) {

Я чувствую, что я не достаточно жестким, и я буду очень признателен, если кто-то может помочь мне улучшить мой код и в конечном итоге ускорить, насколько возможно, функции.



Комментарии
2 ответа

В целом ваш код-это не плохо, но мы можем сделать лучше. Начнем со стилем кодирования.

Ваши переменные будут четко названы и ваш код прост и понятен, однако использование & Не надо.

const NumericVector &x

R объекты всегда передаются по ссылке Rcpp даже без &. Дирк имеет некоторые слайды (см. 29,30) с более подробной информацией об этой теме здесь. Также в Rcpp часто задаваемые вопросы. Дело в том, что мы можем смело отбросить эти амперсанды.

Теперь, давайте поговорим о том, как мы можем улучшить производительность вашего кода. Сначала заметим, что есть много ненужных проверок во внутренних for loop.

if((ind_x >= 0) & (ind_x < n))

Этого можно избежать за пределами петли вот так :

// [[Rcpp::export]]
NumericVector roll_meanReduceChecks(const NumericVector x,
const NumericVector w) {

int n = x.size();
int w_size = w.size();
int size = (w_size - 1) / 2;

NumericVector res(n);
int i, ind_x, ind_w, strt, endx;

double tmp_wsum, tmp_xwsum, tmp_w;

for (i = 0; i < n; i++) {
tmp_xwsum = 0;
tmp_wsum = 0;

if ((i - size) <= 0) {
strt = 0;
ind_w = size - i;
} else {
strt = i - size;
ind_w = 0;
}

endx = ((i + size) >= n) ? n : (i + size);

for (ind_x = strt; ind_x < endx; ind_x++, ind_w++) {
// This check is no longer necessary
// if((ind_x >= 0) & (ind_x < n)){
tmp_w = w(ind_w);
tmp_xwsum += x(ind_x) * tmp_w;
tmp_wsum += tmp_w;
}

res[i] = tmp_xwsum / tmp_wsum;
}

return res;
}

С этой модификацией, мы получаем около 10% быстрее тайминги:

set.seed(42)
x <- sample(10^3, 10^4, TRUE)
w <- sample(100, 10^3 , TRUE) / 100

## Gives the same result
all.equal(roll_mean(x, w), roll_meanReduceChecks(x, w))
[1] TRUE

library(microbenchmark)
microbenchmark(roll_mean(x, w), roll_meanReduceChecks(x, w),
times = 50, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
roll_mean(x, w) 1.138508 1.135057 1.119695 1.132882 1.123308 1.149893 50
roll_meanReduceChecks(x, w) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 50

Мы все еще можем сделать лучше. Реальная экономия вступает в игру с использованием iterators. Есть много хорошей информации там по этой теме, и вот только один ресурс написал Хэдли Уикхэм, который является хорошей отправной точкой.

Собираем все вместе мы:

// [[Rcpp::export]]
NumericVector roll_meanIterator(const NumericVector x,
const NumericVector w) {

int n = x.size();
int w_size = w.size();
int size = (w_size - 1) / 2;

NumericVector res(n);
NumericVector::const_iterator itw, itx, itEnd;
double tmp_wsum, tmp_xwsum;

for (int i = 0; i < n; i++) {
tmp_xwsum = 0;
tmp_wsum = 0;

if ((i - size) <= 0) {
itx = x.begin();
itw = w.begin() + size - i;
} else {
itx = x.begin() + i - size;
itw = w.begin();
}

itEnd = ((i + size) >= n) ? x.end() : x.begin() + i + size;

for (; itx < itEnd; itx++, itw++) {
tmp_xwsum += (*itx) * (*itw);
tmp_wsum += (*itw);
}

res[i] = tmp_xwsum / tmp_wsum;
}

return res;
}

Вот это проверка на вменяемость:

all.equal(roll_mean(x, w), roll_meanIterator(x, w))
[1] TRUE

И вот несколько критериев:

microbenchmark(roll_mean(x, w), roll_meanIterator(x, w), 
times = 50, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
roll_mean(x, w) 8.480686 8.458674 8.493975 8.894612 8.548443 8.453865 50
roll_meanIterator(x, w) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 50

Это почти 9x быстрее. Не плохо, учитывая, что мы не должны изменить код, что очень.

Обновление

Ниже приведен пример по оптимизации и урок по теме "когда нужно остановиться" , так как код является более эффективным, но это гораздо сложнее, чем код выше с практически нет компромисса в эффективности.

Если мы рассмотрим tmp_wsum более пристально, мы обнаружим, что на каждой итерации, он только меняется на концах. В roll_meanReduceChecks выше, если мы вставим std::cout << ind_w << ' '; просто перед внутренней for loop и std::cout << ind_w << std::endl; только после этого, вот результат за первые 5 итераций:

499 998
498 998
497 998
496 998
495 998

И последние 5 итераций:

0 504
0 503
0 502
0 501
0 500

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

NumericVector roll_meanOverOptimized(const NumericVector x,
const NumericVector w) {

int n = x.size();
int w_size = w.size();
int size = (w_size - 1) / 2;

NumericVector res(n); // must add variables for checks below
NumericVector::const_iterator itw, itwBeg, itwEnd, itx, itxEnd;
double tmp_wsum = 0, tmp_xwsum = 0;

// We must first initialize tmp_wsum, itwBeg,
// itwEnd, as well as populate res[0]
itx = x.begin();
itwBeg = itw = w.begin() + size;
itxEnd = (size >= n) ? x.end() : x.begin() + size;

for (; itx < itxEnd; itx++, itw++) {
tmp_xwsum += (*itx) * (*itw);
tmp_wsum += (*itw);
}

res[0] = tmp_xwsum / tmp_wsum;
itwEnd = itw;

// Start i @ 1 instead of 0 as the first
// iteration was taken care of above
for (int i = 1; i < n; i++) {
tmp_xwsum = 0;

if ((i - size) < 0) {
itx = x.begin();
itw = w.begin() + size - i;
} else {
itx = x.begin() + i - size;
itw = w.begin();
}

// first check
if (itw != itwBeg)
tmp_wsum += (*itw);
itwBeg = itw;

itxEnd = ((i + size) > n) ? x.end() : x.begin() + i + size;

// N.B. only one variable is being updated now
for (; itx < itxEnd; itx++, itw++)
tmp_xwsum += (*itx) * (*itw);

// second check
if (itw != itwEnd)
tmp_wsum -= (*itw);
itwEnd = itw;

res[i] = tmp_xwsum / tmp_wsum;
}

return res;
}

Давайте посмотрим, если эта дополнительная оптимизация окупились:

## Gives the same results
all.equal(roll_meanIterator(x, w), roll_meanOverOptimized(x, w))
[1] TRUE

microbenchmark(roll_mean(x, w),
roll_meanIterator(x, w),
roll_meanOverOptimized(x, w),
times = 50, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
roll_mean(x, w) 8.679318 8.751559 8.605352 8.796965 8.509649 7.683454 50
roll_meanIterator(x, w) 1.015420 1.028507 1.041876 1.019177 1.017208 1.133470 50
roll_meanOverOptimized(x, w) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000 50

Практически идентичные результаты, как и roll_meanIterator. ИМО, ~2% прирост эффективности не стоит, так как мы сделали наш код более трудным для чтения.

Обновление 2 (Большое спасибо @HongOoi)

В комментариях @HongOoi государства, "ускорение не потому, что итераторы как таковой, а потому, что ОП использует с контролем границ массива доступ через оператора (). Переключение в обычный [] должно дать вам тот же результат, как итераторы. Давайте тест!!!

// [[Rcpp::export]]
NumericVector roll_meanBrackets(const NumericVector x,
const NumericVector w) {

int n = x.size();
int w_size = w.size();
int size = (w_size - 1) / 2;

NumericVector res(n);
double tmp_wsum, tmp_xwsum;
unsigned long int ind_x, ind_w, xEnd;

for (int i = 0; i < n; i++) {
tmp_xwsum = 0;
tmp_wsum = 0;

if ((i - size) < 0) {
ind_x = 0;
ind_w = size - i;
} else {
ind_x = i - size;
ind_w = 0;
}

xEnd = ((i + size) > n) ? x.size() : i + size;

for (; ind_x < xEnd; ind_x++, ind_w++) {
tmp_xwsum += x[ind_x] * w[ind_w];
tmp_wsum += w[ind_w];
}

res[i] = tmp_xwsum / tmp_wsum;
}

return res;
}

И вот несколько критериев:

microbenchmark(roll_mean(x, w),
roll_meanIterator(x, w),
roll_meanBrackets(x, w),
times = 50, unit = "relative")
Unit: relative
expr min lq mean median uq max neval
roll_mean(x, w) 8.596342 8.425527 8.374452 8.414136 8.3626927 7.479161 50
roll_meanIterator(x, w) 1.012213 1.007219 1.007606 1.000564 0.9973372 1.177148 50
roll_meanBrackets(x, w) 1.000000 1.000000 1.000000 1.000000 1.0000000 1.000000 50

Как и предсказывал @HongOoi, эффективность практически идентична. Напрашивается вопрос


Какой метод следует использовать? Итераторы или C-стиль индексации?

К счастью для нас, есть некоторые большие информацию там, обращаясь к этой самой теме. На самом деле, я нашел ответ на сайте StackOverflow вопрос итератора цикла против индекса цикла очень полезно. Вот резюме по @TemplateRex:


"если вам действительно нужен индекс (например, доступ к предыдущему или следующему элементу, печать/лесозаготовки индекс внутри цикла и т. д.) или вам нужен Страйд разных, чем 1, то я хотел бы пойти на явно проиндексированных контур, иначе я пошел бы по серии-для петли.

Для обобщенных алгоритмов на универсальных контейнеров, я бы за прямого итератора цикла, если код не содержит никаких регулирование потока внутри цикла и необходимый шаг 1, в этом случае я бы пошел на стл for_each + лямбда."

6
ответ дан 26 марта 2018 в 02:03 Источник Поделиться

Как правило, средневзвешенная веса должен добавить к 1, а в образце он добавляет к 2

w=c(1/2,1,1/2)

sum(w)=2

может, так и должно быть?

w=c(1/4,1/2,1/4) 

и вы можете получите скользящую среднюю с функцией фильтра

F <- filter(x, filter = w, method = c("convolution"), sides = 2)

Если ваши ядра свертки w слишком большой, и вы хотите проверить его на скорость, я бы попробовал быстрого преобразования Фурье свертки. БПФ свертки должны быть уже где-то на всех языках.

convolve(x, w, conj = TRUE, type = c( "open"))

БПФ обладает свойством

#pseudocode
FFT(F) = FFT(x) * FFT(w)

так, чтобы получить F вы

#pseudocode
F <- inverseFFT( FFT(x) * FFT(w) )

Некоторые недостатки с БПФ, что


  • обычно он должен длина(х) == 2^n (где N целое)

  • оно имеет периодический характер.

Эти проблемы, как правило, должны быть решены путем заполнения данных X с 0

Кроме того, в коде Rcpp, часто средневзвешенная ядра (ж) является симметричной и имеет неоднократные Весов, так что вы можете воспользоваться его сохранения умножения, чтобы не повторять их. Я не делаю из кода C++, но кто-то другой может использовать его, чтобы улучшить ваш код.

2
ответ дан 10 апреля 2018 в 03:04 Источник Поделиться