Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Rfim.cpp
1 #include "cheetah/rfim/sum_threshold/Rfim.h"
2 #include "panda/Buffer.h"
3 #include "panda/Log.h"
4 #include "panda/TypeTraits.h"
5 #pragma GCC diagnostic push
6 #pragma GCC diagnostic ignored "-Wunused-parameter"
7 #include <boost/accumulators/accumulators.hpp>
8 #include <boost/accumulators/statistics.hpp>
9 #pragma GCC diagnostic pop
10 #include <cmath>
11 
12 
13 namespace ska {
14 namespace cheetah {
15 namespace rfim {
16 namespace sum_threshold {
17 
18 namespace {
19 typedef data::TimeFrequencyFlags<Cpu> FlagType;
20 
21 inline float calcThresholdI(float threshold1, unsigned window, float p)
22 {
23  if (p <= 0.0f) {
24  p = 1.5f; // according to Andre's RFI paper, this is a good default value
25  }
26  return (float) (threshold1 * std::pow(p, std::log2(window)) / window);
27 }
28 
29 template<typename IteratorType>
30 static void calculateNormalStatistics(
31  IteratorType&& begin, IteratorType const&& end, float& mean, float& median, float& std_dev)
32 {
33  boost::accumulators::accumulator_set<
34  double, boost::accumulators::stats<boost::accumulators::tag::variance,
35  boost::accumulators::tag::median>> acc;
36  for_each(begin, end, std::bind<void>(std::ref(acc), std::placeholders::_1));
37 
38  mean = boost::accumulators::mean(acc);
39  std_dev = sqrtf(boost::accumulators::variance(acc));
40  median = boost::accumulators::median(acc);
41 }
42 
43 template<typename IteratorType>
44 static void calculateStatistics(IteratorType&& begin, IteratorType const&& end, float& mean, float& median, float& std_dev) {
45  // In Rob's code, there is also an option for calculateWinsorizedStatistics. I left it out for now.
46  calculateNormalStatistics(std::forward<IteratorType>(begin), std::forward<IteratorType const>(end), mean, median, std_dev);
47 }
48 
49 
50 
51 template<typename DataType>
52 static unsigned sumThreshold2DHorizontal(DataType const & powers, FlagType &flags, const unsigned window, const float threshold) {
53  unsigned extra_flagged = 0;
54 
55  for(data::DimensionIndex<data::Frequency> channel(0); channel<powers.number_of_channels(); ++channel) {
56  for (unsigned base = 0; base + window < powers.number_of_spectra(); ++base) {
57  float sum = 0.0f;
58  for(data::DimensionIndex<data::Time> time(base); time < base + window; ++time) {
59  if (flags[channel][time]) { // If it was flagged in a previous iteration, replace sample with current threshold
60  sum += threshold;
61  } else {
62  sum += powers.spectrum(time)[channel];
63  }
64  }
65 
66  if (sum >= window * threshold) {
67  // flag all samples in the sequence!
68  for(data::DimensionIndex<data::Time> time(base); time < base + window; ++time) {
69  if(!flags[channel][time]) {
70  ++extra_flagged;
71  flags[channel][time] = true;
72  }
73  }
74  }
75  }
76  }
77 
78  return extra_flagged;
79 }
80 
81 template<typename DataType>
82 static unsigned sumThreshold2DVertical(DataType const & powers
83  , FlagType &flags
84  , const unsigned window
85  , const float threshold)
86 {
87  unsigned extra_flagged = 0;
88 
89  for(data::DimensionIndex<data::Frequency> channel(0); channel<powers.number_of_channels(); ++channel) {
90  for (unsigned base = 0; base + window < powers.number_of_spectra(); base++) {
91  float sum = 0.0f;
92 
93  for (data::DimensionIndex<data::Time> time(base); time < base + window; ++time) {
94  if(flags[channel][time]) { // If it was flagged in a previous iteration, replace sample with current threshold
95  sum += threshold;
96  } else {
97  sum += powers.spectrum(time)[channel];
98  }
99  }
100 
101  if (sum >= window * threshold) {
102  // flag all samples in the sequence!
103  for (data::DimensionIndex<data::Time> time(base); time < base + window; ++time) {
104  if(!flags[channel][time]) {
105  ++extra_flagged;
106  flags[channel][time] = true;
107  }
108  }
109  }
110  }
111  }
112 
113  return extra_flagged;
114 }
115 } // namespace
116 
117 template<typename RfimTraits>
118 Rfim<RfimTraits>::Rfim(const Config& config)
119  : _config(config)
120 {
121 }
122 
123 template<typename RfimTraits>
124 Rfim<RfimTraits>::~Rfim()
125 {
126 }
127 
128 namespace {
129  // Produces a suitable type for recording rfim flags
130  // default is to use a new structure
131  template<typename DataAdapterType, typename Enable=void>
132  struct FlagTypeHelper {
133  typedef FlagType ReturnType;
134  template<typename DataType>
135  FlagType flags(DataAdapterType&, DataType const& data) {
136  return FlagType(data.template dimension<data::Time>(), data.template dimension<data::Frequency>());
137  }
138  };
139 
140  // unless a flags structure already exists
141  template<typename DataAdapterType>
142  struct FlagTypeHelper<DataAdapterType, typename std::enable_if<std::is_convertible<typename panda::is_pointer_wrapper<typename DataAdapterType::ReturnType>::type, FlagType&>::value>::type> {
143  typedef FlagType& ReturnType;
144  template<typename DataType>
145  FlagType& flags(DataAdapterType& adapter, DataType const&) {
146  return static_cast<FlagType&>(*adapter.data());
147  }
148  };
149 
150 } // namespace
151 
152 template<typename RfimTraits>
153 template<typename DataType>
154 void Rfim<RfimTraits>::operator()(DataType& data, DataAdapter& adapter)
155 {
156  float mean, std_dev, median;
157  calculateStatistics(data.begin(), data.end(), mean, median, std_dev);
158  float factor;
159  if (std_dev == 0.0f) {
160  factor = _config.base_sensitivity();
161  } else {
162  factor = std_dev * _config.base_sensitivity();
163  }
164 
165  unsigned extra_flagged = 0;
166 
167  FlagTypeHelper<DataAdapter> flag_helper;
168  typename FlagTypeHelper<DataAdapter>::ReturnType flags = flag_helper.flags(adapter, data);
169  for (unsigned window : _config.thresholding_data_sizes()) {
170  float thresholdI = median + calcThresholdI(_config.its_cutoff_threshold(), window, 1.5f) * factor;
171  extra_flagged += sumThreshold2DHorizontal(data, flags, window, thresholdI);
172  extra_flagged += sumThreshold2DVertical(data, flags, window, thresholdI);
173  }
174 
175 
176  // Ensure adpater interface is called //TODO this should really happen during rfim discovery rather than as a final step
177  auto it=data.cbegin();
178  auto flag_it =flags.cbegin();
179  for ( data::DimensionIndex<data::Time> spectrum_num(0); spectrum_num < data.number_of_spectra(); ++spectrum_num ) {
180  std::size_t bad_channel=0;
181  for ( data::DimensionIndex<data::Frequency> channel_num(0); channel_num < data.number_of_channels(); ++channel_num ) {
182  if(*flag_it) {
183  adapter.mark_bad(*it, channel_num);
184  ++bad_channel;
185  }
186  ++it;
187  ++flag_it;
188  }
189  // we have a good spectrum
190  if(bad_channel==0) adapter.mark_good(data.spectrum(spectrum_num));
191  }
192 }
193 
194 
195 } // namespace sum_threshold
196 } // namespace rfim
197 } // namespace cheetah
198 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35