Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
SummaryStatistics.cu
1 /*
2  * The MIT License (MIT)
3  *
4  * Copyright (c) 2016 The SKA organisation
5  *
6  * Permission is hereby granted, free of charge, to any person obtaining a copy
7  * of this software and associated documentation files (the "Software"), to deal
8  * in the Software without restriction, including without limitation the rights
9  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
10  * copies of the Software, and to permit persons to whom the Software is
11  * furnished to do so, subject to the following conditions:
12  *
13  * The above copyright notice and this permission notice shall be included in all
14  * copies or substantial portions of the Software.
15  *
16  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
17  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
18  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
19  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
20  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
21  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
22  * SOFTWARE.
23  */
24 
25 #include "cheetah/cuda_utils/cuda_thrust.h"
26 #include <cmath>
27 #include <limits>
28 
29 //
30 // This example computes several statistical properties of a data
31 // series in a single reduction. The algorithm is described in detail here:
32 // http://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
33 //
34 // Adapted from Thrust example code
35 //
36 
37 // structure used to accumulate the moments and other
38 // statistical properties encountered so far.
39 namespace ska {
40 namespace cheetah {
41 namespace fldo {
42 namespace cuda {
43 namespace util {
44 
45 template <typename T>
47 {
48  T n;
49  T mean;
50  T M2;
51 
52  // initialize to the identity element
53  void initialize()
54  {
55  n = mean = M2 = 0;
56  }
57 
58  T variance() { return M2 / (n - 1); }
59  T variance_n() { return M2 / n; }
60 };
61 
62 // stats_unary_op is a functor that takes in a value x and
63 // returns a variace_data whose mean value is initialized to x.
64 template <typename T>
66 {
67  __host__ __device__
68  summary_stats_data<T> operator()(const unsigned char& x) const
69  {
70  summary_stats_data<T> result;
71  result.n = 1;
72  result.mean = (float)x - 128.;
73  result.M2 = 0;
74 
75  return result;
76  }
77 };
78 
79 template <typename T>
81 {
82  __host__ __device__
83  summary_stats_data<T> operator()(const T& x) const
84  {
85  summary_stats_data<T> result;
86  result.n = 1;
87  result.mean = (float)x ;
88  result.M2 = 0;
89 
90  return result;
91  }
92 };
93 
94 // summary_stats_binary_op is a functor that accepts two summary_stats_data
95 // structs and returns a new summary_stats_data which are an
96 // approximation to the summary_stats for
97 // all values that have been agregated so far
98 template <typename T>
100  : public thrust::binary_function<const summary_stats_data<T>&,
101  const summary_stats_data<T>&,
102  summary_stats_data<T> >
103 {
104  __host__ __device__
105  summary_stats_data<T> operator()(const summary_stats_data<T>& x, const summary_stats_data <T>& y) const
106  {
107  summary_stats_data<T> result;
108 
109  // precompute some common subexpressions
110  T n = x.n + y.n;
111  //T n2 = n * n;
112 
113  T delta = y.mean - x.mean;
114  //printf("y.mean: %f x.mean: %f x.n: %d\n", y.mean, x.mean, x.n);
115  T delta2 = delta * delta;
116 
117  //Basic number of samples n
118  result.n = n;
119  result.mean = x.mean + delta * y.n / n;
120 
121  result.M2 = x.M2 + y.M2;
122  result.M2 += delta2 * x.n * y.n / n;
123  return result;
124  }
125 };
126 
127 /*
128  * void statistics(unsigned char *raw_in, int N, float *mean, float *rms, cudaStream_t stream)
129  *
130  * @brief Calculates data statistics (mean and variance) of the input data
131  * (unsigned char).
132  *
133  * @param raw_in device pointer to input data
134  * @param N the number of data to process
135  * @param mean the calculcated mean value
136  * @param rms the calculcated rms value
137  * @param stream the cuda stream to run kernels
138  *
139  */
140 void statistics(unsigned char *raw_in, int N, float *mean, float *rms, cudaStream_t stream)
141 {
142  typedef float T;
143 
144  // initialize host array
145 
146  // transfer to device
147  thrust::cuda::par.on(stream);
148  thrust::device_ptr<unsigned char> dev_in = thrust::device_pointer_cast(raw_in);
149 
150  // setup arguments
151  summary_stats_unary_op<T> unary_op;
152  summary_stats_binary_op<T> binary_op;
154 
155  init.initialize();
156 
157  // compute summary statistics
158  summary_stats_data<T> result = thrust::transform_reduce(dev_in, dev_in + N, unary_op, init, binary_op);
159  *mean = result.mean;
160  *rms = sqrt(result.variance_n());
161  return;
162 }
163 
164 /*
165  * void statistics_float(float *raw_in, int N, float *mean, float *rms, cudaStream_t stream)
166  *
167  * @brief Calculates data statistics (mean and variance) of the input data
168  * (unsigned char).
169  *
170  * @param raw_in device pointer to input data
171  * @param N the number of data to process
172  * @param mean the calculcated mean value
173  * @param rms the calculcated rms value
174  * @param stream the cuda stream to run kernels
175  *
176  */
177 void statistics_float(float *raw_in, int N, float *mean, float *rms, cudaStream_t stream)
178 {
179  typedef float T;
180 
181  // initialize host array
182 
183  // transfer to device
184  thrust::cuda::par.on(stream);
185  thrust::device_ptr<T> dev_in = thrust::device_pointer_cast(raw_in);
186 
187  // setup arguments
189  summary_stats_binary_op<T> binary_op;
191 
192  init.initialize();
193 
194  // compute summary statistics
195  summary_stats_data<T> result = thrust::transform_reduce(dev_in, dev_in + N, unary_op, init, binary_op);
196  *mean = result.mean;
197  *rms = sqrt(result.variance_n());
198  return;
199 }
200 
201 void find_max(float *raw_in, int start, int N, float &sn_max, float &pulse_width, int &index)
202 {
203  typedef float T;
204  thrust::device_ptr<float> d_in = thrust::device_pointer_cast(raw_in);
205 
206  //thrust::device_vector<float>::iterator iter = thrust::max_element(d_in, d_in + N);
207  thrust::device_ptr<float> max_sn_pos;
208  thrust::device_ptr<float> width_pos;
209  max_sn_pos = thrust::max_element(d_in + start, d_in + start + N);
210 
211  index = max_sn_pos - (d_in + start);
212  sn_max = *max_sn_pos;
213  width_pos = d_in + start + index + N;
214  pulse_width = *width_pos;
215 }
216 } // util
217 } // namespace cuda
218 } // namespace fldo
219 } // namespace cheetah
220 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35