Line data Source code
1 : //# FlagAgentRFlag.cc: This file contains the implementation of the FlagAgentAntennaIntegrations class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2017, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2017, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 :
23 : #include <flagging/Flagging/FlagAgentAntennaIntegrations.h>
24 :
25 : #include <sstream>
26 :
27 : #include <casacore/casa/Quanta/MVTime.h>
28 : #include <casacore/casa/Utilities/DataType.h>
29 :
30 : namespace casa { //# NAMESPACE CASA - BEGIN
31 :
32 : /**
33 : * Will use the per-row iteration approach, with vis-buffer preprocessing.
34 : *
35 : * @param dh Data handler, just passed to the base class
36 : * @param config Flag agent configuration object, antint specific parameters are processed
37 : * @param writePrivateFlagCube Write option for the base class
38 : * @param flag Flag for the base class
39 : */
40 8 : FlagAgentAntennaIntegrations::FlagAgentAntennaIntegrations(
41 : FlagDataHandler *dh, casacore::Record config,
42 : casacore::Bool writePrivateFlagCube,
43 8 : casacore::Bool flag):
44 : FlagAgentBase(dh, config, FlagAgentBase::ROWS_PREPROCESS_BUFFER, writePrivateFlagCube,
45 8 : flag)
46 : {
47 8 : logger_p->origin(casacore::LogOrigin(agentName_p,__FUNCTION__,WHERE));
48 :
49 8 : setAgentParameters(config, flagDataHandler_p->antennaNames_p);
50 8 : }
51 :
52 : /**
53 : * Deal with parameters specific to 'antint'.
54 : *
55 : * minchanfrac: fraction of channels with flagged antenna to initiate flagging.
56 : * verbose: boolean to print the timestamps of the integrations that are flagged
57 : *
58 : * @param config A flagdata configuration/parameters object
59 : */
60 8 : void FlagAgentAntennaIntegrations::setAgentParameters(casacore::Record config,
61 : casacore::Vector<casacore::String> *antennaNames)
62 : {
63 8 : const auto fields = config.nfields();
64 8 : *logger_p << casacore::LogIO::NORMAL << "The configuration received by this agent has "
65 8 : << fields << " fields with the following values:" << casacore::LogIO::POST;
66 8 : std::ostringstream ostr;
67 8 : config.print(ostr);
68 8 : *logger_p << casacore::LogIO::NORMAL << ostr.str() << casacore::LogIO::POST;
69 :
70 :
71 8 : const auto minChanOpt = "minchanfrac";
72 8 : int found = config.fieldNumber(minChanOpt);
73 8 : if (found >= 0) {
74 8 : minChanThreshold_p = config.asDouble(minChanOpt);
75 : }
76 :
77 8 : const auto antOpt = "antint_ref_antenna";
78 8 : found = config.fieldNumber(antOpt);
79 8 : if (found >= 0) {
80 8 : casacore::String antintAnt = config.asString(antOpt);
81 8 : antIdx_p = findAntennaID(antintAnt, antennaNames);
82 8 : *logger_p << casacore::LogIO::NORMAL << "Found requested antenna of interest, with "
83 8 : "name: " << antintAnt << ", and index: " << antIdx_p << casacore::LogIO::POST;
84 :
85 8 : } else {
86 0 : throw casacore::AipsError(casacore::String("The parameter ") + antOpt + " was not "
87 0 : "given. This parameter is mandatory in this flagging mode.");
88 : }
89 :
90 8 : const auto verboseOpt = "verbose";
91 8 : found = config.fieldNumber(verboseOpt);
92 8 : if (found >= 0) {
93 7 : verbose_p = config.asBool(verboseOpt);
94 : }
95 8 : }
96 :
97 : /**
98 : * Find the antenna ID for an antenna name
99 : *
100 : * @param name the name of an antenna
101 : * @param antennaNames gives antena names by idx
102 : *
103 : * @return the antenna id if found for the given name
104 : *
105 : * @throws AipsError if no antenna is found by the name given..
106 : */
107 : casacore::uInt
108 8 : FlagAgentAntennaIntegrations::findAntennaID(const casacore::String &name,
109 : const casacore::Vector<casacore::String> *antennaNames) {
110 8 : casacore::uInt idx = 0;
111 8 : for( ; idx < antennaNames->nelements(); ++idx) {
112 8 : if (name == antennaNames->operator()(idx)) {
113 8 : break;
114 : }
115 : }
116 :
117 8 : if (idx >= antennaNames->nelements()) {
118 0 : throw casacore::AipsError("The antenna requested (" + name + ") could not be found in "
119 0 : "the input dataset");
120 : }
121 8 : return idx;
122 : }
123 :
124 : /*
125 : * Processes blocks of rows by every point in time,
126 : * Checking all the baselines to the antenna of interest
127 : *
128 : * @param visBuffer Buffer to pre-process
129 : */
130 : void
131 1285 : FlagAgentAntennaIntegrations::preProcessBuffer(const vi::VisBuffer2 &visBuffer)
132 : {
133 1285 : const auto &time_p = visBuffer.time();
134 1285 : auto rowCnt = time_p.size();
135 1285 : if (rowCnt <= 0)
136 0 : return;
137 :
138 1285 : auto antennasCnt = visBuffer.nAntennas();
139 1285 : auto channelsCnt = visBuffer.nChannels();
140 1285 : *logger_p << casacore::LogIO::DEBUG1 << "Pre-processing visibility buffer, with "
141 1285 : << " nRows: " << visBuffer.nRows()
142 1285 : << " nCorrelations: " << visBuffer.nCorrelations()
143 : << " nChannels: " << channelsCnt << ", nAntennas:: "
144 1285 : << antennasCnt << casacore::LogIO::POST;
145 :
146 : // TODO: use a casacore::Matrix here?
147 1285 : TableFlagPerBaselinePerChannel flagPerBaselinePerChannel;
148 1285 : flagPerBaselinePerChannel.resize(antennasCnt - 1);
149 5140 : for(auto &row : flagPerBaselinePerChannel) {
150 3855 : row.assign(channelsCnt, false);
151 : }
152 :
153 1285 : const auto ant1 = visBuffer.antenna1();
154 1285 : const auto ant2 = visBuffer.antenna2();
155 1285 : const auto &polChanRowFlags = visBuffer.flagCube();
156 1285 : auto currentTime = time_p[0];
157 1285 : auto baselineIdx = 0;
158 17443 : for (casacore::uInt rowIdx = 0; rowIdx < rowCnt; ++rowIdx) {
159 16158 : if (ant1[rowIdx] != antIdx_p and ant2[rowIdx] != antIdx_p)
160 8079 : continue;
161 :
162 8079 : auto rowTime = time_p[rowIdx];
163 8079 : if (rowTime != currentTime) {
164 1408 : doPreProcessingTimePoint(doFlagTime_p, currentTime, flagPerBaselinePerChannel);
165 1408 : currentTime = rowTime;
166 1408 : baselineIdx = 0;
167 : }
168 :
169 8079 : checkAnyPolarizationFlagged(polChanRowFlags, flagPerBaselinePerChannel,
170 8079 : rowIdx, baselineIdx++);
171 :
172 : }
173 : // last time
174 1285 : doPreProcessingTimePoint(doFlagTime_p, currentTime, flagPerBaselinePerChannel);
175 1285 : }
176 :
177 : /**
178 : * Checks channel by channel. Count if there are any polarizations flagged per channel
179 : * once all the rows for a time point are scanned, doPreProcessingTimePoint() will
180 : * decide whether to flag all the integrations of the antenna of interest
181 : *
182 : * @param polChanRowFlags flag cube from a row, with [polarizations, channels, rows]
183 : * @param flagPerBaselinePerChannel
184 : * @param row row index for
185 : * @param baselineIdx index for the table of baseline-channel flags. Counting from 0.
186 : */
187 8079 : void FlagAgentAntennaIntegrations::checkAnyPolarizationFlagged(const casacore::Cube<casacore::Bool> &polChanRowFlags,
188 : TableFlagPerBaselinePerChannel &flagPerBaselinePerChannel,
189 : casacore::uInt row, casacore::uInt baselineIdx) {
190 :
191 8079 : const auto channelCnt = polChanRowFlags.ncolumn();
192 8079 : const auto polCnt = polChanRowFlags.nrow();
193 283215 : for (casacore::uInt chanIdx = 0; chanIdx < channelCnt; ++chanIdx) {
194 1375680 : for (casacore::uInt polIdx = 0; polIdx < polCnt; ++polIdx) {
195 : // assume that all polarization products must be unflagged for a baseline to be
196 : // deemed unflagged
197 1100544 : if (polChanRowFlags(polIdx, chanIdx, row)) {
198 0 : flagPerBaselinePerChannel[baselineIdx][chanIdx] = true;
199 0 : break;
200 : }
201 : }
202 : }
203 8079 : }
204 :
205 : /**
206 : * Once all the rows for a point in time have been processed, we know
207 : * what baselines to the antenna of interest are flagged. If all these
208 : * baselines are flagged, note down that all integrations for this
209 : * point in time should be flagged.
210 : *
211 : * @param flaggedTimes Structure to record what times must be flagged
212 : * @param rowTime value of the TIME column
213 : * @param flagPerBaselinePerChannel table with flag values of every baseline for the
214 : * different channels
215 : */
216 : void
217 2693 : FlagAgentAntennaIntegrations::doPreProcessingTimePoint(FlaggedTimesMap &flaggedTimes,
218 : casacore::Double rowTime,
219 : TableFlagPerBaselinePerChannel &flagPerBaselinePerChannel)
220 : {
221 :
222 2693 : const auto channelCnt = flagPerBaselinePerChannel.front().size();
223 2693 : if (1 == channelCnt) {
224 0 : doPreProcessingTimePointSingleChannel(flaggedTimes, rowTime, flagPerBaselinePerChannel);
225 : } else {
226 2693 : doPreProcessingTimePointMultiChannel(flaggedTimes, rowTime, flagPerBaselinePerChannel);
227 : }
228 :
229 : // clear flag counting table
230 10772 : for (auto row : flagPerBaselinePerChannel) {
231 8079 : row.assign(row.size(), false);
232 8079 : }
233 2693 : }
234 :
235 : /**
236 : * Pre-process with multiple channels (requires calculating the
237 : * fraction of channels flagged).
238 : *
239 : * Same parameters as doPreProcessingTimePoint()
240 : */
241 : void
242 2693 : FlagAgentAntennaIntegrations::doPreProcessingTimePointMultiChannel(FlaggedTimesMap &flaggedTimes,
243 : casacore::Double rowTime,
244 : const TableFlagPerBaselinePerChannel &flagPerBaselinePerChannel)
245 : {
246 :
247 2693 : const auto baseCnt = flagPerBaselinePerChannel.size();
248 2693 : const auto channelCnt = flagPerBaselinePerChannel.front().size();
249 :
250 2693 : casacore::uInt flaggedChannelsCnt = 0;
251 94405 : for (casacore::uInt chanIdx = 0; chanIdx < channelCnt; ++chanIdx) {
252 91712 : auto allBaselinesFlagged = true;
253 91712 : for (casacore::uInt baseIdx = 0; baseIdx < baseCnt; ++baseIdx) {
254 91712 : if (!flagPerBaselinePerChannel[baseIdx][chanIdx]) {
255 91712 : allBaselinesFlagged = false;
256 91712 : break;
257 : }
258 : }
259 :
260 91712 : if (allBaselinesFlagged)
261 0 : flaggedChannelsCnt++;
262 : }
263 : // Baselines will be considered flagged if a fraction greater than a nominated fraction
264 : // of channels is flagged.
265 2693 : auto frac = static_cast<double>(flaggedChannelsCnt) / channelCnt;
266 2693 : if (frac <= minChanThreshold_p) {
267 2514 : return;
268 : }
269 :
270 : // all baselines are flagged for a fraction of channels greater than the threshold
271 179 : flaggedTimes[rowTime] = true;
272 179 : if (verbose_p) {
273 0 : casacore::MVTime time(rowTime/casacore::C::day);
274 0 : auto fracAlreadyFlagged = static_cast<float>(flaggedChannelsCnt) / channelCnt;
275 0 : *logger_p << casacore::LogIO::NORMAL << "Flagging integration at time: "
276 0 : << time.string(casacore::MVTime::YMD,7)
277 : << " (fraction of flagged channels found: " << fracAlreadyFlagged << ")"
278 0 : << casacore::LogIO::POST;
279 0 : }
280 : }
281 :
282 :
283 : /**
284 : * If the selected channel is flagged for all the baselines, mark this time for flagging.
285 : *
286 : * This method Is not exactly the same as the MultiChannel version. When minchanfrac=1
287 : * MultiChannel would never flag, as the comparison requires that the 'fraction of channels
288 : * flagged' be > minchanfrac.
289 : *
290 : * Same parameters as doPreProcessingTimePoint()
291 : */
292 : void
293 0 : FlagAgentAntennaIntegrations::doPreProcessingTimePointSingleChannel(FlaggedTimesMap &flaggedTimes,
294 : casacore::Double rowTime,
295 : const TableFlagPerBaselinePerChannel &flagPerBaselinePerChannel)
296 : {
297 :
298 0 : const auto baseCnt = flagPerBaselinePerChannel.size();
299 0 : for (casacore::uInt baseIdx = 0; baseIdx < baseCnt; ++baseIdx) {
300 : // A baseline is flagged if the only channel selected is flagged
301 : // Found a baseline that is not flagged:
302 0 : if (!flagPerBaselinePerChannel[baseIdx].front())
303 0 : return;
304 : }
305 :
306 : // all baselines are flagged for a fraction of channels greater than the threshold
307 0 : flaggedTimes[rowTime] = true;
308 0 : if (verbose_p) {
309 0 : casacore::MVTime time(rowTime/casacore::C::day);
310 0 : *logger_p << casacore::LogIO::NORMAL << "Flagging integration at time:"
311 0 : << time.string(casacore::MVTime::YMD,7) << casacore::LogIO::POST;
312 0 : }
313 : }
314 :
315 : bool
316 24702 : FlagAgentAntennaIntegrations::computeRowFlags(const vi::VisBuffer2 &visBuffer,
317 : FlagMapper &/*flags*/, casacore::uInt row)
318 : {
319 24702 : auto flag = false;
320 : // As per preprocessBuffer(), all rows for this point in time have to be flagged
321 24702 : const auto rowTime = visBuffer.time()[row];
322 24702 : if (doFlagTime_p.cend() != doFlagTime_p.find(rowTime)) {
323 1074 : flag = true;
324 : }
325 :
326 24702 : return flag;
327 : }
328 :
329 :
330 : } //# NAMESPACE CASA - END
|
