partial_sum(部分和)とadjacent_difference(隣接差)
partial_sumはaccumulateのバリエーション、accumulateがX[N]の総和を求めるのに対し、partian_sumはY[i]=X[0]からX[i]までの和となるよう、Y[i]を埋めてくれます。
template<class InputIterator, OutputIterator>
OutputIterator partial_sum(InputIterator first, InputIterator last,
OutputIterator result);
template<class InputIterator, OutputIterator,
class BinaryOperation>
OutputIterator partial_sum(InputIterator first, InputIterator last,
OutputIterator result
BinaryOperation binary_op);
adjacent_differenceはpartial_sumの逆、Y[i] = X[i] - X[i-1]、つまり直前の値との差を求めます(ただしY[0] = X[0])。
ではサンプル、有理数列の和1/2 + 1/4 + 1/8 + ……が次第に1に近づく様子を観察します。
#include <numeric>
#include <iostream>
#include <vector>
#include <iterator>
#include "sequence.h"
#include "rational.h"
template<typename T>
std::ostream& operator<<(std::ostream& stream, const rational<T>& r)
{ return stream << right << r.first << '/' << left << r.second; }
using namespace std;
// 1/2 + 1/4 + 1/8 + 1/16 …… は 次第に 1 に近づく
int main() {
const size_t N = 16;
vector<double> c(N);
iota(begin(c), end(c), sequence<double>(0.5, [](double x) { return x*0.5;}));
// 第i項までの和を dout に求める
vector<double> dout;
partial_sum(begin(c), end(c), back_inserter(dout));
for ( auto item : dout ) {
cout << item << endl;
}
// 同じことを有理数について行う
vector<rational<>> r(N, make_rational(0L));
// r[] = 1/2, 1/4, 1/8 ……
iota(begin(r), end(r),
sequence<rational<>>(make_rational(1L,2L),
[](const rational<>& x)
{ return make_rational(numerator(x),denominator(x)*2L); }));
// 第i項までの和を rout に求める
vector<rational<>> rout;
partial_sum(begin(r), end(r), back_inserter(rout),
[](const rational<>& a, const rational<>& b) // 有理数の和
{ return make_rational(numerator(a)*denominator(b) + numerator(b)*denominator(a),
denominator(a)*denominator(b)); });
for ( auto item : rout ) {
cout << left << setw(10) << dout[i] << setw(10) << rout[i] << endl;
}
}
以前紹介した並列STLにも、accumulate,partial_sumに相当する並列アルゴリズムが導入されています。それが reduceとinclusive_scan、さらにinclusive_scanの姉妹品(?)exclusive_scan。inclusive_scanとexclusive_scanとはちょっとだけ違っていて、入力をX[],得られた部分和(の列)をY[]とすると、
inclusive_scanだと、
Y[0] = X[0] Y[1] = X[0] + X[1] ... Y[i] = X[0] + X[1] + .... + X[i]
対してexclusive_scanは、
Y[0] = 0 Y[1] = 0 + X[0] Y[2] = 0 + X[0] + X[1] ... Y[i] = 0 + X[0] + X[1] + .... + X[i]
つまり、exclusive_scanのi番目の結果Y[i]はX[i]を含まず、X[i-1]までの和になってます。
reduce, inclusive_scan(exclusive_scan)は複数のスレッドで並列に処理を行います。複数の演算が同時に行われるため、演算の順序が一意に定まりません。従ってreduce, inclusive_scan(exclusive_scan)内で呼び出されるoperator+, operator*および二項関数オブジェクト:binary_opが、
交換則:
x + y = y + x x * y = y * x binary_op(x,y) = binary_op(y,x)
結合則:
(x + y) + z = x + (y + z) (x * y) * z = x * (y * z) binary_op(binary_op(x,y),z) = binary_op(x,binary_op(y,z))
を満たさないと結果の一意性が保障されません。
exclusive_scanを使ったちょっと面白いサンプル:移動平均(simple moving average)を紹介します。
exclusive_scanの結果Y[i]はX[0]からX[i-1]までの和、Y[i]からwだけ離れたY[i+w]にはX[0]からX[i+w-1]までの和ですから、その差:Y[i+w]-Y[i]はX[i]からX[i+w-1]までの和となります。これをwで割るとX[i]からw個分の平均が求まるってスンポーです。
#include <experimental/algorithm>
#include <experimental/numeric>
#include <iostream>
#include <algorithm>
#include <random>
#include <vector>
#include <iterator>
using namespace std;
using namespace std::experimental;
int main() {
mt19937 gen; // メルセンヌ・ツイスタ
normal_distribution<float> dist; // (平均0,標準偏差1の)気温の揺れ
float t = 25.0f;
const size_t N = 31;
vector<float> temp; // 日々の気温データ
generate_n(back_inserter(temp), N, [&]() { return t += dist(gen);});
vector<float> sum; // 部分和
parallel::exclusive_scan(parallel::par, begin(temp), end(temp), back_inserter(sum), 0.0f);
vector<float> sma; // simple-moving-average
size_t w = 3; // window size
parallel::transform(parallel::par,
begin(sum)+w, end(sum), begin(sum), back_inserter(sma),
[=](float x, float y) { return (x-y)/w; });
for ( unsigned int i = 0U; i < sma.size(); ++i ) {
cout << temp[i] << '\t' << sma[i] << endl;
}
}
……<numeric>が提供してくれている関数はたったの5つながら、従来for-loopでくるくる回していた計算が関数呼び出しイッパツで済むのは快感です。サンプルで示したとおり、統計がらみの計算にはうってつけですね。

