Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Dred.cu
1 #include "cheetah/dred/cuda/Dred.cuh"
2 #include "cheetah/cuda_utils/cuda_thrust.h"
3 #include "cheetah/cuda_utils/cuda_errorhandling.h"
4 #include "cheetah/data/Units.h"
5 
6 namespace ska {
7 namespace cheetah {
8 namespace dred {
9 namespace cuda {
10 
11 template <typename T>
12 Dred<T>::Dred(Config const& config, dred::Config const& algo_config)
13  : utils::AlgorithmBase<Config,dred::Config>(config, algo_config)
14  , _pwft(_algo_config.pwft_config())
15 {
16 }
17 
18 template <typename T>
20 {
21 }
22 
23 template <typename T>
25 {
26  size_t count = in.size();
27  const T* ptr = thrust::raw_pointer_cast(in.data());
28  if( count == 1 )
29  {
30  out[0] = in[0];
31  }
32  else if( count == 2 )
33  {
34  out[0] = 0.5f*(in[0] + in[1]);
35  }
36  else if( count == 3 )
37  {
38  out[0] = detail::median3<T>(in[0],in[1],in[2]);
39  }
40  else if( count == 4 )
41  {
42  out[0] = detail::median4<T>(in[0],in[1],in[2],in[3]);
43  }
44  else
45  {
46  thrust::transform(thrust::cuda::par,
47  thrust::make_counting_iterator<unsigned int>(0),
48  thrust::make_counting_iterator<unsigned int>(in.size()/5),
49  out.begin(), detail::Median5Functor<T>(ptr));
50  }
51 }
52 
53 template <typename T>
54 void Dred<T>::linear_stretch(const SeriesType& in, SeriesType& out, float step)
55 {
56  T const* ptr = thrust::raw_pointer_cast(in.data());
57  thrust::transform(thrust::cuda::par,
58  thrust::make_counting_iterator<unsigned int>(0),
59  thrust::make_counting_iterator<unsigned int>(out.size()),
60  out.begin(), detail::LinearStretchFunctor<T>(ptr, in.size(), step));
61 }
62 
63 template <typename T>
64 template <typename Alloc>
67  data::AccelerationType maximum_acceleration)
68 {
69  output.resize(input.size());
70  output.frequency_step(input.frequency_step());
71 
72  //observation length in seconds
73  float tobs = 1/input.frequency_step().value();
74  //nyquist frequency in hz
75  float nyquist = input.frequency_step().value() * input.size();
76  //factor by which median window should be bigger than the Z value at a given frequency
77  std::size_t k = _algo_config.oversmoothing_factor();
78  //The power of 5 to smooth by
79  std::size_t power = 1;
80 
81  //clear the vector of boundaries and windows
82  _boundaries.clear();
83 
84  //Get the acceleration value that will determine the smoothing windows
85  auto accel_max = maximum_acceleration.value();
86 
87  //calculate the boundaries
88  //this will be done for every execution but should incur only minor performance penalties
89  while (true) {
90  Boundary boundary;
91  boundary.window = std::pow((std::size_t)5,(std::size_t)power);
92  auto c = boost::units::si::constants::codata::c.value().value();
93  boundary.frequency = (c * boundary.window) / (k * accel_max * tobs * tobs);
94  boundary.idx = std::min((std::size_t) (boundary.frequency/input.frequency_step().value()), input.size());
95  _boundaries.push_back(boundary);
96  if (boundary.frequency>nyquist)
97  break;
98  power+=1;
99  }
100 
101  //This array will only be resized if the number of boundaries
102  //has changed since the last execution
103  _medians.resize(_boundaries.size());
104  for (int ii=0; ii<_boundaries.size(); ++ii)
105  {
106  //similarly these array will only be resized if either the input
107  //size or the actual boundaries have changed.
108  _medians[ii].resize(input.size()/_boundaries[ii].window);
109  }
110 }
111 
112 
113 template <typename T>
114 template <typename Alloc>
115 void Dred<T>::process(ResourceType& gpu,
118  data::AccelerationType maximum_acceleration)
119 {
120  prepare(input, output, maximum_acceleration);
121  SeriesType power_spectrum(input.size());
122  power_spectrum.frequency_step(input.frequency_step());
123  SeriesType intermediate(input.size());
124  SeriesType baseline(input.size());
125 
126  _pwft.process_direct(gpu,input, power_spectrum);
127  SeriesType* in = &(power_spectrum);
128  SeriesType* out;
129  std::size_t offset = 0;
130  for (int ii=0;ii<_boundaries.size();ii++)
131  {
132  out = &(_medians)[ii];
133  median_scrunch5(*in, *out);
134  linear_stretch(*out, intermediate, (float)(_boundaries[ii].window));
135  thrust::copy_n(intermediate.begin()+offset,
136  _boundaries[ii].idx-offset,
137  baseline.begin()+offset);
138  offset = _boundaries[ii].idx;
139  in = out;
140  }
141  thrust::transform(thrust::cuda::par, input.begin(),input.end(),
142  baseline.begin(),output.begin(),detail::PowerNormalisingFunctor<T>());
143  CUDA_ERROR_CHECK(cudaDeviceSynchronize());
144 }
145 
146 } //cuda
147 } //dred
148 } //cheetah
149 } //ska
Functor for normalising spectral data.
Definition: Normaliser.cuh:18
FourierFrequencyType const & frequency_step() const
Retrieve the frequency step of the series.
void process_direct(panda::PoolResource< Arch > &resource, data::FrequencySeries< Arch, typename data::ComplexTypeTraits< Arch, T >::type, InputAlloc >const &input, data::PowerSeries< Arch, T, OutputAlloc > &output, Args &&... args)
Form power spectrum using absolute squared.
Definition: Pwft.cpp:33
A container of Fourier series data.
CUDA/Thrust implementation of the Dred algorithm.
Definition: Dred.cuh:33
Configuration object for Dred CUDA implementation.
Definition: Config.h:43
Some limits and constants for FLDO.
Definition: Brdz.h:35
ConstIterator begin() const
Iterators to device memory.
Definition: Series.cpp:67
Dred(Config const &impl_config, dred::Config const &algo_config)
Create a new Dred instance.
Definition: Dred.cu:12
Class for power series (detected FrequencySeries).
Definition: PowerSeries.h:58
std::size_t oversmoothing_factor() const
The oversmoothing factor for median window calculation.
Definition: Config.cpp:61
Algorithm configuration for the Dred module.
Definition: Config.h:47
void process(panda::PoolResource< cheetah::Cuda > &resource, data::FrequencySeries< cheetah::Cuda, ComplexType, Alloc >const &input, data::FrequencySeries< cheetah::Cuda, ComplexType, Alloc > &output, data::AccelerationType maximum_acceleration)
Deredden a complex fourier spectrum.
void resize(std::size_t size)
resize the data
Definition: Series.cpp:115
std::size_t size() const
the size of the series
Definition: Series.cpp:109