23 #include "cheetah/rfim/iqrmcpu/Rfim.h" 24 #include "panda/Log.h" 25 #include "panda/TypeTraits.h" 39 static inline double Lerp(T v0, T v1, T t)
41 return (1 - t)*v0 + t*v1;
45 static inline std::vector<T> quantile(
const std::vector<T>& in_data,
const std::vector<T>& probs)
49 return std::vector<T>();
52 if (1 == in_data.size())
54 return std::vector<T>(1, in_data[0]);
57 std::vector<T> data = in_data;
58 std::sort(data.begin(), data.end());
59 std::vector<T> quantiles;
61 for (
size_t i = 0; i < probs.size(); ++i)
63 T poi = Lerp<T>(-0.5, data.size() - 0.5, probs[i]);
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));
68 T datLeft = data.at(left);
69 T datRight = data.at(right);
71 T quantile = Lerp<T>(datLeft, datRight, poi - left);
73 quantiles.push_back(quantile);
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);
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());
107 template<
typename RfimTraits>
108 Rfim<RfimTraits>::Rfim(Config
const& config)
113 template<
typename RfimTraits>
114 Rfim<RfimTraits>::~Rfim()
118 template<
typename RfimTraits>
119 template<
typename DataType>
120 inline void Rfim<RfimTraits>::operator()(DataType
const& data, DataAdapter& adapter)
122 std::size_t data_size=data.number_of_channels() * data.number_of_spectra();
123 if(data_size == 0)
return;
126 FlaggedDataTypeHelper<DataAdapter, DataType> flagged_data_helper;
127 typename FlaggedDataTypeHelper<DataAdapter, DataType>::ReturnType flagged_data = flagged_data_helper.flagged_data(adapter, data);
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);
135 for (std::size_t ii= 0; ii<= 2*_config.max_lag(); ++ii)
141 auto const& channel_stats = flagged_data.channel_stats();
142 if (ii <= _config.max_lag())
145 for (std::size_t kk=0; kk < channel_stats.size(); ++kk)
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);
151 diff_vector[kk] = std::sqrt(channel_stats[kk].variance) - std::sqrt(channel_stats[ind].variance);
159 for (std::size_t kk=0; kk < channel_stats.size(); ++kk)
163 diff_vector[kk] = std::sqrt(channel_stats[kk].variance) - std::sqrt(channel_stats[channel_stats.size() - 1 -ind].variance);
167 diff_vector[kk] = std::sqrt(channel_stats[kk].variance) - std::sqrt(channel_stats[kk - ii].variance);
171 std::vector<float> quant = quantile<float>(diff_vector, {0.25, 0.5, 0.75});
174 float std_dev = std::abs((quant[2] - quant[0]) / 1.349);
176 std::generate(flags_per_lag.begin(), flags_per_lag.end(),
179 bool flag = diff_vector[jj] - quant[1] > _config.nsigma()*std_dev;
185 for (std::size_t kk=0; kk < flags.size(); ++kk)
187 flags[kk] = flags_per_lag[kk] | flags[kk];
192 float fraction = 0.0;
194 for (std::size_t ii=0; ii < flags.size(); ++ii)
196 if (ii < _config.edge_channels() || ii > (flags.size()- _config.edge_channels()))
198 adapter.mark_bad(data.channel(ii));
201 if (flags[ii] ==
true)
203 adapter.mark_bad(data.channel(ii));
208 PANDA_LOG <<
"rfim::iqrm : Number of channels flagged: " << fraction/(float)flags.size();
Some limits and constants for FLDO.