Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
SpCcl.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/data/SpCcl.h"
25 #include "cheetah/utils/ModifiedJulianClock.h"
26 #include "pss/astrotypes/units/Time.h"
27 #include <algorithm>
28 #include <limits>
29 
30 namespace ska {
31 namespace cheetah {
32 namespace data {
33 
34 static const boost::units::quantity<MilliSeconds, double> dm_const_factor = (4.15 * 1000000.0 * boost::units::si::milli * boost::units::si::seconds);
35 
36 template<typename NumericalRep>
37 SpCcl<NumericalRep>::SpCcl()
38  : _dm_range(Dm(std::numeric_limits<typename Dm::value_type>::max() * parsecs_per_cube_cm),
39  Dm(std::numeric_limits<typename Dm::value_type>::min() * parsecs_per_cube_cm))
40 {
41 }
42 
43 template<typename NumericalRep>
44 SpCcl<NumericalRep>::SpCcl(SpCcl<NumericalRep>::BlocksType const& blocks, data::DimensionIndex<data::Time> offset)
45  : _dm_range(Dm(std::numeric_limits<typename Dm::value_type>::max() * parsecs_per_cube_cm),
46  Dm(std::numeric_limits<typename Dm::value_type>::min() * parsecs_per_cube_cm))
47  , _blocks(blocks)
48 {
49  if(!blocks.empty()) {
50  auto const& first_block = _blocks.front();
51  _offset_time = pss::astrotypes::units::duration_cast<boost::units::quantity<data::MilliSeconds,double>>((double)offset * first_block->sample_interval());
52  _start_time=first_block->start_time() + _offset_time;
53  std::pair<typename TimeFrequencyType::FrequencyType, typename TimeFrequencyType::FrequencyType> freq = first_block->low_high_frequencies();
54  _dm_delay_factor = ((1.0/(freq.first * freq.first)) - (1.0/(freq.second * freq.second)));
55  }
56 }
57 
58 template<typename NumericalRep>
59 SpCcl<NumericalRep>::~SpCcl()
60 {
61 }
62 
63 template<typename NumericalRep>
64 typename SpCcl<NumericalRep>::BlocksType const& SpCcl<NumericalRep>::tf_blocks() const
65 {
66  return _blocks;
67 }
68 
69 template<typename NumericalRep>
70 std::pair<typename SpCcl<NumericalRep>::Dm, typename SpCcl<NumericalRep>::Dm> const& SpCcl<NumericalRep>::dm_range() const
71 {
72  return _dm_range;
73 }
74 
75 template<typename NumericalRep>
77 {
78  if(dm > _dm_range.second) _dm_range.second = dm;
79  if(dm < _dm_range.first) _dm_range.first = dm;
80 }
81 
82 template<typename NumericalRep>
84 {
85  dm_max_min(cand.dm());
86  auto start_time = cand.tstart();
87  auto it = std::find_if(BaseT::rbegin()
88  , BaseT::rend()
89  , [&](SpCandidateType const& cand)
90  {
91  return start_time >= cand.tstart();
92  });
93  auto r_it = BaseT::insert(it.base(), cand);
94  r_it->tend(start_time + dm_const_factor * (cand.dm() * _dm_delay_factor).value());
95 }
96 
97 template<typename NumericalRep>
99 {
100  dm_max_min(cand.dm());
101  auto start_time = cand.tstart();
102  auto it = std::find_if(BaseT::rbegin()
103  , BaseT::rend()
104  , [&](SpCandidateType const& cand)
105  {
106  return start_time >= cand.tstart();
107  });
108  BaseT::insert(it.base(), cand);
109 }
110 
111 template<typename NumericalRep>
113 {
114  dm_max_min(cand.dm());
115  if(BaseT::empty() || cand.tstart() >= BaseT::back().tstart()) {
116  BaseT::emplace_back(std::move(cand));
117  }
118  else {
119  auto start_time = cand.tstart();
120  auto it = std::lower_bound(BaseT::begin()
121  , --BaseT::end()
122  , start_time
123  , [&](SpCandidateType const& cand, boost::units::quantity<MilliSeconds, double> start_time)
124  {
125  return start_time > cand.tstart();
126  });
127  BaseT::insert(it,std::move(cand));
128  }
129 }
130 
131 template<typename NumericalRep>
133 {
134  dm_max_min(cand.dm());
135  auto start_time = cand.tstart();
136  auto it = std::lower_bound(BaseT::begin()
137  , BaseT::end()
138  , start_time
139  , [&](SpCandidateType const& cand, boost::units::quantity<MilliSeconds, double> start_time)
140  {
141  return start_time > cand.tstart();
142  });
143 
144  BaseT::insert(it, std::move(cand));
145 }
146 
147 template<typename NumericalRep>
149 {
150  dm_max_min(cand.dm());
151  auto start_time = cand.tstart();
152  cand.tend(start_time + dm_const_factor * (_dm_delay_factor * cand.dm()).value());
153  auto it = std::lower_bound(BaseT::begin()
154  , BaseT::end()
155  , start_time
156  , [&](SpCandidateType const& cand, boost::units::quantity<MilliSeconds, double> start_time)
157  {
158  return start_time > cand.tstart();
159  });
160  BaseT::insert(it,std::move(cand));
161 }
162 
163 template<typename NumericalRep>
164 void SpCcl<NumericalRep>::emplace(typename SpCandidateType::Dm dm,
165  boost::units::quantity<MilliSeconds, double> start_time,
166  boost::units::quantity<MilliSeconds, double> pulse_width,
167  float sigma,
168  std::size_t ident)
169 {
170  dm_max_min(dm);
171  auto it = std::lower_bound(BaseT::begin()
172  , BaseT::end()
173  , start_time
174  , [&](SpCandidateType const& cand, boost::units::quantity<MilliSeconds, double> start_time)
175  {
176  return start_time > cand.tstart();
177  });
178  BaseT::insert(it,SpCandidateType(dm, start_time, pulse_width
179  ,start_time + dm_const_factor * (_dm_delay_factor * dm).value()
180  ,sigma, ident));
181 }
182 
183 template<typename NumericalRep>
184 void SpCcl<NumericalRep>::emplace_back(typename SpCandidateType::Dm dm,
185  boost::units::quantity<MilliSeconds, double> start_time,
186  boost::units::quantity<MilliSeconds, double> pulse_width,
187  float sigma,
188  std::size_t ident)
189 {
190  dm_max_min(dm);
191  BaseT::emplace_back(dm, start_time, pulse_width
192  ,start_time + dm_const_factor * (_dm_delay_factor * dm).value()
193  ,sigma, ident);
194 }
195 
196 template<typename NumericalRep>
198 {
199  if (BaseT::begin() == BaseT::end())
200  {
201  return ConstDataIterator(CandidateWindow(), BaseT::begin(), BaseT::end(), _blocks, _start_time);
202  }
203  else
204  {
205  return ConstDataIterator(window, BaseT::begin(), BaseT::end(), _blocks, _start_time);
206  }
207 }
208 
209 template<typename NumericalRep>
211 {
212  return ConstDataIterator(CandidateWindow(), BaseT::end(), BaseT::end(), _blocks, _start_time);
213 }
214 
215 template<typename NumericalRep>
216 utils::ModifiedJulianClock::time_point const& SpCcl<NumericalRep>::start_time() const
217 {
218  return _start_time;
219 }
220 
221 template<typename NumericalRep>
222 utils::ModifiedJulianClock::time_point SpCcl<NumericalRep>::start_time(SpCandidateType const& cand) const
223 {
224  return _start_time + pss::astrotypes::units::duration_cast<utils::julian_day>(cand.tstart());
225 }
226 
227 template<typename NumericalRep>
228 boost::units::quantity<data::MilliSeconds,double> SpCcl<NumericalRep>::offset_time() const
229 {
230  return _offset_time;
231 }
232 
233 template<typename NumericalRep>
234 template<typename PredicateT>
235 void SpCcl<NumericalRep>::remove_if(PredicateT const& pred)
236 {
237  Dm low = _dm_range.second;
238  Dm high = _dm_range.first;
239  auto it = std::remove_if(begin(), end(),
240  [&pred, &low, &high](SpCandidateType const& cand) {
241  bool v = pred(cand);
242  if(!v) {
243  if(cand.dm() > high) high = cand.dm();
244  if(cand.dm() < low) low = cand.dm();
245  }
246  return v;
247  });
248 
249  if (it != BaseT::end())
250  {
251  erase(it, BaseT::end());
252  }
253  _dm_range.first = low;
254  _dm_range.second = high;
255 }
256 
258 template<typename NumericalRep>
260  , ConstCandidateIterator cand_it_begin
261  , ConstCandidateIterator const& cand_it_end
262  , SpCcl::BlocksType const& tf_blocks
263  , utils::ModifiedJulianClock::time_point const& epoch
264  )
265  : _window(window)
266  , _cand_it(cand_it_begin)
267  , _cand_it_end(cand_it_end)
268  , _tf_it(tf_blocks.begin())
269  , _tf_end(tf_blocks.end())
270  , _tf_start_it(tf_blocks.begin())
271  , _epoch(epoch)
272 {
273  if (_cand_it != _cand_it_end)
274  {
275  _span = ((*_tf_it)->number_of_spectra() * static_cast<boost::units::quantity<MilliSeconds,double>>((*_tf_it)->sample_interval()));
276  _duration = _span;
277  this->next_candidate();
278  }
279 }
280 
281 template<typename NumericalRep>
283 {
284  if( _cand_it == _cand_it_end )
285  {
286  return;
287  }
288  auto increment_block = [this]()
289  {
290  ++_tf_it;
291  if (_tf_it == _tf_end)
292  {
293  return false;
294  }
295  _span = ((*_tf_it)->number_of_spectra() * static_cast<boost::units::quantity<MilliSeconds,double>>((*_tf_it)->sample_interval()));
296  _duration += _span;
297  return true;
298  };
299 
300  boost::units::quantity<MilliSeconds, double> start_time = _cand_it->tstart() - _window.ms_before();
301 
302  if (start_time < _duration - _span)
303  {
304  start_time = _duration - _span;
305  }
306 
307  auto end_time = _cand_it->tend() + _window.ms_after();
308  auto start_cand = _cand_it;
309 
310  // determine overlapping candidates
311  while(++_cand_it != _cand_it_end)
312  {
313  if(_cand_it->tstart() - _window.ms_before() <= end_time) {
314  auto next_end_time = _cand_it->tend() + _window.ms_after();
315  if(next_end_time > end_time) end_time = next_end_time;
316  }
317  else {
318  break;
319  }
320  }
321  // _cand_it is now in another data block
322  _cands.emplace_back(start_cand, _cand_it); // mark which candidate overlap taking into account windowing
323 
324  while(start_time >= _duration)
325  {
326  if (!increment_block())
327  {
328  throw panda::Error("Candidate Start time larger than end of data!!");
329  }
330  }
331  // set start time of the candidate block in MJD
332  _start_time = _epoch + pss::astrotypes::units::duration_cast<utils::julian_day>(start_time);
333 
334  data::DimensionIndex<data::Time> start_offset((*_tf_it)->number_of_spectra() - (std::size_t)(((_duration - start_time)/static_cast<boost::units::quantity<MilliSeconds,double>>((*_tf_it)->sample_interval()))));
335 
336  PANDA_LOG_DEBUG << "Start offset:" << start_offset;
337  while (end_time >= _duration)
338  {
339  data::DimensionIndex<data::Time> end_offset((std::size_t)((*_tf_it)->number_of_spectra()-1));
340  PANDA_LOG_DEBUG << "End Offset:" << end_offset;
341  _slices.push_back((*_tf_it)->slice(data::DimensionSpan<data::Time>(start_offset, end_offset)));
342  start_offset = data::DimensionIndex<data::Time>(0);
343  if (!increment_block())
344  {
345  _slices.back()._last_block = true;
346  return;
347  }
348  }
349  data::DimensionIndex<data::Time> end_offset((*_tf_it)->number_of_spectra() - 1 - (std::size_t)(((_duration - end_time)/static_cast<boost::units::quantity<MilliSeconds,double>>((*_tf_it)->sample_interval()))));
350  PANDA_LOG_DEBUG << "Start Offset=" << start_offset << " End Offset:" << end_offset;
351  _slices.emplace_back((*_tf_it)->slice(data::DimensionSpan<data::Time>(start_offset,end_offset)), true);
352 
353 }
354 
355 template<typename NumericalRep>
357 {
358  return _slices.front().last_block();
359 }
360 
361 template<typename NumericalRep>
362 utils::ModifiedJulianClock::time_point const& SpCcl<NumericalRep>::ConstDataIterator::start_time() const
363 {
364  return _start_time;
365 }
366 
367 template<typename NumericalRep>
369 {
370  auto& slice = _slices.front();
371  _start_time += pss::astrotypes::units::duration_cast<utils::julian_day>(slice.slice().template dimension<data::Time>() * (*_tf_start_it)->sample_interval());
372 
373  if(slice.last_block()) {
374  // keep candidate info up to date
375  _cands.pop_front();
376  }
377  _slices.pop_front();
378  if(_slices.size() == 0)
379  {
380  next_candidate();
381  }
382  return *this;
383 }
384 
385 template<typename NumericalRep>
386 typename SpCcl<NumericalRep>::ConstCandidateIterator SpCcl<NumericalRep>::ConstDataIterator::candidate_begin() const
387 {
388  return _cands.front().first;
389 }
390 
391 template<typename NumericalRep>
392 typename SpCcl<NumericalRep>::ConstCandidateIterator const& SpCcl<NumericalRep>::ConstDataIterator::candidate_end() const
393 {
394  return _cands.front().second;
395 }
396 
397 template<typename NumericalRep>
399 {
400  ConstDataIterator tmp(*this);
401  ++*this;
402  return tmp;
403 }
404 
405 template<typename NumericalRep>
406 typename SpCcl<NumericalRep>::ConstDataIterator::const_reference SpCcl<NumericalRep>::ConstDataIterator::operator*() const
407 {
408  return _slices.front().slice();
409 }
410 
411 template<typename NumericalRep>
413 {
414  if(d._slices.size() != this->_slices.size())
415  {
416  return false;
417  }
418 
419  if(d._slices.size() == 0) return true;
420 
421  return d._slices.front().slice() == this->_slices.front().slice();
422 }
423 
424 template<typename NumericalRep>
426 {
427  return !(d == *this);
428 }
429 
430 } // namespace data
431 } // namespace cheetah
432 } // namespace ska
ValueType const & front() const
the first emelment
Definition: VectorLike.cpp:232
ConstDataIterator data_begin(CandidateWindow const &window) const
Extract data corresponding to the candidates start and end times.
Definition: SpCcl.cpp:197
A simple record to hold &#39;candidate&#39; proprerties.
Definition: SpCandidate.h:58
void push_back(SpCandidateType const &cand)
Push back method for inserting Single pulse Candidates but without calculating tend They will be ins...
Definition: SpCcl.cpp:98
void push_back_calculate_duration(SpCandidateType const &cand)
Push back method for inserting Single pulse Candidates calculating temd based on dispersion and freq ...
Definition: SpCcl.cpp:83
std::pair< Dm, Dm > const & dm_range() const
the max and minimum Dispersion Measure values of the candidates
Definition: SpCcl.cpp:70
BlocksType const & tf_blocks() const
Access the TimeFrequency blocks associated with the single pulse candidtes.
Definition: SpCcl.cpp:64
void emplace(SpCandidateType &&cand)
Emplace back method for inserting Single pulse Candidates, ensuring ordering by start_time. They will be inserted in a sorted order with respect to the start time of each candidate instead of insert to the end.
Definition: SpCcl.cpp:132
TimePointType const & start_time() const
Get the start time of the candidate.
bool end_block() const
retunrs true is the current block is the end of contiguous series
Definition: SpCcl.cpp:356
void emplace_back(SpCandidateType &&cand)
Add a Single pulse Candidate to the end of the cand list, without calculating the duration Checking ...
Definition: SpCcl.cpp:112
Dm const & dm() const
access a reference to dm.
Some limits and constants for FLDO.
Definition: Brdz.h:35
ConstCandidateIterator candidate_begin() const
retunrs true is the current block is the end of contiguous series
Definition: SpCcl.cpp:386
MsecTimeType const & tstart() const
Get the start time of the candidate.
Definition: SpCandidate.cpp:74
utils::ModifiedJulianClock::time_point const & start_time() const
the start time of the first relevant data (taking into account the offset)
Definition: SpCcl.cpp:216
SpCandidate list.
Definition: SpCcl.h:52
ConstDataIterator & operator++()
pre-increment operator
Definition: SpCcl.cpp:368
Dm const & dm() const
access a reference to dm.
define a window of data surrounding a single pulse candidate
boost::units::quantity< data::MilliSeconds, double > offset_time() const
the time duration from the beginning of the first block to start_time()
Definition: SpCcl.cpp:228
Sigma const & sigma() const
access a reference to standard deviation (sigma).
void remove_if(PredicateT const &)
remove candidates that meet the specified condition
Definition: SpCcl.cpp:235
utils::ModifiedJulianClock::time_point const & start_time() const
: Get the start time
Definition: SpCcl.cpp:362