1 #include <thrust/device_vector.h>
2 #include <thrust/host_vector.h>
3 #include <thrust/transform_reduce.h>
4 #include <thrust/functional.h>
5 #include <thrust/extrema.h>
6 #include <cmath>
7 #include <limits>
8 #include <iostream>
9
10 // This example computes several statistical properties of a data
11 // series in a single reduction. The algorithm is described in detail here:
12 // http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
13 //
14 // Thanks to Joseph Rhoads for contributing this example
15
16
17 // structure used to accumulate the moments and other
18 // statistical properties encountered so far.
19 template <typename T>
20 struct summary_stats_data
21 {
22 T n;
23 T min;
24 T max;
25 T mean;
26 T M2;
27 T M3;
28 T M4;
29
30 // initialize to the identity element
initializesummary_stats_data31 void initialize()
32 {
33 n = mean = M2 = M3 = M4 = 0;
34 min = std::numeric_limits<T>::max();
35 max = std::numeric_limits<T>::min();
36 }
37
variancesummary_stats_data38 T variance() { return M2 / (n - 1); }
variance_nsummary_stats_data39 T variance_n() { return M2 / n; }
skewnesssummary_stats_data40 T skewness() { return std::sqrt(n) * M3 / std::pow(M2, (T) 1.5); }
kurtosissummary_stats_data41 T kurtosis() { return n * M4 / (M2 * M2); }
42 };
43
44 // stats_unary_op is a functor that takes in a value x and
45 // returns a variace_data whose mean value is initialized to x.
46 template <typename T>
47 struct summary_stats_unary_op
48 {
49 __host__ __device__
operator ()summary_stats_unary_op50 summary_stats_data<T> operator()(const T& x) const
51 {
52 summary_stats_data<T> result;
53 result.n = 1;
54 result.min = x;
55 result.max = x;
56 result.mean = x;
57 result.M2 = 0;
58 result.M3 = 0;
59 result.M4 = 0;
60
61 return result;
62 }
63 };
64
65 // summary_stats_binary_op is a functor that accepts two summary_stats_data
66 // structs and returns a new summary_stats_data which are an
67 // approximation to the summary_stats for
68 // all values that have been agregated so far
69 template <typename T>
70 struct summary_stats_binary_op
71 : public thrust::binary_function<const summary_stats_data<T>&,
72 const summary_stats_data<T>&,
73 summary_stats_data<T> >
74 {
75 __host__ __device__
operator ()summary_stats_binary_op76 summary_stats_data<T> operator()(const summary_stats_data<T>& x, const summary_stats_data <T>& y) const
77 {
78 summary_stats_data<T> result;
79
80 // precompute some common subexpressions
81 T n = x.n + y.n;
82 T n2 = n * n;
83 T n3 = n2 * n;
84
85 T delta = y.mean - x.mean;
86 T delta2 = delta * delta;
87 T delta3 = delta2 * delta;
88 T delta4 = delta3 * delta;
89
90 //Basic number of samples (n), min, and max
91 result.n = n;
92 result.min = thrust::min(x.min, y.min);
93 result.max = thrust::max(x.max, y.max);
94
95 result.mean = x.mean + delta * y.n / n;
96
97 result.M2 = x.M2 + y.M2;
98 result.M2 += delta2 * x.n * y.n / n;
99
100 result.M3 = x.M3 + y.M3;
101 result.M3 += delta3 * x.n * y.n * (x.n - y.n) / n2;
102 result.M3 += (T) 3.0 * delta * (x.n * y.M2 - y.n * x.M2) / n;
103
104 result.M4 = x.M4 + y.M4;
105 result.M4 += delta4 * x.n * y.n * (x.n * x.n - x.n * y.n + y.n * y.n) / n3;
106 result.M4 += (T) 6.0 * delta2 * (x.n * x.n * y.M2 + y.n * y.n * x.M2) / n2;
107 result.M4 += (T) 4.0 * delta * (x.n * y.M3 - y.n * x.M3) / n;
108
109 return result;
110 }
111 };
112
113 template <typename Iterator>
print_range(const std::string & name,Iterator first,Iterator last)114 void print_range(const std::string& name, Iterator first, Iterator last)
115 {
116 typedef typename std::iterator_traits<Iterator>::value_type T;
117
118 std::cout << name << ": ";
119 thrust::copy(first, last, std::ostream_iterator<T>(std::cout, " "));
120 std::cout << "\n";
121 }
122
123
main(void)124 int main(void)
125 {
126 typedef float T;
127
128 // initialize host array
129 T h_x[] = {4, 7, 13, 16};
130
131 // transfer to device
132 thrust::device_vector<T> d_x(h_x, h_x + sizeof(h_x) / sizeof(T));
133
134 // setup arguments
135 summary_stats_unary_op<T> unary_op;
136 summary_stats_binary_op<T> binary_op;
137 summary_stats_data<T> init;
138
139 init.initialize();
140
141 // compute summary statistics
142 summary_stats_data<T> result = thrust::transform_reduce(d_x.begin(), d_x.end(), unary_op, init, binary_op);
143
144 std::cout <<"******Summary Statistics Example*****"<<std::endl;
145 print_range("The data", d_x.begin(), d_x.end());
146
147 std::cout <<"Count : "<< result.n << std::endl;
148 std::cout <<"Minimum : "<< result.min <<std::endl;
149 std::cout <<"Maximum : "<< result.max <<std::endl;
150 std::cout <<"Mean : "<< result.mean << std::endl;
151 std::cout <<"Variance : "<< result.variance() << std::endl;
152 std::cout <<"Standard Deviation : "<< std::sqrt(result.variance_n()) << std::endl;
153 std::cout <<"Skewness : "<< result.skewness() << std::endl;
154 std::cout <<"Kurtosis : "<< result.kurtosis() << std::endl;
155
156 return 0;
157 }
158
159