1 #include <thrust/device_vector.h>
2 #include <thrust/scan.h>
3 #include <thrust/transform.h>
4 #include <thrust/functional.h>
5 #include <thrust/sequence.h>
6 #include <thrust/random.h>
7 #include <iostream>
8 #include <iomanip>
9
10 // Efficiently computes the simple moving average (SMA) [1] of a data series
11 // using a parallel prefix-sum or "scan" operation.
12 //
13 // Note: additional numerical precision should be used in the cumulative summing
14 // stage when computing the SMA of large data series. The most straightforward
15 // remedy is to replace 'float' with 'double'. Alternatively a Kahan or
16 // "compensated" summation algorithm could be applied [2].
17 //
18 // [1] http://en.wikipedia.org/wiki/Moving_average#Simple_moving_average
19 // [2] http://en.wikipedia.org/wiki/Kahan_summation_algorithm
20
21
22 // compute the difference of two positions in the cumumulative sum and
23 // divide by the SMA window size w.
24 template <typename T>
25 struct minus_and_divide : public thrust::binary_function<T,T,T>
26 {
27 T w;
28
minus_and_divideminus_and_divide29 minus_and_divide(T w) : w(w) {}
30
31 __host__ __device__
operator ()minus_and_divide32 T operator()(const T& a, const T& b) const
33 {
34 return (a - b) / w;
35 }
36 };
37
38 template <typename InputVector, typename OutputVector>
simple_moving_average(const InputVector & data,size_t w,OutputVector & output)39 void simple_moving_average(const InputVector& data, size_t w, OutputVector& output)
40 {
41 typedef typename InputVector::value_type T;
42
43 if (data.size() < w)
44 return;
45
46 // allocate storage for cumulative sum
47 thrust::device_vector<T> temp(data.size() + 1);
48
49 // compute cumulative sum
50 thrust::exclusive_scan(data.begin(), data.end(), temp.begin());
51 temp[data.size()] = data.back() + temp[data.size() - 1];
52
53 // compute moving averages from cumulative sum
54 thrust::transform(temp.begin() + w, temp.end(), temp.begin(), output.begin(), minus_and_divide<T>(T(w)));
55 }
56
main(void)57 int main(void)
58 {
59 // length of data series
60 size_t n = 30;
61
62 // window size of the moving average
63 size_t w = 4;
64
65 // generate random data series
66 thrust::device_vector<float> data(n);
67 thrust::default_random_engine rng;
68 thrust::uniform_int_distribution<int> dist(0, 10);
69 for (size_t i = 0; i < n; i++)
70 data[i] = static_cast<float>(dist(rng));
71
72 // allocate storage for averages
73 thrust::device_vector<float> averages(data.size() - (w - 1));
74
75 // compute SMA using standard summation
76 simple_moving_average(data, w, averages);
77
78 // print data series
79 std::cout << "data series: [ ";
80 for (size_t i = 0; i < data.size(); i++)
81 std::cout << data[i] << " ";
82 std::cout << "]" << std::endl;
83
84 // print moving averages
85 std::cout << "simple moving averages (window = " << w << ")" << std::endl;
86 for (size_t i = 0; i < averages.size(); i++)
87 std::cout << " [" << std::setw(2) << i << "," << std::setw(2) << (i + w) << ") = " << averages[i] << std::endl;
88
89 return 0;
90 }
91
92