Cheetah - SKA - PSS - Prototype Time Domain Search Pipeline
Rfim.cpp
1 /* The MIT License (MIT)
2  *
3  * Copyright (c) 2016 The SKA organisation
4  *
5  * Permission is hereby granted, free of charge, to any person obtaining a copy
6  * of this software and associated documentation files (the "Software"), to deal
7  * in the Software without restriction, including without limitation the rights
8  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
9  * copies of the Software, and to permit persons to whom the Software is
10  * furnished to do so, subject to the following conditions:
11  *
12  * The above copyright notice and this permission notice shall be included in all
13  * copies or substantial portions of the Software.
14  *
15  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
18  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
20  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
21  * SOFTWARE.
22  */
23 #include "cheetah/rfim/iqrmcpu/Rfim.h"
24 #include "panda/Log.h"
25 #include "panda/TypeTraits.h"
26 #include <cmath>
27 #include <algorithm>
28 
29 
30 
31 namespace ska {
32 namespace cheetah {
33 namespace rfim {
34 namespace iqrmcpu {
35 
36 namespace {
37 
38 template<typename T>
39 static inline double Lerp(T v0, T v1, T t)
40 {
41  return (1 - t)*v0 + t*v1;
42 }
43 
44 template<typename T>
45 static inline std::vector<T> quantile(const std::vector<T>& in_data, const std::vector<T>& probs)
46 {
47  if (in_data.empty())
48  {
49  return std::vector<T>();
50  }
51 
52  if (1 == in_data.size())
53  {
54  return std::vector<T>(1, in_data[0]);
55  }
56 
57  std::vector<T> data = in_data;
58  std::sort(data.begin(), data.end());
59  std::vector<T> quantiles;
60 
61  for (size_t i = 0; i < probs.size(); ++i)
62  {
63  T poi = Lerp<T>(-0.5, data.size() - 0.5, probs[i]);
64 
65  size_t left = std::max(int64_t(std::floor(poi)), int64_t(0));
66  size_t right = std::min(int64_t(std::ceil(poi)), int64_t(data.size() - 1));
67 
68  T datLeft = data.at(left);
69  T datRight = data.at(right);
70 
71  T quantile = Lerp<T>(datLeft, datRight, poi - left);
72 
73  quantiles.push_back(quantile);
74  }
75 
76  return quantiles;
77 }
78 
79 // Produces a suitable type for recording rfim flags
80 // default is to use a new structure
81 
82 template<typename DataAdapterType, typename TimeFrequencyType, typename Enable=void>
83 struct FlaggedDataTypeHelper {
84  typedef data::RfimFlaggedData< const TimeFrequencyType> FlaggedDataType;
85  typedef FlaggedDataType ReturnType;
86  template<typename DataType>
87  ReturnType flagged_data(DataAdapterType& , DataType const& data) {
88  return FlaggedDataType(data);
89  }
90 };
91 
92 // unless a flags structure already exists
93 template<typename DataAdapterType, typename TimeFrequencyType>
94 struct FlaggedDataTypeHelper<DataAdapterType, TimeFrequencyType, typename std::enable_if<std::is_convertible<typename panda::is_pointer_wrapper<typename DataAdapterType::ReturnType>::type, data::RfimFlaggedData<TimeFrequencyType>>::value>::type> {
95  typedef data::RfimFlaggedData<TimeFrequencyType> FlaggedDataType;
96  typedef FlaggedDataType& ReturnType;
97  template<typename DataType>
98  ReturnType flagged_data(DataAdapterType& adapter, DataType const&) {
99  return static_cast<FlaggedDataType&>(*adapter.data());
100  }
101 };
102 
103 } // namespace
104 
105 
106 
107 template<typename RfimTraits>
108 Rfim<RfimTraits>::Rfim(Config const& config)
109  : _config(config)
110 {
111 }
112 
113 template<typename RfimTraits>
114 Rfim<RfimTraits>::~Rfim()
115 {
116 }
117 
118 template<typename RfimTraits>
119 template<typename DataType>
120 inline void Rfim<RfimTraits>::operator()(DataType const& data, DataAdapter& adapter)
121 {
122  std::size_t data_size=data.number_of_channels() * data.number_of_spectra();
123  if(data_size == 0) return;
124 
125  // compute per channel std for the TF block
126  FlaggedDataTypeHelper<DataAdapter, DataType> flagged_data_helper;
127  typename FlaggedDataTypeHelper<DataAdapter, DataType>::ReturnType flagged_data = flagged_data_helper.flagged_data(adapter, data);
128 
129  // Compute IQR range
130  std::size_t flagged_data_channels = flagged_data.template dimension<data::Frequency>();
131  std::vector<bool> flags(flagged_data_channels, false);
132  std::vector<bool> flags_per_lag(flagged_data_channels, false);
133  std::vector<float> diff_vector(flagged_data_channels);
134 
135  for (std::size_t ii= 0; ii<= 2*_config.max_lag(); ++ii)
136  {
137 
138  //Rotate the vector
139  std::size_t jj = 0;
140 
141  auto const& channel_stats = flagged_data.channel_stats();
142  if (ii <= _config.max_lag())
143  {
144  std::size_t ind=0;
145  for (std::size_t kk=0; kk < channel_stats.size(); ++kk)
146  {
147  if (kk <= flagged_data_channels -( _config.max_lag() - ii) - 1)
148  diff_vector[kk] = std::sqrt(channel_stats[kk].variance) - std::sqrt(channel_stats[kk + (_config.max_lag() - ii)].variance);
149  else
150  {
151  diff_vector[kk] = std::sqrt(channel_stats[kk].variance) - std::sqrt(channel_stats[ind].variance);
152  ++ind;
153  }
154  }
155  }
156  else
157  {
158  std::size_t ind=0;
159  for (std::size_t kk=0; kk < channel_stats.size(); ++kk)
160  {
161  if (kk < ii)
162  {
163  diff_vector[kk] = std::sqrt(channel_stats[kk].variance) - std::sqrt(channel_stats[channel_stats.size() - 1 -ind].variance);
164  ++ind;
165  }
166  else
167  diff_vector[kk] = std::sqrt(channel_stats[kk].variance) - std::sqrt(channel_stats[kk - ii].variance);
168  }
169  }
170 
171  std::vector<float> quant = quantile<float>(diff_vector, {0.25, 0.5, 0.75});
172 
173  // Generate a difference vector for different lags
174  float std_dev = std::abs((quant[2] - quant[0]) / 1.349);
175 
176  std::generate(flags_per_lag.begin(), flags_per_lag.end(),
177  [&]()
178  {
179  bool flag = diff_vector[jj] - quant[1] > _config.nsigma()*std_dev;
180  ++jj;
181  return flag;
182  });
183 
184  //Logical OR on the an initial mask
185  for (std::size_t kk=0; kk < flags.size(); ++kk)
186  {
187  flags[kk] = flags_per_lag[kk] | flags[kk];
188  }
189  }
190 
191  // Write new flags via the adapter
192  float fraction = 0.0;
193 
194  for (std::size_t ii=0; ii < flags.size(); ++ii)
195  {
196  if (ii < _config.edge_channels() || ii > (flags.size()- _config.edge_channels()))
197  {
198  adapter.mark_bad(data.channel(ii));
199  continue;
200  }
201  if (flags[ii] == true)
202  {
203  adapter.mark_bad(data.channel(ii));
204  ++fraction;
205  }
206  }
207 
208  PANDA_LOG << "rfim::iqrm : Number of channels flagged: " << fraction/(float)flags.size();
209 
210 }
211 } // namespace iqrmcpu
212 } // namespace rfim
213 } // namespace cheetah
214 } // namespace ska
Some limits and constants for FLDO.
Definition: Brdz.h:35