Line data Source code
1 : //# FlagAgentRFlag.cc: This file contains the implementation of the FlagAgentRFlag class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2012, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2012, 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/FlagAgentRFlag.h>
24 :
25 : using std::pair;
26 : using namespace casacore;
27 :
28 : namespace casa { //# NAMESPACE CASA - BEGIN
29 :
30 :
31 53 : FlagAgentRFlag::FlagAgentRFlag(FlagDataHandler *dh, Record config, Bool writePrivateFlagCube, Bool flag):
32 53 : FlagAgentBase(dh,config,ANTENNA_PAIRS,writePrivateFlagCube,flag)
33 : {
34 53 : setAgentParameters(config);
35 :
36 : // Request pre-loading spw
37 53 : flagDataHandler_p->preLoadColumn(VisBufferComponent2::SpectralWindows);
38 : // flagDataHandler_p->preLoadColumn(vi::Freq);
39 :
40 : // Initialize parameters for robust stats (spectral analysis)
41 53 : nIterationsRobust_p = 12;
42 53 : thresholdRobust_p = vector<Double>(nIterationsRobust_p);
43 53 : thresholdRobust_p[0] = 6.0;
44 53 : thresholdRobust_p[1] = 5.0;
45 53 : thresholdRobust_p[2] = 4.0;
46 53 : thresholdRobust_p[3] = 3.6;
47 53 : thresholdRobust_p[4] = 3.2;
48 53 : thresholdRobust_p[5] = 3.0;
49 53 : thresholdRobust_p[6] = 2.7;
50 53 : thresholdRobust_p[7] = 2.6;
51 53 : thresholdRobust_p[8] = 2.5;
52 53 : thresholdRobust_p[9] = 2.5;
53 53 : thresholdRobust_p[10] = 2.5;
54 53 : thresholdRobust_p[11] = 3.5;
55 53 : }
56 :
57 106 : FlagAgentRFlag::~FlagAgentRFlag()
58 : {
59 : // Compiler automagically calls FlagAgentBase::~FlagAgentBase()
60 106 : }
61 :
62 53 : void FlagAgentRFlag::setAgentParameters(Record config)
63 : {
64 53 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__));
65 :
66 : int exists;
67 :
68 : // Determine if we have to apply the flags in the modified flag cube
69 106 : String display("none");
70 53 : exists = config.fieldNumber ("display");
71 53 : if (exists >= 0)
72 : {
73 53 : if( config.type(exists) != TpString )
74 : {
75 0 : throw( AipsError ( "Parameter 'display' must be of type 'string'" ) );
76 : }
77 53 : display = config.asString("display");
78 53 : *logger_p << LogIO::NORMAL << " display is: " << display << LogIO::POST;
79 : }
80 :
81 : // Determine if we have to generate plot reports.
82 53 : if ( (display == String("report")) or (display == String("both")) )
83 : {
84 0 : doplot_p = true;
85 : }
86 : else
87 : {
88 53 : doplot_p = false;
89 : }
90 :
91 : // Do we actually flag, or just calculate thresholds without applying flags?
92 53 : Bool writeflags(display != "report");
93 53 : exists = config.fieldNumber ("writeflags");
94 53 : if (exists >= 0)
95 : {
96 53 : if( config.type(exists) != TpBool )
97 : {
98 0 : throw( AipsError ( "Parameter 'writeflags' must be of type 'bool'" ) );
99 : }
100 53 : writeflags = config.asBool("writeflags");
101 53 : *logger_p << LogIO::NORMAL << " writeflags is: " << writeflags << LogIO::POST;
102 : }
103 :
104 53 : if( (writeflags == true) or (display == String("data")) or (display == String("both")) )
105 : {
106 38 : doflag_p = true;
107 38 : *logger_p << LogIO::NORMAL << " (writeflags OR display)=(" << writeflags << "," << display << "), will apply flags on modified flag cube " << LogIO::POST;
108 : }
109 : else
110 : {
111 15 : doflag_p = false;
112 : }
113 :
114 : // AIPS RFlag FPARM(1)
115 53 : exists = config.fieldNumber ("winsize");
116 53 : if (exists >= 0)
117 : {
118 44 : if( config.type(exists) != TpInt )
119 : {
120 0 : throw( AipsError ( "Parameter 'winsize' must be of type 'Int'" ) );
121 : }
122 :
123 44 : nTimeSteps_p = config.asuInt("winsize");
124 : }
125 : else
126 : {
127 9 : nTimeSteps_p = 3;
128 : }
129 :
130 53 : *logger_p << logLevel_p << " winsize is " << nTimeSteps_p << LogIO::POST;
131 :
132 : // AIPS RFlag FPARM(5)
133 53 : exists = config.fieldNumber ("spectralmax");
134 53 : if (exists >= 0)
135 : {
136 44 : if( config.type(exists) != TpDouble && config.type(exists) != TpFloat && config.type(exists) != TpInt)
137 : {
138 0 : throw( AipsError ( "Parameter 'spectralmax' must be of type 'double'" ) );
139 : }
140 :
141 44 : spectralmax_p = config.asDouble("spectralmax");
142 : }
143 : else
144 : {
145 9 : spectralmax_p = 1E6;
146 : }
147 :
148 53 : *logger_p << logLevel_p << " spectralmax is " << spectralmax_p << LogIO::POST;
149 :
150 : // AIPS RFlag FPARM(6)
151 53 : exists = config.fieldNumber ("spectralmin");
152 53 : if (exists >= 0)
153 : {
154 44 : if( config.type(exists) != TpDouble && config.type(exists) != TpFloat && config.type(exists) != TpInt )
155 : {
156 0 : throw( AipsError ( "Parameter 'spectralmin' must be of type 'double'" ) );
157 : }
158 :
159 44 : spectralmin_p = config.asDouble("spectralmin");
160 : }
161 : else
162 : {
163 9 : spectralmin_p = 0;
164 : }
165 :
166 53 : *logger_p << logLevel_p << " spectralmin is " << spectralmin_p << LogIO::POST;
167 :
168 : // AIPS RFlag OPTYPE (type of operation for spectral analysis)
169 106 : String optype("MEDIAN");
170 53 : optype_p = MEDIAN;
171 53 : spectralAnalysis_p = &FlagAgentRFlag::simpleMedian;
172 53 : exists = config.fieldNumber ("optype");
173 53 : if (exists >= 0)
174 : {
175 0 : if( config.type(exists) != TpString )
176 : {
177 0 : throw( AipsError ( "Parameter 'optype' must be of type 'string'" ) );
178 : }
179 :
180 0 : optype = config.asString("optype");
181 0 : optype.upcase();
182 :
183 0 : if (optype == String("MEDIAN"))
184 : {
185 0 : optype_p = MEDIAN;
186 0 : spectralAnalysis_p = &FlagAgentRFlag::simpleMedian;
187 0 : *logger_p << LogIO::NORMAL << " optype for spectral analysis is simple median" << LogIO::POST;
188 : }
189 0 : else if (optype == String("RMEDIAN"))
190 : {
191 0 : optype_p = ROBUST_MEDIAN;
192 0 : *logger_p << LogIO::NORMAL << " optype for spectral analysis is robust median" << LogIO::POST;
193 : }
194 0 : else if (optype == String("MEAN"))
195 : {
196 0 : optype_p = MEAN;
197 0 : *logger_p << LogIO::NORMAL << " optype for spectral analysis is simple mean" << LogIO::POST;
198 : }
199 0 : else if (optype == String("RMEAN"))
200 : {
201 0 : optype_p = ROBUST_MEAN;
202 0 : spectralAnalysis_p = &FlagAgentRFlag::robustMean;
203 0 : *logger_p << LogIO::NORMAL << " optype for spectral analysis is robust mean" << LogIO::POST;
204 : }
205 : else
206 : {
207 0 : optype_p = MEDIAN;
208 0 : spectralAnalysis_p = &FlagAgentRFlag::simpleMedian;
209 0 : *logger_p << LogIO::NORMAL << " optype for spectral analysis is simple median" << LogIO::POST;
210 : }
211 : }
212 :
213 : // AIPS RFlag FPARM(9)
214 53 : exists = config.fieldNumber ("timedevscale");
215 53 : if (exists >= 0)
216 : {
217 46 : if( config.type(exists) != TpDouble && config.type(exists) != TpFloat && config.type(exists) != TpInt)
218 : {
219 0 : throw( AipsError ( "Parameter 'timedevscale' must be of type 'double'" ) );
220 : }
221 :
222 46 : noiseScale_p = config.asDouble("timedevscale");
223 : }
224 : else
225 : {
226 7 : noiseScale_p = 5;
227 : }
228 :
229 53 : *logger_p << logLevel_p << " timedevscale is " << noiseScale_p << LogIO::POST;
230 53 : const auto thresholds_note = "Note: not applying flags, simply calculating "
231 : "thresholds (action not 'apply')";
232 53 : if (not writeflags) {
233 15 : *logger_p << logLevel_p << " " << thresholds_note <<
234 15 : " so timedevscale is not used." << LogIO::POST;
235 : }
236 :
237 : // AIPS RFlag FPARM(10)
238 53 : exists = config.fieldNumber ("freqdevscale");
239 53 : if (exists >= 0)
240 : {
241 48 : if( config.type(exists) != TpDouble && config.type(exists) != TpFloat && config.type(exists) != TpInt )
242 : {
243 0 : throw( AipsError ( "Parameter 'freqdevscale' must be of type 'double'" ) );
244 : }
245 :
246 48 : scutoffScale_p = config.asDouble("freqdevscale");
247 : }
248 : else
249 : {
250 5 : scutoffScale_p = 5;
251 : }
252 :
253 53 : *logger_p << logLevel_p << " freqdevscale is " << scutoffScale_p << LogIO::POST;
254 53 : if (not writeflags) {
255 15 : *logger_p << logLevel_p << " " << thresholds_note <<
256 15 : "so freqdevscale is not used." << LogIO::POST;
257 : }
258 :
259 : // timedev - Matrix for time analysis deviation thresholds - (old AIPS RFlag FPARM(3)/NOISE)
260 53 : noise_p = 0;
261 53 : exists = config.fieldNumber ("timedev");
262 53 : bool nonempty = exists >= 0;
263 : // Special empty values
264 53 : if (nonempty) {
265 : // when no value is given to timedev, it 'exists' but as an empty string
266 47 : if (casacore::TpString == config.type(exists)) {
267 26 : auto value = config.asString(exists);
268 13 : if (0 == value.length()) {
269 12 : nonempty = false;
270 : }
271 : }
272 : // "timedev=[]" arrives here as an empty array of bool!
273 34 : else if (casacore::TpArrayBool == config.type(exists)) {
274 42 : auto tdev = config.asArrayBool(RecordFieldId("timedev"));
275 21 : if (0 == tdev.size()) {
276 21 : nonempty = false;
277 : }
278 : }
279 : }
280 53 : if (nonempty)
281 : {
282 14 : if ( config.type( exists ) == casacore::TpFloat || config.type( exists ) == casacore::TpDouble || config.type(exists) == casacore::TpInt )
283 : {
284 2 : noise_p = config.asDouble("timedev");
285 2 : *logger_p << logLevel_p << " timedev (same for all fields and spws) "
286 2 : "is " << noise_p << LogIO::POST;
287 : }
288 12 : else if( config.type(exists) == casacore::TpArrayDouble || config.type(exists) == casacore::TpArrayFloat || config.type(exists) == casacore::TpArrayInt)
289 : {
290 22 : Matrix<Double> timedev = config.asArrayDouble( RecordFieldId("timedev") );
291 11 : if(timedev.ncolumn()==3)
292 : {
293 11 : *logger_p << logLevel_p << " timedev [field,spw,dev] is "
294 11 : << timedev << LogIO::POST;
295 :
296 22 : IPosition shape = timedev.shape();
297 11 : uInt nDevs = shape[0];
298 33 : for(uInt dev_i=0; dev_i<nDevs; ++dev_i)
299 : {
300 0 : auto field_spw = std::make_pair((Int)timedev(dev_i, 0),
301 22 : (Int)timedev(dev_i, 1));
302 22 : field_spw_noise_map_p[field_spw] = timedev(dev_i, 2);
303 22 : user_field_spw_noise_map_p[field_spw] = true;
304 22 : *logger_p << LogIO::DEBUG1 << "timedev matrix - field=" << timedev(dev_i,0) << " spw=" << timedev(dev_i,1) << " dev=" << timedev(dev_i,2) << LogIO::POST;
305 : }
306 : }
307 : else
308 : {
309 0 : throw( AipsError( " Matrix for time analysis deviation thresholds (timedev) must have 3 columns [[field,spw,dev]]. Set to [] to use automatically-computed thresholds." ) );
310 : }
311 : }
312 : else
313 : {
314 1 : if (writeflags) {
315 : // action='calculate' (not writeflags) will pass the file name if
316 : // given. But we still want to accept the values to compute stats,
317 : // if they're valid
318 : // (see for example test_rflag_CAS_5037 / test_rflag_return_dict1)
319 : throw AipsError("The timedev value given cannot be interpreted as a "
320 : "numerical value or array of numerical values. Refusing "
321 0 : "to run RFlag!");
322 : }
323 : }
324 : }
325 : else
326 : {
327 39 : *logger_p << logLevel_p << "No timedev value given. Will Use automatically "
328 39 : "computed values if applying flags." << LogIO::POST;
329 : // noise_p initialized to 0 above.
330 : }
331 :
332 :
333 : // freqdev - Matrix for time analysis deviation thresholds (freqdev) - (old AIPS RFlag FPARM(4)/SCUTOFF)
334 53 : scutoff_p = 0;
335 53 : exists = config.fieldNumber ("freqdev");
336 53 : nonempty = exists >= 0;
337 : // Special empty values
338 53 : if (nonempty) {
339 : // when no value is given to freqdev, it 'exists' but as an empty string
340 47 : if (casacore::TpString == config.type(exists)) {
341 24 : auto value = config.asString(exists);
342 12 : if (0 == value.length()) {
343 11 : nonempty = false;
344 : }
345 : }
346 : // "freqdev=[]" arrives here as an empty array of bool!
347 35 : else if (casacore::TpArrayBool == config.type(exists)) {
348 48 : auto fdev = config.asArrayBool(RecordFieldId("freqdev"));
349 24 : if (0 == fdev.size()) {
350 24 : nonempty = false;
351 : }
352 : }
353 : }
354 53 : if (nonempty)
355 : {
356 12 : if ( config.type( exists ) == casacore::TpFloat || config.type( exists ) == casacore::TpDouble || config.type(exists) == casacore::TpInt )
357 : {
358 2 : scutoff_p = config.asDouble("freqdev");
359 2 : *logger_p << logLevel_p << " freqdev (same for all fields and spws) "
360 2 : "is " << scutoff_p << LogIO::POST;
361 : }
362 10 : else if( config.type(exists) == casacore::TpArrayDouble || config.type(exists) == casacore::TpArrayFloat || config.type(exists) == casacore::TpArrayInt)
363 : {
364 18 : Matrix<Double> freqdev = config.asArrayDouble( RecordFieldId("freqdev") );
365 9 : if(freqdev.ncolumn()==3)
366 : {
367 9 : *logger_p << logLevel_p << " freqdev [field,spw,dev] is " <<
368 9 : freqdev << LogIO::POST;
369 :
370 18 : IPosition shape = freqdev.shape();
371 9 : uInt nDevs = shape[0];
372 27 : for(uInt dev_i=0; dev_i<nDevs; ++dev_i)
373 : {
374 0 : auto field_spw = std::make_pair(static_cast<Int>(freqdev(dev_i, 0)),
375 18 : static_cast<Int>(freqdev(dev_i, 1)));
376 18 : field_spw_scutoff_map_p[field_spw] = freqdev(dev_i, 2);
377 18 : user_field_spw_scutoff_map_p[field_spw] = true;
378 18 : *logger_p << LogIO::DEBUG1 << "freqdev matrix - field=" << freqdev(dev_i,0) << " spw=" << freqdev(dev_i,1) << " dev=" << freqdev(dev_i,2) << LogIO::POST;
379 : }
380 : }
381 : else
382 : {
383 0 : throw( AipsError( " Matrix for spectral analysis deviation thresholds (freqdev) must have 3 columns [[field,spw,dev]]. Set to [] to use automatically-computed thresholds." ) );
384 : }
385 : }
386 : else
387 : {
388 1 : if (writeflags) {
389 : throw AipsError("The freqdev value given cannot be interpreted as a "
390 : "numerical value or array of numerical values. Refusing "
391 0 : "to run RFlag!");
392 : }
393 : }
394 : }
395 : else
396 : {
397 41 : *logger_p << logLevel_p << "No freqdev value given. Will Use automatically "
398 41 : "computed values if applying flags." << LogIO::POST;
399 : // scutoff initialized to zero above.
400 : }
401 :
402 106 : return;
403 : }
404 :
405 0 : Double FlagAgentRFlag::mean(vector<Double> &data,vector<Double> &counts)
406 : {
407 0 : Double sumAvg = 0;
408 0 : Double avg = 0;
409 :
410 0 : if (data.size() == 0)
411 : {
412 0 : for (size_t index = 0; index < data.size(); ++index)
413 : {
414 0 : if (counts[index] > 0)
415 : {
416 0 : sumAvg += data[index]/counts[index];
417 : }
418 : }
419 0 : avg = sumAvg/data.size();
420 : }
421 : else
422 : {
423 0 : avg = 0;
424 : }
425 :
426 0 : return avg;
427 : }
428 :
429 4602104 : void FlagAgentRFlag::noiseVsRef(vector<Double> &data, Double ref)
430 : {
431 287271728 : for (size_t index = 0; index < data.size(); ++index)
432 : {
433 282669624 : data[index] -= ref;
434 : }
435 :
436 4602104 : return;
437 : }
438 :
439 9204932 : Double FlagAgentRFlag::median(vector<Double> &data)
440 : {
441 : Double med;
442 9204932 : vector<Double> datacopy = data;
443 9204932 : sort(data.begin(),data.end());
444 :
445 9204932 : if (data.size() == 0)
446 : {
447 0 : med = 0;
448 : }
449 9204932 : else if (data.size() % 2 == 1)
450 : {
451 34524 : med = data[(data.size()-1)/2];
452 : }
453 : else
454 : {
455 9170408 : med = 0.5*(data[data.size()/2] + data[(data.size()/2)-1]);
456 : }
457 :
458 18409864 : return med;
459 : }
460 :
461 362 : Double FlagAgentRFlag::computeThreshold(vector<Double> &data,vector<Double> &/*dataSquared*/,vector<Double> &counts)
462 : {
463 : // Declare working variables
464 : Double avg;//,avgSquared,std;
465 :
466 : // Produce samples for median
467 724 : vector<Double> samplesForMedian(data.size(),0);
468 18340 : for (size_t index = 0; index < data.size(); ++index)
469 : {
470 17978 : if (counts[index] > 0)
471 : {
472 17724 : avg = data[index]/counts[index];
473 17724 : samplesForMedian[index] = avg;
474 : }
475 : }
476 :
477 : // Compute median
478 362 : Double med = median(samplesForMedian);
479 :
480 : // Produce samples for median absolute deviation
481 362 : vector<Double> samplesForMad(samplesForMedian.size(),0);
482 18340 : for (size_t index = 0; index < samplesForMedian.size(); ++index)
483 : {
484 17978 : samplesForMad[index] = abs(samplesForMedian[index] - med);
485 : }
486 :
487 : // Compute median absolute deviation
488 362 : Double mad = median(samplesForMad);
489 :
490 724 : return (med + 1.4826*mad);
491 : }
492 :
493 : /**
494 : * This will calculate the timedev and freqdev thresholds, using the
495 : * RMS/variance values calculated in computeAntennaPairFlagsCore(()
496 : * (into the field_spw_noise_histogram_).
497 : *
498 : * @return a FlagReport with the overall thresholds (timedev and
499 : * freqdev thresholds for every field-SPW pair).
500 : */
501 53 : FlagReport FlagAgentRFlag::getReport()
502 : {
503 106 : FlagReport totalRep(String("list"),agentName_p);
504 :
505 53 : if (doflag_p) {
506 38 : return totalRep;
507 : }
508 :
509 45 : FlagReport plotRepCont(String("list"),agentName_p);
510 :
511 : // Calculate thresholds
512 30 : generateThresholds( field_spw_noise_histogram_sum_p,
513 15 : field_spw_noise_histogram_sum_squares_p,
514 15 : field_spw_noise_histogram_counts_p,
515 15 : field_spw_noise_map_p,
516 : "Time analysis");
517 30 : generateThresholds( field_spw_scutoff_histogram_sum_p,
518 15 : field_spw_scutoff_histogram_sum_squares_p,
519 15 : field_spw_scutoff_histogram_counts_p,
520 15 : field_spw_scutoff_map_p,
521 : "Spectral analysis");
522 :
523 : // Threshold reports (should be returned if params were calculated)
524 30 : Record threshList;
525 15 : Int nEntriesNoise = field_spw_noise_map_p.size();
526 15 : Int nEntriesScutoff = field_spw_scutoff_map_p.size();
527 :
528 15 : if(nEntriesNoise > 0 || nEntriesScutoff > 0)
529 : {
530 30 : Matrix<Double> timedev(nEntriesNoise,3), freqdev(nEntriesScutoff,3);
531 15 : Int threshCountTime = 0, threshCountFreq = 0;
532 15 : pair<Int,Int> field_spw;
533 45 : for (auto spw_field_iter = field_spw_noise_map_p.begin();
534 45 : spw_field_iter != field_spw_noise_map_p.end();
535 30 : ++spw_field_iter)
536 : {
537 30 : field_spw = spw_field_iter->first;
538 :
539 30 : timedev(threshCountTime,0) = field_spw.first;
540 30 : timedev(threshCountTime,1) = field_spw.second;
541 30 : timedev(threshCountTime,2) = field_spw_noise_map_p[field_spw];
542 30 : ++threshCountTime;
543 : }
544 :
545 47 : for (auto spw_field_iter = field_spw_scutoff_map_p.begin();
546 47 : spw_field_iter != field_spw_scutoff_map_p.end();
547 32 : ++spw_field_iter)
548 : {
549 32 : field_spw = spw_field_iter->first;
550 :
551 32 : freqdev(threshCountFreq,0) = field_spw.first;
552 32 : freqdev(threshCountFreq,1) = field_spw.second;
553 32 : freqdev(threshCountFreq,2) = field_spw_scutoff_map_p[field_spw];
554 32 : ++threshCountFreq;
555 : }
556 :
557 15 : threshList.define( RecordFieldId("timedev") , timedev );
558 15 : threshList.define( RecordFieldId("freqdev") , freqdev );
559 :
560 30 : FlagReport returnThresh("rflag",agentName_p, threshList);
561 15 : totalRep.addReport(returnThresh);
562 : }
563 :
564 : // Add the plot-reports last. Display needs plot-reports to be at the end of the list
565 : // Should be fixed in the display agent....
566 15 : if(doplot_p==true)
567 : {
568 : // Plot reports (should be returned if params were calculated and display is activated)
569 0 : getReportCore(field_spw_noise_histogram_sum_p,
570 0 : field_spw_noise_histogram_sum_squares_p,
571 0 : field_spw_noise_histogram_counts_p,
572 0 : field_spw_noise_map_p,
573 : plotRepCont,
574 : "Time analysis",
575 : noiseScale_p);
576 0 : getReportCore(field_spw_scutoff_histogram_sum_p,
577 0 : field_spw_scutoff_histogram_sum_squares_p,
578 0 : field_spw_scutoff_histogram_counts_p,
579 0 : field_spw_scutoff_map_p,
580 : plotRepCont,
581 : "Spectral analysis",
582 : scutoffScale_p);
583 :
584 0 : Int nReports = plotRepCont.nReport();
585 0 : for (Int report_idx=0; report_idx<nReports; ++report_idx)
586 : {
587 0 : FlagReport report_i;
588 0 : Bool valid = plotRepCont.accessReport(report_idx,report_i);
589 0 : if (valid) totalRep.addReport(report_i);
590 : }
591 : }
592 :
593 15 : return totalRep;
594 : }
595 :
596 0 : FlagReport FlagAgentRFlag::getReportCore( map< pair<Int,Int>,vector<Double> > &data,
597 : map< pair<Int,Int>,vector<Double> > &dataSquared,
598 : map< pair<Int,Int>,vector<Double> > &counts,
599 : map< pair<Int,Int>,Double > &threshold,
600 : FlagReport &totalReport,
601 : string label,
602 : Double scale)
603 : {
604 : // Set logger origin
605 0 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__,WHERE));
606 :
607 : // First of all determine each SPW frequency in order to produce ordered vectors
608 0 : pair<Int,Int> current_field_spw;
609 0 : map< Int, vector<pair<Double,Int> > > field_spw_order;
610 0 : for (auto field_spw_iter = data.begin();
611 0 : field_spw_iter != data.end();
612 0 : ++field_spw_iter)
613 : {
614 0 : current_field_spw = field_spw_iter->first;
615 0 : auto freq_spw = std::make_pair(field_spw_frequencies_p[current_field_spw],current_field_spw.second);
616 0 : field_spw_order[current_field_spw.first].push_back(freq_spw);
617 : }
618 :
619 :
620 : // Now iterate over the field-frequency-spw map, and sort it on the fly
621 : Int field,spw;
622 : Double spwStd;
623 0 : String fieldName;
624 0 : for (auto field_freq_spw_iter = field_spw_order.begin();
625 0 : field_freq_spw_iter != field_spw_order.end();
626 0 : ++field_freq_spw_iter)
627 : {
628 : // Get field
629 0 : field = field_freq_spw_iter->first;
630 0 : fieldName = flagDataHandler_p->fieldNames_p->operator()(field);
631 :
632 : // Create specific report per field
633 0 : FlagReport fieldReport = FlagReport("plotpoints",agentName_p,fieldName + String(" - ") + String(label), "Frequency (GHz)", "Deviation");
634 :
635 : // Sort SPWs by ascending frequency
636 0 : sort(field_freq_spw_iter->second.begin(),field_freq_spw_iter->second.end());
637 :
638 : // Initialize tmp vectors
639 0 : vector<Double> total_frequency;
640 0 : vector<Double> total_threshold;
641 0 : vector<Double> total_threshold_squared;
642 0 : vector<Double> total_threshold_counts;
643 0 : vector<Float> total_threshold_spw_average;
644 0 : vector<Float> spw_separator_frequency;
645 0 : vector<Float> spw_separator_central;
646 0 : vector<Float> spw_separator_bar;
647 :
648 : // Log info
649 0 : *logger_p << logLevel_p << label.c_str() << " gathering stats for field "
650 0 : << field << " (" << fieldName << ")" << LogIO::POST;
651 :
652 : // Now iterate over SPWs
653 0 : for (uInt spw_i=0; spw_i<field_freq_spw_iter->second.size(); ++spw_i)
654 : {
655 : // Get SPW
656 0 : spw = field_freq_spw_iter->second[spw_i].second;
657 :
658 : // Create field-spw pair for accessing the data
659 0 : current_field_spw = std::make_pair(field,spw);
660 :
661 : // Access the data for this specific field-spw
662 0 : vector<Double> ¤t_spw_frequency = field_spw_frequency_p[current_field_spw];
663 0 : vector<Double> ¤t_spw_threshold = data[current_field_spw];
664 0 : vector<Double> ¤t_spw_threshold_squared = dataSquared[current_field_spw];
665 0 : vector<Double> ¤t_spw_threshold_counts = counts[current_field_spw];
666 :
667 : // Insert this field-spw data in total arrays
668 0 : total_frequency.insert(total_frequency.end(),current_spw_frequency.begin(),current_spw_frequency.end());
669 0 : total_threshold.insert(total_threshold.end(),current_spw_threshold.begin(),current_spw_threshold.end());
670 0 : total_threshold_squared.insert(total_threshold_squared.end(),current_spw_threshold_squared.begin(),current_spw_threshold_squared.end());
671 0 : total_threshold_counts.insert(total_threshold_counts.end(),current_spw_threshold_counts.begin(),current_spw_threshold_counts.end());
672 :
673 : // Prepare threshold array for plotting (note: display using scale factor)
674 0 : spwStd = scale * threshold[current_field_spw];
675 0 : vector<Float> aux(current_spw_threshold.size(),spwStd);
676 0 : total_threshold_spw_average.insert(total_threshold_spw_average.end(),aux.begin(),aux.end());
677 :
678 : // Prepare bars for separating SPWs
679 0 : spw_separator_frequency.push_back(current_spw_frequency[0]);
680 0 : spw_separator_central.push_back(0);
681 0 : spw_separator_bar.push_back(spwStd);
682 0 : spw_separator_frequency.push_back(current_spw_frequency[current_spw_frequency.size()-1]);
683 0 : spw_separator_central.push_back(0);
684 0 : spw_separator_bar.push_back(spwStd);
685 : }
686 :
687 : // Now copy values from std::vector to casa::vector
688 0 : size_t idx = 0;
689 0 : Double avg,sumSquare,variance = 0;
690 : // Vector<Float> threshold_index(total_threshold_counts.size(),0);
691 0 : Vector<Float> threshold_frequency(total_threshold_counts.size(),0);
692 0 : Vector<Float> threshold_avg(total_threshold_counts.size(),0);
693 0 : Vector<Float> threshold_up(total_threshold_counts.size(),0);
694 : //Vector<Float> threshold_down(total_threshold_counts.size(),0);
695 : //Vector<Float> threshold_variance(total_threshold_counts.size(),0);
696 0 : Vector<Float> threshold_dev(total_threshold_counts.size(),0);
697 0 : for (auto iter = total_threshold.begin(); iter != total_threshold.end(); ++iter)
698 : {
699 : // threshold_index(idx) = idx;
700 0 : threshold_frequency(idx) = total_frequency[idx];
701 0 : threshold_dev(idx) = total_threshold_spw_average[idx];
702 0 : if (total_threshold_counts[idx] > 0)
703 : {
704 0 : avg = total_threshold[idx]/total_threshold_counts[idx];
705 0 : sumSquare = total_threshold_squared[idx]/total_threshold_counts[idx];
706 : }
707 : else
708 : {
709 0 : avg = 0;
710 0 : sumSquare = 0;
711 : }
712 :
713 0 : variance = sumSquare - avg*avg;
714 0 : variance = sqrt(variance > 0? variance:0);
715 :
716 0 : threshold_avg(idx) = avg;
717 0 : threshold_up(idx) = avg+variance;
718 : //threshold_down(idx) = avg-variance;
719 : //threshold_variance(idx) = variance; // New
720 0 : ++idx;
721 : }
722 :
723 : // Finally add plots to field report
724 0 : fieldReport.addData("line", threshold_frequency,threshold_dev,"",Vector<Float>(),"rflag threshold");
725 0 : fieldReport.addData("line",threshold_frequency,threshold_up,"",Vector<Float>(),"median_deviation + variance");
726 0 : fieldReport.addData("scatter",threshold_frequency,threshold_avg,"",Vector<Float>(),"median_deviation");
727 : // fieldReport.addData("line", threshold_frequency,threshold_down,"",Vector<Float>(),"median_deviation - variance");
728 0 : Vector<Float> const xData(spw_separator_frequency);
729 0 : Vector<Float> const yData(spw_separator_central);
730 0 : Vector<Float> const error(spw_separator_bar);
731 0 : fieldReport.addData("scatter", xData, yData, "separator", error, "SPW separator");
732 :
733 : // Add bars to separate SPWs
734 0 : totalReport.addReport(fieldReport);
735 : }
736 :
737 0 : return totalReport;
738 : }
739 :
740 30 : void FlagAgentRFlag::generateThresholds(map< pair<Int,Int>,vector<Double> > &data,
741 : map< pair<Int,Int>,vector<Double> > &dataSquared,
742 : map< pair<Int,Int>,vector<Double> > &counts,
743 : map< pair<Int,Int>,Double > &threshold,
744 : string label)
745 : {
746 : // Set logger origin
747 30 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__,WHERE));
748 :
749 : // First of all determine each SPW frequency in order to produce ordered vectors
750 30 : pair<Int,Int> current_field_spw;
751 92 : for (auto field_spw_iter = data.begin(); field_spw_iter != data.end();
752 62 : ++field_spw_iter)
753 : {
754 62 : current_field_spw = field_spw_iter->first;
755 :
756 : // Compute threshold
757 62 : threshold[current_field_spw] = computeThreshold(data[current_field_spw],
758 62 : dataSquared[current_field_spw],
759 62 : counts[current_field_spw]);
760 :
761 : // Log info
762 62 : auto nfreq = field_spw_frequency_p[current_field_spw].size();
763 62 : auto freqStart = field_spw_frequency_p[current_field_spw][0];
764 62 : auto freqEnd = field_spw_frequency_p[current_field_spw][nfreq-1];
765 62 : *logger_p << logLevel_p << label.c_str() << " - Field " << current_field_spw.first
766 : << " - Spw " << current_field_spw.second
767 : << " - Frequency " << freqStart << "~" << freqEnd << "GHz"
768 : << " threshold (over baselines/timesteps) avg: "
769 62 : << threshold[current_field_spw] << LogIO::POST;
770 : }
771 :
772 60 : return;
773 : }
774 :
775 597208 : void FlagAgentRFlag::computeAntennaPairFlagsCore(pair<Int,Int> spw_field,
776 : Double noise,
777 : Double scutoff,
778 : uInt timeStart,
779 : uInt timeStop,
780 : uInt centralTime,
781 : VisMapper &visibilities,
782 : FlagMapper &flags)
783 : {
784 : // Get flag cube size
785 : Int nPols,nChannels,nTimesteps;
786 597208 : visibilities.shape(nPols, nChannels, nTimesteps);
787 :
788 : // Declare variables
789 597208 : Complex visibility;
790 597208 : Double SumWeight = 0;
791 597208 : Double SumWeightReal = 0;
792 597208 : Double SumWeightImag = 0;
793 597208 : Double StdTotal = 0;
794 597208 : Double SumReal = 0;
795 597208 : Double SumRealSquare = 0;
796 597208 : Double AverageReal = 0;
797 597208 : Double StdReal = 0;
798 597208 : Double SumImag = 0;
799 597208 : Double SumImagSquare = 0;
800 597208 : Double AverageImag = 0;
801 597208 : Double StdImag = 0;
802 597208 : Double deviationReal = 0;
803 597208 : Double deviationImag = 0;
804 :
805 : // Time-Direction analysis: Fix channel/polarization and compute stats with all time-steps
806 : // NOTE: It is better to operate in channel/polarization sequence for data contiguity
807 597208 : if ( (noise <= 0) or ((noise >= 0) and (doflag_p == true) and (prepass_p == false)) )
808 : {
809 37192752 : for (uInt chan_j=0; chan_j<(uInt)nChannels; ++chan_j)
810 : {
811 : // Compute variance
812 178313548 : for (uInt pol_k=0;pol_k<(uInt)nPols; ++pol_k)
813 : {
814 141715856 : SumWeight = 0;
815 141715856 : StdTotal = 0;
816 141715856 : SumReal = 0;
817 141715856 : SumRealSquare = 0;
818 141715856 : AverageReal = 0;
819 141715856 : StdReal = 0;
820 141715856 : SumImag = 0;
821 141715856 : SumImagSquare = 0;
822 141715856 : AverageImag = 0;
823 141715856 : StdImag = 0;
824 :
825 559332816 : for (uInt timestep_i=timeStart; timestep_i<=(uInt) timeStop; ++timestep_i)
826 : {
827 : // Ignore data point if it is already flagged
828 : // NOTE: In our case visibilities come w/o weights, so we check vs flags instead
829 417616960 : if (flags.getOriginalFlags(pol_k,chan_j,timestep_i)) continue;
830 :
831 417235916 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
832 417235916 : SumWeight += 1;
833 417235916 : SumReal += visibility.real();
834 417235916 : SumRealSquare += visibility.real()*visibility.real();
835 417235916 : SumImag += visibility.imag();
836 417235916 : SumImagSquare += visibility.imag()*visibility.imag();
837 : }
838 :
839 : // Now flag all timesteps if variance is greater than threshold
840 : // NOTE: In our case, SumWeight is zero when all the data is already flagged so we don't need to re-flag it
841 141715856 : if (SumWeight > 0)
842 : {
843 139082380 : AverageReal = SumReal / SumWeight;
844 139082380 : SumRealSquare = SumRealSquare / SumWeight;
845 139082380 : StdReal = SumRealSquare - AverageReal * AverageReal;
846 :
847 139082380 : AverageImag = SumImag / SumWeight;
848 139082380 : SumImagSquare = SumImagSquare / SumWeight;
849 139082380 : StdImag = SumImagSquare - AverageImag * AverageImag;
850 :
851 : // NOTE: It it not necessary to extract the square root if then we are going to pow2
852 : // StdReal = sqrt (StdReal > 0? StdReal:0);
853 : // StdImag = sqrt (StdImag > 0? StdImag:0);
854 : // StdTotal = sqrt (StdReal*StdReal + StdImag*StdImag);
855 139082380 : StdReal = StdReal > 0? StdReal:0;
856 139082380 : StdImag = StdImag > 0? StdImag:0;
857 139082380 : StdTotal = sqrt(StdReal + StdImag);
858 :
859 : // Apply flags or generate histogram?
860 : // NOTE: AIPS RFlag has the previous code duplicated in two separated
861 : // routines, but I don't see a reason to do this, performance-wise
862 139082380 : if (noise <= 0)
863 : {
864 69137708 : field_spw_noise_histogram_counts_p[spw_field][chan_j] += 1;
865 69137708 : field_spw_noise_histogram_sum_p[spw_field][chan_j] += StdTotal;
866 69137708 : field_spw_noise_histogram_sum_squares_p[spw_field][chan_j] += StdTotal*StdTotal;
867 : }
868 :
869 139082380 : if (noise >=0 and StdTotal > noise)
870 : {
871 9843518 : for (uInt timestep_i=timeStart; timestep_i<=timeStop; ++timestep_i)
872 : {
873 7382609 : if (!flags.getModifiedFlags(pol_k,chan_j,timestep_i))
874 : {
875 3077522 : flags.setModifiedFlags(pol_k,chan_j,timestep_i);
876 3077522 : visBufferFlags_p += 1;
877 : }
878 : }
879 : }
880 : }
881 : }
882 : }
883 : }
884 :
885 : // Spectral analysis: Fix timestep/polarization and compute stats with all channels
886 597208 : if ( (scutoff <= 0) or ((scutoff >= 0) and (doflag_p == true) and (prepass_p == false)) )
887 : {
888 1190120 : for (uInt timestep_i=centralTime; timestep_i<=centralTime; ++timestep_i)
889 : {
890 2902160 : for (uInt pol_k=0; pol_k<(uInt) nPols; ++pol_k)
891 : {
892 :
893 2307100 : (*this.*spectralAnalysis_p)( timestep_i,
894 : pol_k,
895 : nChannels,
896 : AverageReal,
897 : AverageImag,
898 : StdReal,
899 : StdImag,
900 : SumWeightReal,
901 : SumWeightImag,
902 : visibilities,
903 2307100 : flags);
904 :
905 2307100 : if (scutoff <= 0)
906 : {
907 71817222 : for (uInt chan_j=0; chan_j<(uInt) nChannels; ++chan_j)
908 : {
909 : // Ignore data point if it is already flagged
910 : // NOTE: In our case visibilities come w/o weights, so we check vs flags instead
911 70652456 : if (flags.getOriginalFlags(pol_k,chan_j,timestep_i)) continue;
912 :
913 70271412 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
914 :
915 70271412 : if (SumWeightReal > 0)
916 : {
917 70271412 : deviationReal = abs(visibility.real()-AverageReal);
918 70271412 : field_spw_scutoff_histogram_counts_p[spw_field][chan_j] += 1;
919 70271412 : field_spw_scutoff_histogram_sum_p[spw_field][chan_j] += deviationReal;
920 70271412 : field_spw_scutoff_histogram_sum_squares_p[spw_field][chan_j] += deviationReal*deviationReal;
921 : }
922 :
923 70271412 : if (SumWeightImag > 0)
924 : {
925 70271412 : deviationImag = abs(visibility.imag()-AverageImag);
926 70271412 : field_spw_scutoff_histogram_counts_p[spw_field][chan_j] += 1;
927 70271412 : field_spw_scutoff_histogram_sum_p[spw_field][chan_j] += deviationImag;
928 70271412 : field_spw_scutoff_histogram_sum_squares_p[spw_field][chan_j] += deviationImag*deviationImag;
929 : }
930 : }
931 : }
932 :
933 2307100 : if (scutoff >=0)
934 : {
935 : // Flag all channels?
936 1142334 : if ( (StdReal > spectralmax_p) or
937 1142334 : (StdImag > spectralmax_p) or
938 1142334 : (StdReal < spectralmin_p) or
939 1142334 : (StdImag < spectralmin_p) )
940 : {
941 0 : for (uInt chan_j=0; chan_j<(uInt) nChannels; ++chan_j)
942 : {
943 0 : if (!flags.getModifiedFlags(pol_k,chan_j,timestep_i))
944 : {
945 0 : flags.setModifiedFlags(pol_k,chan_j,timestep_i);
946 0 : visBufferFlags_p += 1;
947 : }
948 0 : }
949 : }
950 : // Check each channel separately vs the scutoff level
951 : else
952 : {
953 72205734 : for (uInt chan_j=0; chan_j<(uInt) nChannels; ++chan_j)
954 : {
955 71063400 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
956 140095465 : if ( (abs(visibility.real()-AverageReal)>scutoff) or
957 69032065 : (abs(visibility.imag()-AverageImag)>scutoff) )
958 : {
959 2884459 : if (!flags.getModifiedFlags(pol_k,chan_j,timestep_i))
960 : {
961 969415 : flags.setModifiedFlags(pol_k,chan_j,timestep_i);
962 969415 : visBufferFlags_p += 1;
963 : }
964 : }
965 : }
966 : }
967 : }
968 : }
969 : }
970 : }
971 :
972 1194416 : return;
973 : }
974 :
975 0 : void FlagAgentRFlag::robustMean( uInt timestep_i,
976 : uInt pol_k,
977 : uInt nChannels,
978 : Double &AverageReal,
979 : Double &AverageImag,
980 : Double &StdReal,
981 : Double &StdImag,
982 : Double &SumWeightReal,
983 : Double &SumWeightImag,
984 : VisMapper &visibilities,
985 : FlagMapper &flags)
986 : {
987 : // NOTE: To apply the robust coefficients we need some initial values of avg/std
988 : // In AIPS they simply use Std=1000 for the first iteration
989 : // I'm not very convinced with this but if we have to cross-validate...
990 0 : AverageReal = 0;
991 0 : AverageImag = 0;
992 0 : StdReal = 1000.0;
993 0 : StdImag = 1000.0;
994 0 : SumWeightReal = 0;
995 0 : SumWeightImag = 0;
996 :
997 : // Temporal working variables
998 0 : Complex visibility;
999 0 : Double SumReal = 0;
1000 0 : Double SumRealSquare = 0;
1001 0 : Double SumImag = 0;
1002 0 : Double SumImagSquare = 0;
1003 :
1004 0 : for (uInt robustIter=0; robustIter<nIterationsRobust_p; ++robustIter)
1005 : {
1006 0 : SumWeightReal = 0;
1007 0 : SumWeightImag = 0;
1008 0 : SumReal = 0;
1009 0 : SumRealSquare = 0;
1010 0 : SumImag = 0;
1011 0 : SumImagSquare = 0;
1012 :
1013 0 : for (uInt chan_j=0; chan_j<(uInt) nChannels; ++chan_j)
1014 : {
1015 : // Ignore data point if it is already flagged or weight is <= 0
1016 : // NOTE: In our case visibilities come w/o weights, so we check only vs flags
1017 0 : if (flags.getOriginalFlags(pol_k,chan_j,timestep_i)) continue;
1018 :
1019 0 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
1020 0 : if (abs(visibility.real()-AverageReal)<thresholdRobust_p[robustIter]*StdReal)
1021 : {
1022 0 : SumWeightReal += 1;
1023 0 : SumReal += visibility.real();
1024 0 : SumRealSquare += visibility.real()*visibility.real();
1025 : }
1026 :
1027 0 : if (abs(visibility.imag()-AverageImag)<thresholdRobust_p[robustIter]*StdImag)
1028 : {
1029 0 : SumWeightImag += 1;
1030 0 : SumImag += visibility.imag();
1031 0 : SumImagSquare += visibility.imag()*visibility.imag();
1032 : }
1033 : }
1034 :
1035 0 : if (SumWeightReal > 0)
1036 : {
1037 0 : AverageReal = SumReal / SumWeightReal;
1038 0 : SumRealSquare = SumRealSquare / SumWeightReal;
1039 0 : StdReal = SumRealSquare - AverageReal * AverageReal;
1040 0 : StdReal = sqrt(StdReal > 0? StdReal:0);
1041 : }
1042 : else
1043 : {
1044 0 : AverageReal = 0;
1045 0 : StdReal = 1000.0;
1046 : }
1047 :
1048 0 : if (SumWeightImag > 0)
1049 : {
1050 0 : AverageImag = SumImag / SumWeightImag;
1051 0 : SumImagSquare = SumImagSquare / SumWeightImag;
1052 0 : StdImag = SumImagSquare - AverageImag * AverageImag;
1053 0 : StdImag = sqrt(StdImag > 0? StdImag:0);
1054 : }
1055 : else
1056 : {
1057 0 : AverageImag = 0;
1058 0 : StdImag = 1000.0;
1059 : }
1060 : }
1061 :
1062 0 : return;
1063 : }
1064 :
1065 2307100 : void FlagAgentRFlag::simpleMedian( uInt timestep_i,
1066 : uInt pol_k,
1067 : uInt nChannels,
1068 : Double &AverageReal,
1069 : Double &AverageImag,
1070 : Double &StdReal,
1071 : Double &StdImag,
1072 : Double &SumWeightReal,
1073 : Double &SumWeightImag,
1074 : VisMapper &visibilities,
1075 : FlagMapper &flags)
1076 : {
1077 : // NOTE: To apply the robust coefficients we need some initial values of avg/std
1078 : // In AIPS they simply use Std=1000 for the first iteration
1079 : // I'm not very convinced with this but if we have to cross-validate...
1080 2307100 : AverageReal = 0;
1081 2307100 : AverageImag = 0;
1082 2307100 : StdReal = 1000.0;
1083 2307100 : StdImag = 1000.0;
1084 2307100 : SumWeightReal = 0;
1085 2307100 : SumWeightImag = 0;
1086 :
1087 : // Temporal working variables
1088 2307100 : Complex visibility = Complex(0,0);
1089 4614200 : vector<Double> realVisForMedian;
1090 4614200 : vector<Double> imagVisForMedian;
1091 : // Double realMedian = 0;
1092 : // Double imagMedian = 0;
1093 :
1094 144022956 : for (uInt chan_j=0; chan_j<(uInt) nChannels; ++chan_j)
1095 : {
1096 : // Ignore data point if it is already flagged or weight is <= 0
1097 : // NOTE: In our case visibilities come w/o weights, so we check only vs flags
1098 141715856 : if (flags.getOriginalFlags(pol_k,chan_j,timestep_i)) continue;
1099 :
1100 141334812 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
1101 141334812 : realVisForMedian.push_back(visibility.real());
1102 141334812 : imagVisForMedian.push_back(visibility.imag());
1103 141334812 : SumWeightReal += 1;
1104 141334812 : SumWeightImag += 1;
1105 : }
1106 :
1107 2307100 : if (SumWeightReal)
1108 : {
1109 2301052 : AverageReal = median(realVisForMedian);
1110 2301052 : AverageImag = median(imagVisForMedian);
1111 :
1112 2301052 : noiseVsRef(realVisForMedian,AverageReal);
1113 2301052 : noiseVsRef(imagVisForMedian,AverageImag);
1114 :
1115 2301052 : StdReal = 1.4826 * median(realVisForMedian);
1116 2301052 : StdImag = 1.4826 * median(imagVisForMedian);
1117 : }
1118 :
1119 4614200 : return;
1120 : }
1121 :
1122 : bool
1123 8746 : FlagAgentRFlag::computeAntennaPairFlags(const vi::VisBuffer2 &visBuffer, VisMapper &visibilities,
1124 : FlagMapper &flags,Int /*antenna1*/,Int /*antenna2*/,vector<uInt> &/*rows*/)
1125 : {
1126 : // Set logger origin
1127 8746 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__,WHERE));
1128 :
1129 : // Get flag cube size
1130 : Int nPols,nChannels,nTimesteps;
1131 8746 : visibilities.shape(nPols, nChannels, nTimesteps);
1132 :
1133 : // Make field-spw pair
1134 8746 : auto field = visBuffer.fieldId()(0);
1135 8746 : auto spw = visBuffer.spectralWindows()(0);
1136 8746 : auto field_spw = std::make_pair(field,spw);
1137 :
1138 : // Check if frequency array has to be initialized
1139 8746 : Bool initFreq = false;
1140 :
1141 : // Get noise and scutofff levels
1142 8746 : Double noise = -1;
1143 13058 : if ( (field_spw_noise_map_p.find(field_spw) != field_spw_noise_map_p.end()) and
1144 4312 : field_spw_noise_map_p[field_spw] > 0)
1145 : {
1146 3544 : noise = field_spw_noise_map_p[field_spw];
1147 : }
1148 5202 : else if (noise_p > 0)
1149 : {
1150 24 : noise = noise_p;
1151 : }
1152 5178 : else if (field_spw_noise_histogram_sum_p.find(field_spw) == field_spw_noise_histogram_sum_p.end())
1153 : {
1154 182 : field_spw_noise_histogram_sum_p[field_spw] = vector<Double>(nChannels,0);
1155 182 : field_spw_noise_histogram_counts_p[field_spw] = vector<Double>(nChannels,0);
1156 182 : field_spw_noise_histogram_sum_squares_p[field_spw] = vector<Double>(nChannels,0);
1157 182 : field_spw_frequency_p[field_spw] = vector<Double>(nChannels,0);
1158 182 : if (doflag_p) prepass_p = true;
1159 182 : initFreq = true;
1160 : }
1161 : // Apply scale factor to the flagging threshold
1162 8746 : noise *= noiseScale_p;
1163 :
1164 : // Get cutoff level
1165 8746 : Double scutoff = -1;
1166 13010 : if ((field_spw_scutoff_map_p.find(field_spw) != field_spw_scutoff_map_p.end()) and
1167 4264 : (field_spw_scutoff_map_p[field_spw] > 0))
1168 : {
1169 3496 : scutoff = field_spw_scutoff_map_p[field_spw];
1170 : }
1171 5250 : else if (scutoff_p > 0)
1172 : {
1173 72 : scutoff = scutoff_p;
1174 : }
1175 5178 : else if (field_spw_scutoff_histogram_sum_p.find(field_spw) == field_spw_scutoff_histogram_sum_p.end())
1176 : {
1177 180 : field_spw_scutoff_histogram_sum_p[field_spw] = vector<Double>(nChannels,0);
1178 180 : field_spw_scutoff_histogram_counts_p[field_spw] = vector<Double>(nChannels,0);
1179 180 : field_spw_scutoff_histogram_sum_squares_p[field_spw] = vector<Double>(nChannels,0);
1180 :
1181 180 : if (field_spw_frequency_p.find(field_spw) == field_spw_frequency_p.end())
1182 : {
1183 2 : field_spw_frequency_p[field_spw] = vector<Double>(nChannels,0);
1184 2 : initFreq = true;
1185 : }
1186 :
1187 180 : if (doflag_p) prepass_p = true;
1188 : }
1189 : // Apply scale factor to the flagging threshold
1190 8746 : scutoff *= scutoffScale_p;
1191 :
1192 : // Initialize frequency array has to be initialized
1193 8746 : if (initFreq)
1194 : {
1195 184 : Vector<Double> freqInHz = visBuffer.getFrequencies(0,MFrequency::TOPO);
1196 : // jagonzal (CAS-4312): We have to take into account channel selection for the frequency mapping
1197 9365 : for (uInt channel_idx=0; channel_idx < channelIndex_p.size(); ++channel_idx)
1198 : {
1199 9181 : field_spw_frequency_p[field_spw][channel_idx] = freqInHz[channelIndex_p[channel_idx]]/1E9;
1200 : }
1201 184 : field_spw_frequencies_p[field_spw] = freqInHz[channelIndex_p[0]]/1E9;
1202 : }
1203 :
1204 :
1205 : uInt effectiveNTimeSteps;
1206 8746 : if (nTimesteps > (Int) nTimeSteps_p)
1207 : {
1208 7220 : effectiveNTimeSteps = nTimeSteps_p;
1209 : }
1210 : else
1211 : {
1212 1526 : effectiveNTimeSteps = nTimesteps;
1213 : }
1214 :
1215 8746 : uInt effectiveNTimeStepsDelta = (effectiveNTimeSteps - 1)/2;
1216 :
1217 : // Beginning time range: Move only central point (only for spectral analysis)
1218 : // We set start/stop time with decreasing values to deactivate time analysis
1219 15966 : for (uInt timestep_i=0; timestep_i<effectiveNTimeStepsDelta; ++timestep_i)
1220 : {
1221 : // computeAntennaPairFlagsCore(field_spw,scutoff,0,effectiveNTimeSteps,timestep_i,visibilities,flags);
1222 7220 : computeAntennaPairFlagsCore(field_spw,noise,scutoff,-1,-2,timestep_i,visibilities,flags);
1223 : }
1224 :
1225 591514 : for (uInt timestep_i=effectiveNTimeStepsDelta; timestep_i<nTimesteps-effectiveNTimeStepsDelta; ++timestep_i)
1226 : {
1227 582768 : computeAntennaPairFlagsCore(field_spw,noise,scutoff,timestep_i-effectiveNTimeStepsDelta,timestep_i+effectiveNTimeStepsDelta,timestep_i,visibilities,flags);
1228 : }
1229 :
1230 : // End time range: Move only central point (only for spectral analysis)
1231 : // We set start/stop time with decreasing values to deactivate time analysis
1232 15966 : for (uInt timestep_i=nTimesteps-effectiveNTimeStepsDelta; timestep_i<(uInt) nTimesteps; ++timestep_i)
1233 : {
1234 : // computeAntennaPairFlagsCore(field_spw,scutoff,nTimesteps-effectiveNTimeSteps,nTimesteps-1,timestep_i,visibilities,flags);
1235 7220 : computeAntennaPairFlagsCore(field_spw,noise,scutoff,-1,-2,timestep_i,visibilities,flags);
1236 : }
1237 :
1238 8746 : return false;
1239 : }
1240 :
1241 : void
1242 152 : FlagAgentRFlag::passIntermediate(const vi::VisBuffer2 &visBuffer)
1243 : {
1244 : // Set logger origin
1245 152 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__,WHERE));
1246 :
1247 : // Make field-spw pair
1248 152 : auto field = visBuffer.fieldId()(0);
1249 152 : auto spw = visBuffer.spectralWindows()(0);
1250 152 : auto field_spw = std::make_pair(field,spw);
1251 :
1252 152 : if (field_spw_noise_map_p.find(field_spw) == field_spw_noise_map_p.end() && !noise_p)
1253 : {
1254 152 : Double noise = computeThreshold(field_spw_noise_histogram_sum_p[field_spw],
1255 152 : field_spw_noise_histogram_sum_squares_p[field_spw],
1256 152 : field_spw_noise_histogram_counts_p[field_spw]);
1257 :
1258 152 : field_spw_noise_map_p[field_spw] = noise;
1259 152 : *logger_p << LogIO::DEBUG1 << " field=" << field << " spw=" << spw << " noise=" << noise << LogIO::POST;
1260 : }
1261 :
1262 152 : if (field_spw_scutoff_map_p.find(field_spw) == field_spw_scutoff_map_p.end() && !scutoff_p)
1263 : {
1264 148 : Double scutoff = computeThreshold(field_spw_scutoff_histogram_sum_p[field_spw],
1265 148 : field_spw_scutoff_histogram_sum_squares_p[field_spw],
1266 148 : field_spw_scutoff_histogram_counts_p[field_spw]);
1267 :
1268 148 : field_spw_scutoff_map_p[field_spw] = scutoff;
1269 148 : *logger_p << LogIO::DEBUG1 << " field=" << field << " spw=" << spw << " scutoff=" << scutoff << LogIO::POST;
1270 : }
1271 :
1272 304 : return;
1273 : }
1274 :
1275 : void
1276 152 : FlagAgentRFlag::passFinal(const vi::VisBuffer2 &visBuffer)
1277 : {
1278 :
1279 : // Make field-spw pair
1280 152 : auto field = visBuffer.fieldId()(0);
1281 152 : auto spw = visBuffer.spectralWindows()(0);
1282 152 : auto field_spw = std::make_pair(field,spw);
1283 152 : if (user_field_spw_noise_map_p.find(field_spw) == user_field_spw_noise_map_p.end())
1284 : {
1285 152 : field_spw_noise_map_p.erase(field_spw);
1286 152 : field_spw_noise_histogram_sum_p.erase(field_spw);
1287 152 : field_spw_noise_histogram_sum_squares_p.erase(field_spw);
1288 152 : field_spw_noise_histogram_counts_p.erase(field_spw);
1289 : }
1290 :
1291 152 : if (user_field_spw_scutoff_map_p.find(field_spw) == user_field_spw_scutoff_map_p.end())
1292 : {
1293 152 : field_spw_scutoff_map_p.erase(field_spw);
1294 152 : field_spw_scutoff_histogram_sum_p.erase(field_spw);
1295 152 : field_spw_scutoff_histogram_sum_squares_p.erase(field_spw);
1296 152 : field_spw_scutoff_histogram_counts_p.erase(field_spw);
1297 : }
1298 :
1299 304 : return;
1300 : }
1301 :
1302 : } //# NAMESPACE CASA - END
1303 :
1304 :
|
