Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
DdtrTester.cpp
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 #include "cheetah/ddtr/test_utils/DdtrTester.h"
25 #include "cheetah/ddtr/Ddtr.h"
26 #include "cheetah/generators/GaussianNoiseConfig.h"
27 #include "cheetah/generators/GaussianNoise.h"
28 #include "cheetah/generators/PulsarInjection.h"
29 #include "cheetah/generators/PulsarInjectionConfig.h"
30 #include "cheetah/data/TimeFrequency.h"
31 #include "cheetah/sigproc/SigProcWriter.h"
32 #include "unistd.h"
33 #include <thread>
34 #include <chrono>
35 #include <sstream>
36 #include <cmath>
37 
38 namespace ska {
39 namespace cheetah {
40 namespace ddtr {
41 namespace test {
42 
43 
44 template <typename TestTraits>
45 DdtrTester<TestTraits>::DdtrTester()
46  : cheetah::utils::test::AlgorithmTester<TestTraits>()
47 {
48 }
49 
50 template <typename TestTraits>
51 DdtrTester<TestTraits>::~DdtrTester()
52 {
53 }
54 
55 template<typename TestTraits>
56 void DdtrTester<TestTraits>::SetUp()
57 {
58 }
59 
60 template<typename TestTraits>
61 void DdtrTester<TestTraits>::TearDown()
62 {
63 }
64 
65 template<typename ArchitectureTag, typename NumericalT, typename TimeFreqT>
66 DdtrTesterTraits<ArchitectureTag, NumericalT, TimeFreqT>::DdtrTesterTraits()
67  : _dm_called(false)
68  , _dm_call_count(0u)
69 {
70 }
71 
72 template<typename ArchitectureTag, typename NumericalT, typename TimeFreqT>
73 typename DdtrTesterTraits<ArchitectureTag, NumericalT, TimeFreqT>::Api& DdtrTesterTraits<ArchitectureTag, NumericalT, TimeFreqT>::api(PoolType& pool)
74 {
75  typedef ddtr::DedispersionConfig::Dm Dm;
76  if(!_api) {
77  configure(_config); // call configuration method
78  _config.add_dm_range(Dm(0.0 * data::parsecs_per_cube_cm),Dm(1000.0 * data::parsecs_per_cube_cm), Dm(100.0 * data::parsecs_per_cube_cm));
79  _config.set_pool(pool);
80  _api.reset(new Api(_config,
81  [this](std::shared_ptr<DmType> data)
82  {
83  PANDA_LOG << "DM handler called";
84  _dm_called = true;
85  _dm_data.push_back(data);
86  ++_dm_call_count;
87  PANDA_LOG << "dm_call_count:"<<_dm_call_count;
88  } ));
89  }
90  return *_api;
91 }
92 
93 template<typename ArchitectureTag, typename NumericalT, typename TimeFreqT>
94 typename DdtrTesterTraits<ArchitectureTag, NumericalT, TimeFreqT>::DmDataContainerType const& DdtrTesterTraits<ArchitectureTag, NumericalT, TimeFreqT>::dm_data() const
95 {
96  return _dm_data;
97 }
98 
99 template <class T>
100 void dump_host_buffer(T* buffer,size_t size, std::string filename)
101 {
102  std::ofstream outfile;
103  outfile.open(filename.c_str(),std::ifstream::out | std::ifstream::binary);
104  outfile.write((char*)buffer, size*sizeof(T));
105  outfile.close();
106 }
107 
108 template<typename ArchitectureTag, typename NumericalT, typename TimeFreqT>
109 ddtr::Config& DdtrTesterTraits<ArchitectureTag, NumericalT, TimeFreqT>::config()
110 {
111  return _config;
112 }
113 
114 template<typename ArchitectureTag, typename NumericalT, typename TimeFreqT>
115 bool DdtrTesterTraits<ArchitectureTag, NumericalT, TimeFreqT>::dm_handler_called() const
116 {
117  return _dm_called;
118 }
119 
120 
121 template<typename ArchitectureTag, typename NumericalT, typename TimeFreqT>
123 {
124  return _dm_call_count;
125 }
126 
127 POOL_ALGORITHM_TYPED_TEST_P(DdtrTester, test_handlers)
128 {
135  TypeParam traits;
136  typedef typename TypeParam::NumericalRep NumericalRep;
137  auto& api = traits.api(pool);
138 
139  // generate random noise data and pass to Ddtr
140  // until we get a callback
141  generators::GaussianNoiseConfig generator_config;
142  generators::GaussianNoise<NumericalRep> generator(generator_config);
143 
144  typedef typename TypeParam::TimeFrequencyType DataType;
145 
146  // keep sending data until we get a callback
147  std::size_t count=0;
148  std::size_t chunk_samples=10000;
149  for(int ii=0; ii < 4; ++ii)
150  {
151  std::shared_ptr<DataType> data = DataType::make_shared(pss::astrotypes::DimensionSize<pss::astrotypes::units::Time>(chunk_samples),
152  pss::astrotypes::DimensionSize<pss::astrotypes::units::Frequency>(10));
153  data->sample_interval(typename DataType::TimeType(150.0 * boost::units::si::milli * boost::units::si::seconds));
154  data->set_channel_frequencies_const_width(data::FrequencyType(1000.0 * boost::units::si::mega * data::hz), data::FrequencyType(-30.0 * boost::units::si::mega * data::hz));
155  PANDA_LOG_DEBUG << "Populating test vector with generator data";
156  generator.next(*data);
157  api(*data);
158  ++count;
159  PANDA_LOG_DEBUG << "In loop on count " << count;
160  }
161 
162  while (!traits.dm_handler_called())
163  {
164  std::this_thread::sleep_for(std::chrono::seconds(1));
165  }
166  ASSERT_TRUE(traits.dm_handler_called());
167 
168  // now we expect the buffer to flush data as it is destroyed. so the pool has to stay alive whilst this happns
169 }
170 
171 
172 POOL_ALGORITHM_TYPED_TEST_P(DdtrTester, test_algorithm)
173 {
182  TypeParam traits;
183  auto& api = traits.api(pool);
184 
185  // generate Pulsar and pass to Ddtr
186  // until we get a callback
187 
188 
189  typedef typename TypeParam::TimeFrequencyType DataType;
190  // keep sending data until we get a callback
191 
192  std::size_t count=0;
193  std::size_t chunk_samples=10000;
194  std::size_t no_of_chans=10;
195  for (int ii=0; ii<4; ++ii)
196  {
197  std::shared_ptr<DataType> data = DataType::make_shared(pss::astrotypes::DimensionSize<pss::astrotypes::units::Time>(chunk_samples),
198  pss::astrotypes::DimensionSize<pss::astrotypes::units::Frequency>(no_of_chans));
199  data->sample_interval(typename DataType::TimeType(150.0 * boost::units::si::milli * boost::units::si::seconds));
200  data->set_channel_frequencies_const_width(data::FrequencyType(1000.0 * boost::units::si::mega * data::hz), data::FrequencyType(-30.0 * boost::units::si::mega * data::hz));
201  PANDA_LOG_DEBUG << "Populating test vector with dispersed pulsed data";
202 
203  //instantiate data with zeros
204  std::fill(data->begin(), data->end(), 0);
205 
206  // Calculate delays and accordingly add 3 bin wide pulse at those positions
207  float dm_constant = 4.1493775933609e3;
208  std::vector<data::FrequencyType> const& channel_freqs = data->channel_frequencies();
209  data::FrequencyType freq_top = *std::max_element(channel_freqs.begin(), channel_freqs.end());
210  std::vector<float> delays;
211  for (auto freq: channel_freqs)
212  {
213  delays.push_back(dm_constant * 500.0 * (1.0/(freq.value()*freq.value()) - 1.0/(freq_top.value()*freq_top.value())) / data->sample_interval().value());
214  }
215 
216  for(std::size_t ii=0; ii<no_of_chans; ++ii)
217  {
218  typename DataType::Channel ts = data->channel(ii);
219  std::size_t delay = std::size_t(round(delays[ii]));
220  ts[pss::astrotypes::DimensionIndex<pss::astrotypes::units::Time>(delay)] = 1;
221  ts[pss::astrotypes::DimensionIndex<pss::astrotypes::units::Time>(delay + 1)] = 1;
222  ts[pss::astrotypes::DimensionIndex<pss::astrotypes::units::Time>(delay + 2)] = 1;
223  }
224 
225  PANDA_LOG_DEBUG << "Calling API";
226  api(*data);
227  ++count;
228  PANDA_LOG << "In loop on count " << count;
229  }
230 
231  while(!traits.dm_handler_called())
232  {
233  std::this_thread::sleep_for(std::chrono::seconds(1));
234  }
235  sync();
236  pool.wait();
237  auto const& dm_data = traits.dm_data();
238  auto const& dm_trials = (*dm_data[0]);
239  std::vector<float> max_values;
240  std::vector<float> trial_nos;
241  std::size_t cnt=0;
242 
243  for (auto it=dm_trials.cbegin();it!=dm_trials.cend();++it)
244  {
245  std::size_t max_value = *std::max_element(it->cbegin(),it->cend());
246  trial_nos.push_back(cnt);
247  max_values.push_back(max_value);
248  ++cnt;
249  }
250  ASSERT_TRUE(*std::max_element(max_values.begin(),max_values.end())==10);
251  ASSERT_TRUE(std::distance(max_values.begin(), std::max_element(max_values.begin(),max_values.end()))==5);
252 }
253 
254 POOL_ALGORITHM_TYPED_TEST_P(DdtrTester, test_pulse_with_noise)
255 {
264  TypeParam traits;
265  typedef typename TypeParam::NumericalRep NumericalRep;
266  auto& api = traits.api(pool);
267 
268  // generate Pulsar and pass to Ddtr
269  // until we get a callback
270  generators::GaussianNoiseConfig generator_config;
271  generators::GaussianNoise<NumericalRep> generator(generator_config);
272 
273  typedef typename TypeParam::TimeFrequencyType DataType;
274 
275  // keep sending data until we get a callback
276 
277  std::size_t count=0;
278  std::size_t chunk_samples=10000;
279  std::size_t no_of_chans=10;
280  for (int ii=0; ii<4; ++ii)
281  {
282  std::shared_ptr<DataType> data = DataType::make_shared(pss::astrotypes::DimensionSize<pss::astrotypes::units::Time>(chunk_samples),
283  pss::astrotypes::DimensionSize<pss::astrotypes::units::Frequency>(no_of_chans));
284  data->sample_interval(typename DataType::TimeType(50.0 * boost::units::si::milli * boost::units::si::seconds));
285  data->set_channel_frequencies_const_width(data::FrequencyType(1000.0 * boost::units::si::mega * data::hz), data::FrequencyType(-30.0 * boost::units::si::mega * data::hz));
286  PANDA_LOG_DEBUG << "Populating test vector with dispersed pulsed data";
287  //instantiate data with noise
288  generator.next(*data);
289  // Calculate delays and accordingly add 3 bin wide pulse at those positions
290  float dm_constant = 4.1493775933609e3;
291  std::vector<data::FrequencyType> const& channel_freqs = data->channel_frequencies();
292  data::FrequencyType freq_top = *std::max_element(channel_freqs.begin(), channel_freqs.end());
293  std::vector<float> delays;
294  for (auto freq: channel_freqs)
295  {
296  delays.push_back(dm_constant * 500.0 * (1.0/(freq.value()*freq.value()) - 1.0/(freq_top.value()*freq_top.value())) / data->sample_interval().value());
297  }
298 
299  for(std::size_t ii=0;ii<no_of_chans;++ii)
300  {
301  typename DataType::Channel ts = data->channel(ii);
302  std::size_t delay = std::size_t(round(delays[ii]));
303  ts[pss::astrotypes::DimensionIndex<pss::astrotypes::units::Time>(delay)] += 100; //Approx 4 sigma since sigma=24 ;
304  ts[pss::astrotypes::DimensionIndex<pss::astrotypes::units::Time>(delay + 1)] += 100;
305  ts[pss::astrotypes::DimensionIndex<pss::astrotypes::units::Time>(delay + 2)] += 100;
306  }
307 
308  PANDA_LOG_DEBUG << "Calling API";
309  api(*data);
310  ++count;
311  PANDA_LOG << "In loop on count " << count;
312  }
313 
314  while(!traits.dm_handler_called())
315  {
316  std::this_thread::sleep_for(std::chrono::seconds(1));
317  }
318  sync();
319  pool.wait();
320  auto const& dm_data = traits.dm_data();
321  auto const& dm_trials = (*dm_data[0]);
322  std::vector<float> max_values;
323  std::vector<float> trial_nos;
324  std::size_t cnt=0;
325 
326  for (auto it=dm_trials.cbegin(); it!=dm_trials.cend(); ++it)
327  {
328  std::size_t max_value = *std::max_element(it->cbegin(),it->cend());
329  trial_nos.push_back(cnt);
330  max_values.push_back(max_value);
331  ++cnt;
332  }
333  ASSERT_TRUE(std::distance(max_values.begin(), std::max_element(max_values.begin(),max_values.end()))==5);
334 
335 }
336 
337 // each test defined by ALGORITHM_TYPED_TEST_P must be added to the
338 // test register (each one as an element of the comma seperated list)
339 REGISTER_TYPED_TEST_CASE_P(DdtrTester, test_handlers, test_algorithm, test_pulse_with_noise);
340 
341 } // namespace test
342 } // namespace ddtr
343 } // namespace cheetah
344 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35
std::size_t dm_handler_call_count() const
return true if the dm_handler has been called
Definition: DdtrTester.cpp:122