Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Sift.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/sift/simple_sift/Sift.h"
25 #include "panda/Log.h"
26 #include <algorithm>
27 
28 namespace ska {
29 namespace cheetah {
30 namespace sift {
31 namespace simple_sift {
32 
33 Sift::Sift(sift::Config const& algo_config)
34  : _config(algo_config.template config<simple_sift::Config>())
35 {
36 }
37 
38 // A comparator to sort a list of candidates by sigma, from highest to lowest
39 static bool sort_vector_by_sigma_asc(data::Ccl::CandidateType const& candidate_one, data::Ccl::CandidateType const& candidate_two) {
40  return (candidate_one.sigma() > candidate_two.sigma());
41 }
42 
43 std::shared_ptr<data::Scl> Sift::operator()(ResourceType& resource, std::shared_ptr<data::Ccl> const& ccl) const
44 {
45  return (*this)(resource, *ccl);
46 }
47 
48 std::shared_ptr<data::Scl> Sift::operator()(ResourceType&, data::Ccl& ccl) const
49 {
50  double candidate_period, test_period, candidate_dm, test_dm, period_ratio;
51  std::vector<bool> identified_as_harmonic(ccl.size(), {false});
52  std::vector<double> harmonic_ratio;
53  std::size_t candidate_index, test_candidate_index;
54  data::Ccl::CandidateType candidate, test_candidate;
55 
56  // Shared pointer to structure in which to store sifted periodicity candidates
57  auto scl = std::make_shared<data::Scl>();
58 
59  // Sort the list by S/N, highest to lowest
60  std::sort(ccl.begin(), ccl.end(), sort_vector_by_sigma_asc);
61 
62  // Calculate harmonic ratios (1/1, 2/2, 3/2, 4/2, ..., 8/7, 8/8)
63  for (std::size_t top = 1; top <= _config.num_candidate_harmonics(); top++) {
64  for (std::size_t bottom = 1; bottom <= top; bottom++) {
65  harmonic_ratio.push_back((double) top/bottom);
66  }
67  }
68 
69  // Sort the vector of ratios
70  std::sort(harmonic_ratio.begin(), harmonic_ratio.end());
71 
72  // Remove duplicate ratios and resize the vector
73  auto iter = unique(harmonic_ratio.begin(), harmonic_ratio.end());
74  harmonic_ratio.resize(distance(harmonic_ratio.begin(), iter));
75 
76  // Loop until all candidates have been dealt with
77  while (true) {
78 
79  // Find the candidate with the highest sigma that has not been sifted
80  candidate_index = std::find(identified_as_harmonic.begin(), identified_as_harmonic.end(), 0) - identified_as_harmonic.begin();
81  if (candidate_index >= identified_as_harmonic.size()) { // If the index is greater than the vector size, it means it has not found any zeros
82  break;
83  }
84 
85  PANDA_LOG_DEBUG << "candidate_index = " << candidate_index;
86 
87  // Start with the highest S/N candidate and check to see if any of the other candidates are related
88  candidate = ccl[candidate_index];
89  candidate_period = boost::units::quantity_cast<double>(candidate.period());
90  candidate_dm = boost::units::quantity_cast<double>(candidate.dm());
91 
92  PANDA_LOG_DEBUG << "candidate_period = " << candidate_period << ", candidate_dm = " << candidate_dm;
93 
94  // Put this candidate into the sifted candidate list
95  scl->push_back(candidate);
96 
97  // Mark that this candidate has been added to the scl
98  identified_as_harmonic[candidate_index] = 1;
99 
100  // Go through every remaining signal in the list to check for related signals
101  for (test_candidate_index = candidate_index + 1; test_candidate_index < ccl.size(); test_candidate_index++) {
102 
103  if (identified_as_harmonic[test_candidate_index]) {
104  continue;
105  } else {
106 
107  PANDA_LOG_DEBUG << "test_cand_index = " << test_candidate_index;
108 
109  test_candidate = ccl[test_candidate_index];
110  test_period = boost::units::quantity_cast<double>(test_candidate.period());
111  test_dm = boost::units::quantity_cast<double>(test_candidate.dm());
112 
113  PANDA_LOG_DEBUG << "test_period = " << test_period << ", test_dm = " << test_dm;
114 
115  // Since all of the ratios generated to test for harmonics are greater than one, make sure the ratio of the periods is greater than one
116  if (candidate_period >= test_period) {
117  period_ratio = candidate_period/test_period;
118  } else {
119  period_ratio = test_period/candidate_period;
120  }
121 
122  // Go through each ratio and test whether this period is some harmonic of the candidate
123  for (auto iter = harmonic_ratio.begin(); iter != harmonic_ratio.end(); ++iter) {
124  if (std::abs(period_ratio - *iter) < (_config.match_factor() * *iter)) {
125  // Test whether this DM is within +-10 of the DM of the candidate
126  if (std::abs(candidate_dm - test_dm) < 10) {
127  // Mark this test_candidate as a harmonic
128  identified_as_harmonic[test_candidate_index] = 2;
129  PANDA_LOG_DEBUG << "Harmonic found!";
130  break;
131  }
132  }
133  }
134 
135  }
136 
137  }
138 
139  }
140 
141  return scl;
142 
143 }
144 
145 } //simple_sift
146 } //sift
147 } //cheetah
148 } //ska
Sift(sift::Config const &config)
Construct a new Sift instance.
Definition: Sift.cpp:33
Configuration for the simple_sift module.
Definition: Config.h:38
std::shared_ptr< data::Scl > operator()(ResourceType &cpu, data::Ccl &ccl) const
Sift a candidate list to aggregate like signals.
Definition: Sift.cpp:48
std::size_t size() const
Retrieve the size of the underlying vector.
Definition: VectorLike.cpp:61
T const & sigma() const
access a reference to _sigma.
Definition: Candidate.cpp:111
std::size_t num_candidate_harmonics() const
The max number of harmonics to search.
Definition: Config.cpp:45
Some limits and constants for FLDO.
Definition: Brdz.h:35
Configuration for the sift module.
Definition: Config.h:39
Candidate list.
Definition: Ccl.h:45
void push_back(ValueType const &value)
Appends element to end of vector.
Definition: VectorLike.cpp:145
Iterator end()
An iterator pointing to the end of the vector.
Definition: VectorLike.cpp:109
A simple record to hold &#39;candidate&#39; proprerties.
Definition: Candidate.h:60
Dm const & dm() const
access a reference to dm.
Definition: Candidate.cpp:85
MsecTimeType const & period() const
access a reference to pulsar period.
Definition: Candidate.cpp:61
Iterator begin()
An iterator pointing to the start of the vector.
Definition: VectorLike.cpp:85