diff --git a/antennaflagger/Flagger.cc b/antennaflagger/Flagger.cc index 0e4f20f1bef3bbcdd57b0190968fb3df37bad1b8..2874f75d95825151c42c4b969ec1fad0787082e1 100644 --- a/antennaflagger/Flagger.cc +++ b/antennaflagger/Flagger.cc @@ -17,10 +17,10 @@ #include "Flagger.h" namespace { -xt::xtensor<int, 1> FindOutliers(float sigma, size_t max_iterations, - xt::xtensor<float, 1>&& data) { +xt::xtensor<char, 1> FindOutliers(float sigma, size_t max_iterations, + xt::xtensor<float, 1>&& data) { // All non finite numbers are flagged as outlier - xt::xtensor<int, 1> flags = !xt::isfinite(data); + xt::xtensor<char, 1> flags = !xt::isfinite(data); // Take the median of the finite numbers only, otherwise the median is NaN. const float median = xt::median(xt::filter(data, !flags)); @@ -101,14 +101,14 @@ xt::xtensor<std::complex<float>, 3> Flagger::GroupStats( return stats_antenna; } -xt::xtensor<int, 2> Flagger::ComputeAntennaFlags( +xt::xtensor<char, 2> Flagger::ComputeAntennaFlags( float sigma, int max_iterations, const xt::xtensor<std::complex<float>, 3>& stats) { const size_t n_stations = stats.shape(0); const size_t n_antennas_per_station = stats.shape(1); const size_t n_correlations = stats.shape(2); - xt::xtensor<int, 2> flags({n_stations, n_antennas_per_station}, false); + xt::xtensor<char, 2> flags({n_stations, n_antennas_per_station}, false); for (size_t station = 0; station < n_stations; ++station) { for (size_t cor = 0; cor < n_correlations; ++cor) { @@ -124,7 +124,7 @@ xt::xtensor<int, 2> Flagger::ComputeAntennaFlags( return flags; } -xt::xtensor<int, 2> Flagger::ComputeStationFlags( +xt::xtensor<char, 2> Flagger::ComputeStationFlags( float sigma, int max_iterations, const xt::xtensor<std::complex<float>, 3>& stats) { const size_t n_stations = stats.shape(0); @@ -133,7 +133,7 @@ xt::xtensor<int, 2> Flagger::ComputeStationFlags( // In case of n_correlations == 4, use only XX and YY (index 0 and 3), // otherwise use only XX (index 0). const size_t n_correlations_out = n_correlations == 4 ? 2 : 1; - xt::xtensor<int, 2> flags({n_stations, n_correlations_out}, false); + xt::xtensor<char, 2> flags({n_stations, n_correlations_out}, false); const xt::xtensor<float, 3> stats_real = xt::real(stats); const xt::xtensor<float, 3> stats_imag = xt::imag(stats); @@ -197,20 +197,20 @@ void Flagger::AssertStatsComputed() const { assert(stats_antenna_stddev_.shape() == stats_antenna_sum_square_.shape()); } -xt::xtensor<int, 1> Flagger::FindBadAntennas(float sigma, int max_iterations) { +xt::xtensor<char, 1> Flagger::FindBadAntennas(float sigma, int max_iterations) { AssertStatsComputed(); find_bad_antennas_timer_.start(); - const xt::xtensor<int, 2> flags_stddev = Flagger::ComputeAntennaFlags( + const xt::xtensor<char, 2> flags_stddev = Flagger::ComputeAntennaFlags( sigma, max_iterations, stats_antenna_stddev_); assert(flags_stddev.shape(0) == n_stations_); assert(flags_stddev.shape(1) == n_antennas_per_station_); - const xt::xtensor<int, 2> flags_sum_square = Flagger::ComputeAntennaFlags( + const xt::xtensor<char, 2> flags_sum_square = Flagger::ComputeAntennaFlags( sigma, max_iterations, stats_antenna_sum_square_); - const xt::xtensor<int, 1> flagged_antennas = + const xt::xtensor<char, 1> flagged_antennas = xt::flatten(flags_stddev) & xt::flatten(flags_sum_square); find_bad_antennas_timer_.stop(); @@ -218,23 +218,23 @@ xt::xtensor<int, 1> Flagger::FindBadAntennas(float sigma, int max_iterations) { return flagged_antennas; } -xt::xtensor<int, 1> Flagger::FindBadStations(float sigma, int max_iterations) { +xt::xtensor<char, 1> Flagger::FindBadStations(float sigma, int max_iterations) { AssertStatsComputed(); find_bad_stations_timer_.start(); - const xt::xtensor<int, 2> flags_stddev = Flagger::ComputeStationFlags( + const xt::xtensor<char, 2> flags_stddev = Flagger::ComputeStationFlags( sigma, max_iterations, stats_antenna_stddev_); assert(flags_stddev.shape(0) == n_stations_); assert(flags_stddev.shape(1) == (n_correlations_ == 4 ? 2 : 1)); - const xt::xtensor<int, 2> flags_sum_square = Flagger::ComputeStationFlags( + const xt::xtensor<char, 2> flags_sum_square = Flagger::ComputeStationFlags( sigma, max_iterations, stats_antenna_sum_square_); const xt::xtensor<size_t, 1> station_indices = xt::flatten_indices( xt::where(xt::prod(flags_stddev & flags_sum_square, {1}))); - xt::xtensor<int, 1> antenna_flags({n_antennas_}, false); + xt::xtensor<char, 1> antenna_flags({n_antennas_}, false); for (size_t station : station_indices) { const size_t first_antenna = station * n_antennas_per_station_; diff --git a/antennaflagger/Flagger.h b/antennaflagger/Flagger.h index 3eada63d0589854c718ede86ee90244f68497de4..e3bde4d8ea4905473940819879d34500c7fe9fcc 100644 --- a/antennaflagger/Flagger.h +++ b/antennaflagger/Flagger.h @@ -63,7 +63,7 @@ class Flagger { * @return Flags as a tensor of booleans of length n_stations * * n_antennas_per_station. */ - xt::xtensor<int, 1> FindBadAntennas(float sigma, int max_iterations); + xt::xtensor<char, 1> FindBadAntennas(float sigma, int max_iterations); /** * Identify stations that are outliers compared to the other stations. @@ -74,7 +74,7 @@ class Flagger { * n_antennas_per_station. The booleans are 'int' instead of 'bool' since * XSimd otherwise converts booleans to integers and back. */ - xt::xtensor<int, 1> FindBadStations(float sigma, int max_iterations); + xt::xtensor<char, 1> FindBadStations(float sigma, int max_iterations); /** * Compute standard deviation for the real and imaginary component seperately @@ -116,7 +116,7 @@ class Flagger { * n_antennas_per_station). The flags are 'int' instead of 'bool' since * XSimd otherwise converts booleans to integers and back. */ - static xt::xtensor<int, 2> ComputeAntennaFlags( + static xt::xtensor<char, 2> ComputeAntennaFlags( float sigma, int max_iterations, const xt::xtensor<std::complex<float>, 3>& stats); @@ -130,7 +130,7 @@ class Flagger { * The booleans are 'int' instead of 'bool' since XSimd otherwise converts * booleans to integers and back. */ - static xt::xtensor<int, 2> ComputeStationFlags( + static xt::xtensor<char, 2> ComputeStationFlags( float sigma, int max_iterations, const xt::xtensor<std::complex<float>, 3>& stats); diff --git a/ddecal/gain_solvers/SolverTools.cc b/ddecal/gain_solvers/SolverTools.cc index 488129eabe2b9348ea80a2016ea6e917af50c579..79ba680b9948253a2b9610a7fbeadcd449eae9ba 100644 --- a/ddecal/gain_solvers/SolverTools.cc +++ b/ddecal/gain_solvers/SolverTools.cc @@ -45,7 +45,7 @@ void AssignAndWeight( // Use 'int' instead of 'bool' for flags since XSimd uses integers when // or'ing flags. Using integers directly prevents unpacking booleans to // integers beforehand and packing integers to booleans afterwards. - xt::xtensor<int, 2> flags( + xt::xtensor<char, 2> flags( {unweighted_data.shape(0), unweighted_data.shape(1)}, false); for (std::size_t correlation = 0; correlation < n_correlations; ++correlation) { diff --git a/steps/AntennaFlagger.cc b/steps/AntennaFlagger.cc index 4ecbf7ac19a12f99440e755fb30ae9ff3852bb59..55b46a1817edd090f009e0a7ff6e5dd6f787c432 100644 --- a/steps/AntennaFlagger.cc +++ b/steps/AntennaFlagger.cc @@ -100,7 +100,7 @@ bool AntennaFlagger::process(std::unique_ptr<base::DPBuffer> buffer) { flagger_->ComputeStats(buffer->GetData(), baseline_order_); // Find outlier antennas - const xt::xtensor<int, 1> antenna_flags = + const xt::xtensor<char, 1> antenna_flags = flagger_->FindBadAntennas(antenna_flagging_sigma_, antenna_flagging_max_iterations_) | flagger_->FindBadStations(station_flagging_sigma_, diff --git a/steps/PreFlagger.cc b/steps/PreFlagger.cc index 2f3574e9742cc38e7acb65d2c9dc85f2e9cbe8aa..db5ca2459de091b6d4d84dc2a768daaff964de71 100644 --- a/steps/PreFlagger.cc +++ b/steps/PreFlagger.cc @@ -131,7 +131,7 @@ bool PreFlagger::process(std::unique_ptr<DPBuffer> buffer) { assert(buffer->GetFlags().size() != 0); // Do the PSet steps and combine the result with the current flags. // Only count if the flag changes. - const xt::xtensor<int, 3>* const flags = + const xt::xtensor<char, 3>* const flags = itsPSet.process(*buffer, itsCount, xt::xtensor<bool, 1>(), itsTimer); switch (itsMode) { case Mode::kSetFlag: @@ -324,7 +324,7 @@ void PreFlagger::PSet::updateInfo(const DPInfo& info) { void PreFlagger::PSet::fillChannels(const DPInfo& info) { unsigned int nrcorr = info.ncorr(); unsigned int nrchan = info.nchan(); - xt::xtensor<int, 1> selChan(std::array<size_t, 1>{nrchan}); + xt::xtensor<char, 1> selChan(std::array<size_t, 1>{nrchan}); if (itsStrChan.empty()) { selChan.fill(true); } else { @@ -449,7 +449,7 @@ void PreFlagger::PSet::show(std::ostream& os, bool showName) const { } } -xt::xtensor<int, 3>* PreFlagger::PSet::process( +xt::xtensor<char, 3>* PreFlagger::PSet::process( DPBuffer& out, unsigned int timeSlot, const xt::xtensor<bool, 1>& matchBL, common::NSTimer& timer) { // No need to process it if the time mismatches or if only time selection. @@ -495,7 +495,7 @@ xt::xtensor<int, 3>* PreFlagger::PSet::process( return &itsFlags; } // Convert each baseline flag to a flag per correlation/channel. - int* flagPtr = itsFlags.data(); + char* flagPtr = itsFlags.data(); for (unsigned int i = 0; i < itsMatchBL.size(); ++i) { if (itsMatchBL[i]) { std::fill(flagPtr, flagPtr + nr, itsMatchBL[i]); @@ -525,17 +525,17 @@ xt::xtensor<int, 3>* PreFlagger::PSet::process( // are reused to AND or OR subexpressions. This can be done harmlessly // and saves the creation of too many arrays. if (!itsPSets.empty()) { - std::stack<xt::xtensor<int, 3>*> results; + std::stack<xt::xtensor<char, 3>*> results; for (int oper : itsRpn) { if (oper >= 0) { results.push(itsPSets[oper]->process(out, timeSlot, itsMatchBL, timer)); } else if (oper == OpNot) { - xt::xtensor<int, 3>* left = results.top(); + xt::xtensor<char, 3>* left = results.top(); *left = !*left; } else if (oper == OpOr || oper == OpAnd) { - xt::xtensor<int, 3>* right = results.top(); + xt::xtensor<char, 3>* right = results.top(); results.pop(); - xt::xtensor<int, 3>* left = results.top(); + xt::xtensor<char, 3>* left = results.top(); if (oper == OpOr) { *left |= *right; } else { @@ -550,7 +550,7 @@ xt::xtensor<int, 3>* PreFlagger::PSet::process( throw std::runtime_error( "Something went wrong while evaluating expression: results.size() != " "1"); - xt::xtensor<int, 3>* mflags = results.top(); + xt::xtensor<char, 3>* mflags = results.top(); itsFlags &= *mflags; } return &itsFlags; @@ -709,7 +709,7 @@ void PreFlagger::PSet::flagAmpl( const xt::xtensor<float, 3> amplitudes = xt::abs(data); const float* valPtr = amplitudes.data(); - int* flagPtr = itsFlags.data(); + char* flagPtr = itsFlags.data(); for (unsigned int i = 0; i < nr; ++i) { bool flag = false; for (unsigned int j = 0; j < nrcorr; ++j) { @@ -733,7 +733,7 @@ void PreFlagger::PSet::flagPhase( const xt::xtensor<float, 3> phases = xt::arg(data); const float* valPtr = phases.data(); - int* flagPtr = itsFlags.data(); + char* flagPtr = itsFlags.data(); for (unsigned int i = 0; i < nr; ++i) { bool flag = false; for (unsigned int j = 0; j < nrcorr; ++j) { @@ -756,7 +756,7 @@ void PreFlagger::PSet::flagReal( const std::size_t nrcorr = data.shape(2); const std::complex<float>* valPtr = data.data(); - int* flagPtr = itsFlags.data(); + char* flagPtr = itsFlags.data(); for (unsigned int i = 0; i < nr; ++i) { bool flag = false; for (unsigned int j = 0; j < nrcorr; ++j) { @@ -780,7 +780,7 @@ void PreFlagger::PSet::flagImag( const std::size_t nrcorr = data.shape(2); const std::complex<float>* valPtr = data.data(); - int* flagPtr = itsFlags.data(); + char* flagPtr = itsFlags.data(); for (unsigned int i = 0; i < nr; ++i) { bool flag = false; for (unsigned int j = 0; j < nrcorr; ++j) { @@ -1058,10 +1058,10 @@ void PreFlagger::PSet::fillBLMatrix() { } } -xt::xtensor<int, 1> PreFlagger::PSet::handleFreqRanges( +xt::xtensor<char, 1> PreFlagger::PSet::handleFreqRanges( const std::vector<double>& chanFreqs) { unsigned int nrchan = chanFreqs.size(); - xt::xtensor<int, 1> selChan({nrchan}, false); + xt::xtensor<char, 1> selChan({nrchan}, false); // A frequency range can be given as value..value or value+-value. // Units can be given for each value; if one is given it applies to both. // Default unit is MHz. diff --git a/steps/PreFlagger.h b/steps/PreFlagger.h index 7fc2d729e93d13d2f9f86793cddc3981871a66a1..afc5b6d188d8fb58457e9fb6f972177f829ad7b3 100644 --- a/steps/PreFlagger.h +++ b/steps/PreFlagger.h @@ -128,9 +128,9 @@ class PreFlagger : public Step { } /// Set and return the flags. - xt::xtensor<int, 3>* process(base::DPBuffer&, unsigned int timeSlot, - const xt::xtensor<bool, 1>& matchBL, - common::NSTimer& timer); + xt::xtensor<char, 3>* process(base::DPBuffer&, unsigned int timeSlot, + const xt::xtensor<bool, 1>& matchBL, + common::NSTimer& timer); /// Update the general info. /// It is used to adjust the parms if needed. @@ -208,7 +208,7 @@ class PreFlagger : public Step { /// Handle the frequency ranges given and determine which channels /// have to be flagged. - xt::xtensor<int, 1> handleFreqRanges(const std::vector<double>& chanFreqs); + xt::xtensor<char, 1> handleFreqRanges(const std::vector<double>& chanFreqs); /// Get the value and possible unit. /// If no unit is given, the argument is left untouched. @@ -269,7 +269,7 @@ class PreFlagger : public Step { xt::xtensor<bool, 2> itsChanFlags; ///< flags for channels to be flagged // Using 'int' instead of 'bool' in itsFlags works better with XSimd, since // XSimd uses an 'int' representation anyway for boolean operations. - xt::xtensor<int, 3> itsFlags; + xt::xtensor<char, 3> itsFlags; xt::xtensor<bool, 1> itsMatchBL; ///< true = baseline in buffer matches };