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 46 : FlagAgentRFlag::FlagAgentRFlag(FlagDataHandler *dh, Record config, Bool writePrivateFlagCube, Bool flag):
32 46 : FlagAgentBase(dh,config,ANTENNA_PAIRS,writePrivateFlagCube,flag)
33 : {
34 46 : setAgentParameters(config);
35 :
36 : // Request pre-loading spw
37 46 : flagDataHandler_p->preLoadColumn(VisBufferComponent2::SpectralWindows);
38 : // flagDataHandler_p->preLoadColumn(vi::Freq);
39 :
40 : // Initialize parameters for robust stats (spectral analysis)
41 46 : nIterationsRobust_p = 12;
42 46 : thresholdRobust_p = vector<Double>(nIterationsRobust_p);
43 46 : thresholdRobust_p[0] = 6.0;
44 46 : thresholdRobust_p[1] = 5.0;
45 46 : thresholdRobust_p[2] = 4.0;
46 46 : thresholdRobust_p[3] = 3.6;
47 46 : thresholdRobust_p[4] = 3.2;
48 46 : thresholdRobust_p[5] = 3.0;
49 46 : thresholdRobust_p[6] = 2.7;
50 46 : thresholdRobust_p[7] = 2.6;
51 46 : thresholdRobust_p[8] = 2.5;
52 46 : thresholdRobust_p[9] = 2.5;
53 46 : thresholdRobust_p[10] = 2.5;
54 46 : thresholdRobust_p[11] = 3.5;
55 46 : }
56 :
57 92 : FlagAgentRFlag::~FlagAgentRFlag()
58 : {
59 : // Compiler automagically calls FlagAgentBase::~FlagAgentBase()
60 92 : }
61 :
62 46 : void FlagAgentRFlag::setAgentParameters(Record config)
63 : {
64 46 : 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 46 : String display("none");
70 46 : exists = config.fieldNumber ("display");
71 46 : if (exists >= 0)
72 : {
73 46 : if( config.type(exists) != TpString )
74 : {
75 0 : throw( AipsError ( "Parameter 'display' must be of type 'string'" ) );
76 : }
77 46 : display = config.asString("display");
78 46 : *logger_p << LogIO::NORMAL << " display is: " << display << LogIO::POST;
79 : }
80 :
81 : // Determine if we have to generate plot reports.
82 46 : if ( (display == String("report")) or (display == String("both")) )
83 : {
84 0 : doplot_p = true;
85 : }
86 : else
87 : {
88 46 : doplot_p = false;
89 : }
90 :
91 : // Do we actually flag, or just calculate thresholds without applying flags?
92 46 : Bool writeflags(display != "report");
93 46 : exists = config.fieldNumber ("writeflags");
94 46 : if (exists >= 0)
95 : {
96 46 : if( config.type(exists) != TpBool )
97 : {
98 0 : throw( AipsError ( "Parameter 'writeflags' must be of type 'bool'" ) );
99 : }
100 46 : writeflags = config.asBool("writeflags");
101 46 : *logger_p << LogIO::NORMAL << " writeflags is: " << writeflags << LogIO::POST;
102 : }
103 :
104 46 : if( (writeflags == true) or (display == String("data")) or (display == String("both")) )
105 : {
106 32 : doflag_p = true;
107 32 : *logger_p << LogIO::NORMAL << " (writeflags OR display)=(" << writeflags << "," << display << "), will apply flags on modified flag cube " << LogIO::POST;
108 : }
109 : else
110 : {
111 14 : doflag_p = false;
112 : }
113 :
114 : // AIPS RFlag FPARM(1)
115 46 : exists = config.fieldNumber ("winsize");
116 46 : if (exists >= 0)
117 : {
118 42 : if( config.type(exists) != TpInt )
119 : {
120 0 : throw( AipsError ( "Parameter 'winsize' must be of type 'Int'" ) );
121 : }
122 :
123 42 : nTimeSteps_p = config.asuInt("winsize");
124 : }
125 : else
126 : {
127 4 : nTimeSteps_p = 3;
128 : }
129 :
130 46 : *logger_p << logLevel_p << " winsize is " << nTimeSteps_p << LogIO::POST;
131 :
132 : // AIPS RFlag FPARM(5)
133 46 : exists = config.fieldNumber ("spectralmax");
134 46 : if (exists >= 0)
135 : {
136 42 : 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 42 : spectralmax_p = config.asDouble("spectralmax");
142 : }
143 : else
144 : {
145 4 : spectralmax_p = 1E6;
146 : }
147 :
148 46 : *logger_p << logLevel_p << " spectralmax is " << spectralmax_p << LogIO::POST;
149 :
150 : // AIPS RFlag FPARM(6)
151 46 : exists = config.fieldNumber ("spectralmin");
152 46 : if (exists >= 0)
153 : {
154 42 : 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 42 : spectralmin_p = config.asDouble("spectralmin");
160 : }
161 : else
162 : {
163 4 : spectralmin_p = 0;
164 : }
165 :
166 46 : *logger_p << logLevel_p << " spectralmin is " << spectralmin_p << LogIO::POST;
167 :
168 : // AIPS RFlag OPTYPE (type of operation for spectral analysis)
169 46 : String optype("MEDIAN");
170 46 : optype_p = MEDIAN;
171 46 : spectralAnalysis_p = &FlagAgentRFlag::simpleMedian;
172 46 : exists = config.fieldNumber ("optype");
173 46 : 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 46 : exists = config.fieldNumber ("timedevscale");
215 46 : if (exists >= 0)
216 : {
217 42 : 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 42 : noiseScale_p = config.asDouble("timedevscale");
223 : }
224 : else
225 : {
226 4 : noiseScale_p = 5;
227 : }
228 :
229 46 : *logger_p << logLevel_p << " timedevscale is " << noiseScale_p << LogIO::POST;
230 46 : const auto thresholds_note = "Note: not applying flags, simply calculating "
231 : "thresholds (action not 'apply')";
232 46 : if (not writeflags) {
233 14 : *logger_p << logLevel_p << " " << thresholds_note <<
234 14 : " so timedevscale is not used." << LogIO::POST;
235 : }
236 :
237 : // AIPS RFlag FPARM(10)
238 46 : exists = config.fieldNumber ("freqdevscale");
239 46 : if (exists >= 0)
240 : {
241 44 : 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 44 : scutoffScale_p = config.asDouble("freqdevscale");
247 : }
248 : else
249 : {
250 2 : scutoffScale_p = 5;
251 : }
252 :
253 46 : *logger_p << logLevel_p << " freqdevscale is " << scutoffScale_p << LogIO::POST;
254 46 : if (not writeflags) {
255 14 : *logger_p << logLevel_p << " " << thresholds_note <<
256 14 : "so freqdevscale is not used." << LogIO::POST;
257 : }
258 :
259 : // timedev - Matrix for time analysis deviation thresholds - (old AIPS RFlag FPARM(3)/NOISE)
260 46 : noise_p = 0;
261 46 : exists = config.fieldNumber ("timedev");
262 46 : bool nonempty = exists >= 0;
263 : // Special empty values
264 46 : if (nonempty) {
265 : // when no value is given to timedev, it 'exists' but as an empty string
266 43 : if (casacore::TpString == config.type(exists)) {
267 12 : auto value = config.asString(exists);
268 12 : if (0 == value.length()) {
269 11 : nonempty = false;
270 : }
271 12 : }
272 : // "timedev=[]" arrives here as an empty array of bool!
273 31 : else if (casacore::TpArrayBool == config.type(exists)) {
274 21 : auto tdev = config.asArrayBool(RecordFieldId("timedev"));
275 21 : if (0 == tdev.size()) {
276 21 : nonempty = false;
277 : }
278 21 : }
279 : }
280 46 : if (nonempty)
281 : {
282 11 : 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 9 : else if( config.type(exists) == casacore::TpArrayDouble || config.type(exists) == casacore::TpArrayFloat || config.type(exists) == casacore::TpArrayInt)
289 : {
290 8 : Matrix<Double> timedev = config.asArrayDouble( RecordFieldId("timedev") );
291 8 : if(timedev.ncolumn()==3)
292 : {
293 8 : *logger_p << logLevel_p << " timedev [field,spw,dev] is "
294 8 : << timedev << LogIO::POST;
295 :
296 8 : IPosition shape = timedev.shape();
297 8 : uInt nDevs = shape[0];
298 24 : for(uInt dev_i=0; dev_i<nDevs; ++dev_i)
299 : {
300 16 : auto field_spw = std::make_pair((Int)timedev(dev_i, 0),
301 16 : (Int)timedev(dev_i, 1));
302 16 : field_spw_noise_map_p[field_spw] = timedev(dev_i, 2);
303 16 : user_field_spw_noise_map_p[field_spw] = true;
304 16 : *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 8 : }
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 8 : }
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 0 : 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 35 : *logger_p << logLevel_p << "No timedev value given. Will Use automatically "
328 35 : "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 46 : scutoff_p = 0;
335 46 : exists = config.fieldNumber ("freqdev");
336 46 : nonempty = exists >= 0;
337 : // Special empty values
338 46 : if (nonempty) {
339 : // when no value is given to freqdev, it 'exists' but as an empty string
340 43 : if (casacore::TpString == config.type(exists)) {
341 12 : auto value = config.asString(exists);
342 12 : if (0 == value.length()) {
343 11 : nonempty = false;
344 : }
345 12 : }
346 : // "freqdev=[]" arrives here as an empty array of bool!
347 31 : else if (casacore::TpArrayBool == config.type(exists)) {
348 23 : auto fdev = config.asArrayBool(RecordFieldId("freqdev"));
349 23 : if (0 == fdev.size()) {
350 23 : nonempty = false;
351 : }
352 23 : }
353 : }
354 46 : if (nonempty)
355 : {
356 9 : 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 7 : else if( config.type(exists) == casacore::TpArrayDouble || config.type(exists) == casacore::TpArrayFloat || config.type(exists) == casacore::TpArrayInt)
363 : {
364 6 : Matrix<Double> freqdev = config.asArrayDouble( RecordFieldId("freqdev") );
365 6 : if(freqdev.ncolumn()==3)
366 : {
367 6 : *logger_p << logLevel_p << " freqdev [field,spw,dev] is " <<
368 6 : freqdev << LogIO::POST;
369 :
370 6 : IPosition shape = freqdev.shape();
371 6 : uInt nDevs = shape[0];
372 18 : for(uInt dev_i=0; dev_i<nDevs; ++dev_i)
373 : {
374 12 : auto field_spw = std::make_pair(static_cast<Int>(freqdev(dev_i, 0)),
375 12 : static_cast<Int>(freqdev(dev_i, 1)));
376 12 : field_spw_scutoff_map_p[field_spw] = freqdev(dev_i, 2);
377 12 : user_field_spw_scutoff_map_p[field_spw] = true;
378 12 : *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 6 : }
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 6 : }
386 : else
387 : {
388 1 : if (writeflags) {
389 0 : 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 37 : *logger_p << logLevel_p << "No freqdev value given. Will Use automatically "
398 37 : "computed values if applying flags." << LogIO::POST;
399 : // scutoff initialized to zero above.
400 : }
401 :
402 92 : return;
403 46 : }
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 4499000 : void FlagAgentRFlag::noiseVsRef(vector<Double> &data, Double ref)
430 : {
431 280569968 : for (size_t index = 0; index < data.size(); ++index)
432 : {
433 276070968 : data[index] -= ref;
434 : }
435 :
436 4499000 : return;
437 : }
438 :
439 8998692 : Double FlagAgentRFlag::median(vector<Double> &data)
440 : {
441 : Double med;
442 8998692 : vector<Double> datacopy = data;
443 8998692 : sort(data.begin(),data.end());
444 :
445 8998692 : if (data.size() == 0)
446 : {
447 0 : med = 0;
448 : }
449 8998692 : else if (data.size() % 2 == 1)
450 : {
451 34516 : med = data[(data.size()-1)/2];
452 : }
453 : else
454 : {
455 8964176 : med = 0.5*(data[data.size()/2] + data[(data.size()/2)-1]);
456 : }
457 :
458 8998692 : return med;
459 8998692 : }
460 :
461 346 : 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 346 : vector<Double> samplesForMedian(data.size(),0);
468 17304 : for (size_t index = 0; index < data.size(); ++index)
469 : {
470 16958 : if (counts[index] > 0)
471 : {
472 16956 : avg = data[index]/counts[index];
473 16956 : samplesForMedian[index] = avg;
474 : }
475 : }
476 :
477 : // Compute median
478 346 : Double med = median(samplesForMedian);
479 :
480 : // Produce samples for median absolute deviation
481 346 : vector<Double> samplesForMad(samplesForMedian.size(),0);
482 17304 : for (size_t index = 0; index < samplesForMedian.size(); ++index)
483 : {
484 16958 : samplesForMad[index] = abs(samplesForMedian[index] - med);
485 : }
486 :
487 : // Compute median absolute deviation
488 346 : Double mad = median(samplesForMad);
489 :
490 346 : return (med + 1.4826*mad);
491 346 : }
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 46 : FlagReport FlagAgentRFlag::getReport()
502 : {
503 92 : FlagReport totalRep(String("list"),agentName_p);
504 :
505 46 : if (doflag_p) {
506 32 : return totalRep;
507 : }
508 :
509 28 : FlagReport plotRepCont(String("list"),agentName_p);
510 :
511 : // Calculate thresholds
512 28 : generateThresholds( field_spw_noise_histogram_sum_p,
513 14 : field_spw_noise_histogram_sum_squares_p,
514 14 : field_spw_noise_histogram_counts_p,
515 14 : field_spw_noise_map_p,
516 : "Time analysis");
517 28 : generateThresholds( field_spw_scutoff_histogram_sum_p,
518 14 : field_spw_scutoff_histogram_sum_squares_p,
519 14 : field_spw_scutoff_histogram_counts_p,
520 14 : field_spw_scutoff_map_p,
521 : "Spectral analysis");
522 :
523 : // Threshold reports (should be returned if params were calculated)
524 14 : Record threshList;
525 14 : Int nEntriesNoise = field_spw_noise_map_p.size();
526 14 : Int nEntriesScutoff = field_spw_scutoff_map_p.size();
527 :
528 14 : if(nEntriesNoise > 0 || nEntriesScutoff > 0)
529 : {
530 14 : Matrix<Double> timedev(nEntriesNoise,3), freqdev(nEntriesScutoff,3);
531 14 : Int threshCountTime = 0, threshCountFreq = 0;
532 14 : pair<Int,Int> field_spw;
533 14 : for (auto spw_field_iter = field_spw_noise_map_p.begin();
534 42 : spw_field_iter != field_spw_noise_map_p.end();
535 28 : ++spw_field_iter)
536 : {
537 28 : field_spw = spw_field_iter->first;
538 :
539 28 : timedev(threshCountTime,0) = field_spw.first;
540 28 : timedev(threshCountTime,1) = field_spw.second;
541 28 : timedev(threshCountTime,2) = field_spw_noise_map_p[field_spw];
542 28 : ++threshCountTime;
543 : }
544 :
545 14 : for (auto spw_field_iter = field_spw_scutoff_map_p.begin();
546 44 : spw_field_iter != field_spw_scutoff_map_p.end();
547 30 : ++spw_field_iter)
548 : {
549 30 : field_spw = spw_field_iter->first;
550 :
551 30 : freqdev(threshCountFreq,0) = field_spw.first;
552 30 : freqdev(threshCountFreq,1) = field_spw.second;
553 30 : freqdev(threshCountFreq,2) = field_spw_scutoff_map_p[field_spw];
554 30 : ++threshCountFreq;
555 : }
556 :
557 14 : threshList.define( RecordFieldId("timedev") , timedev );
558 14 : threshList.define( RecordFieldId("freqdev") , freqdev );
559 :
560 28 : FlagReport returnThresh("rflag",agentName_p, threshList);
561 14 : totalRep.addReport(returnThresh);
562 14 : }
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 14 : 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 0 : }
591 : }
592 :
593 14 : return totalRep;
594 14 : }
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 0 : }
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 0 : }
736 :
737 0 : return totalReport;
738 0 : }
739 :
740 28 : 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 28 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__,WHERE));
748 :
749 : // First of all determine each SPW frequency in order to produce ordered vectors
750 28 : pair<Int,Int> current_field_spw;
751 86 : for (auto field_spw_iter = data.begin(); field_spw_iter != data.end();
752 58 : ++field_spw_iter)
753 : {
754 58 : current_field_spw = field_spw_iter->first;
755 :
756 : // Compute threshold
757 58 : threshold[current_field_spw] = computeThreshold(data[current_field_spw],
758 58 : dataSquared[current_field_spw],
759 58 : counts[current_field_spw]);
760 :
761 : // Log info
762 58 : auto nfreq = field_spw_frequency_p[current_field_spw].size();
763 58 : auto freqStart = field_spw_frequency_p[current_field_spw][0];
764 58 : auto freqEnd = field_spw_frequency_p[current_field_spw][nfreq-1];
765 58 : *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 58 : << threshold[current_field_spw] << LogIO::POST;
770 : }
771 :
772 56 : return;
773 : }
774 :
775 581296 : 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 581296 : visibilities.shape(nPols, nChannels, nTimesteps);
787 :
788 : // Declare variables
789 581296 : Complex visibility;
790 581296 : Double SumWeight = 0;
791 581296 : Double SumWeightReal = 0;
792 581296 : Double SumWeightImag = 0;
793 581296 : Double StdTotal = 0;
794 581296 : Double SumReal = 0;
795 581296 : Double SumRealSquare = 0;
796 581296 : Double AverageReal = 0;
797 581296 : Double StdReal = 0;
798 581296 : Double SumImag = 0;
799 581296 : Double SumImagSquare = 0;
800 581296 : Double AverageImag = 0;
801 581296 : Double StdImag = 0;
802 581296 : Double deviationReal = 0;
803 581296 : 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 581296 : if ( (noise <= 0) or ((noise >= 0) and (doflag_p == true) and (prepass_p == false)) )
808 : {
809 36161496 : for (uInt chan_j=0; chan_j<(uInt)nChannels; ++chan_j)
810 : {
811 : // Compute variance
812 173617852 : for (uInt pol_k=0;pol_k<(uInt)nPols; ++pol_k)
813 : {
814 138035504 : SumWeight = 0;
815 138035504 : StdTotal = 0;
816 138035504 : SumReal = 0;
817 138035504 : SumRealSquare = 0;
818 138035504 : AverageReal = 0;
819 138035504 : StdReal = 0;
820 138035504 : SumImag = 0;
821 138035504 : SumImagSquare = 0;
822 138035504 : AverageImag = 0;
823 138035504 : StdImag = 0;
824 :
825 545594640 : 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 407559136 : if (flags.getOriginalFlags(pol_k,chan_j,timestep_i)) continue;
830 :
831 407559116 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
832 407559116 : SumWeight += 1;
833 407559116 : SumReal += visibility.real();
834 407559116 : SumRealSquare += visibility.real()*visibility.real();
835 407559116 : SumImag += visibility.imag();
836 407559116 : 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 138035504 : if (SumWeight > 0)
842 : {
843 135856780 : AverageReal = SumReal / SumWeight;
844 135856780 : SumRealSquare = SumRealSquare / SumWeight;
845 135856780 : StdReal = SumRealSquare - AverageReal * AverageReal;
846 :
847 135856780 : AverageImag = SumImag / SumWeight;
848 135856780 : SumImagSquare = SumImagSquare / SumWeight;
849 135856780 : 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 135856780 : StdReal = StdReal > 0? StdReal:0;
856 135856780 : StdImag = StdImag > 0? StdImag:0;
857 135856780 : 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 135856780 : if (noise <= 0)
863 : {
864 68062508 : field_spw_noise_histogram_counts_p[spw_field][chan_j] += 1;
865 68062508 : field_spw_noise_histogram_sum_p[spw_field][chan_j] += StdTotal;
866 68062508 : field_spw_noise_histogram_sum_squares_p[spw_field][chan_j] += StdTotal*StdTotal;
867 : }
868 :
869 135856780 : if (noise >=0 and StdTotal > noise)
870 : {
871 9569010 : for (uInt timestep_i=timeStart; timestep_i<=timeStop; ++timestep_i)
872 : {
873 7176728 : if (!flags.getModifiedFlags(pol_k,chan_j,timestep_i))
874 : {
875 2991528 : flags.setModifiedFlags(pol_k,chan_j,timestep_i);
876 2991528 : visBufferFlags_p += 1;
877 : }
878 : }
879 : }
880 : }
881 : }
882 : }
883 : }
884 :
885 : // Spectral analysis: Fix timestep/polarization and compute stats with all channels
886 581296 : if ( (scutoff <= 0) or ((scutoff >= 0) and (doflag_p == true) and (prepass_p == false)) )
887 : {
888 1158296 : for (uInt timestep_i=centralTime; timestep_i<=centralTime; ++timestep_i)
889 : {
890 2828648 : for (uInt pol_k=0; pol_k<(uInt) nPols; ++pol_k)
891 : {
892 :
893 2249500 : (*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 : flags);
904 :
905 2249500 : if (scutoff <= 0)
906 : {
907 70313190 : 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 69171656 : if (flags.getOriginalFlags(pol_k,chan_j,timestep_i)) continue;
912 :
913 69171636 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
914 :
915 69171636 : if (SumWeightReal > 0)
916 : {
917 69171636 : deviationReal = abs(visibility.real()-AverageReal);
918 69171636 : field_spw_scutoff_histogram_counts_p[spw_field][chan_j] += 1;
919 69171636 : field_spw_scutoff_histogram_sum_p[spw_field][chan_j] += deviationReal;
920 69171636 : field_spw_scutoff_histogram_sum_squares_p[spw_field][chan_j] += deviationReal*deviationReal;
921 : }
922 :
923 69171636 : if (SumWeightImag > 0)
924 : {
925 69171636 : deviationImag = abs(visibility.imag()-AverageImag);
926 69171636 : field_spw_scutoff_histogram_counts_p[spw_field][chan_j] += 1;
927 69171636 : field_spw_scutoff_histogram_sum_p[spw_field][chan_j] += deviationImag;
928 69171636 : field_spw_scutoff_histogram_sum_squares_p[spw_field][chan_j] += deviationImag*deviationImag;
929 : }
930 : }
931 : }
932 :
933 2249500 : if (scutoff >=0)
934 : {
935 : // Flag all channels?
936 1107966 : if ( (StdReal > spectralmax_p) or
937 1107966 : (StdImag > spectralmax_p) or
938 1107966 : (StdReal < spectralmin_p) or
939 1107966 : (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 : }
949 0 : }
950 : // Check each channel separately vs the scutoff level
951 : else
952 : {
953 69971814 : for (uInt chan_j=0; chan_j<(uInt) nChannels; ++chan_j)
954 : {
955 68863848 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
956 135729861 : if ( (abs(visibility.real()-AverageReal)>scutoff) or
957 66866013 : (abs(visibility.imag()-AverageImag)>scutoff) )
958 : {
959 2834772 : if (!flags.getModifiedFlags(pol_k,chan_j,timestep_i))
960 : {
961 961235 : flags.setModifiedFlags(pol_k,chan_j,timestep_i);
962 961235 : visBufferFlags_p += 1;
963 : }
964 : }
965 : }
966 : }
967 : }
968 : }
969 : }
970 : }
971 :
972 1162592 : 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 2249500 : 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 2249500 : AverageReal = 0;
1081 2249500 : AverageImag = 0;
1082 2249500 : StdReal = 1000.0;
1083 2249500 : StdImag = 1000.0;
1084 2249500 : SumWeightReal = 0;
1085 2249500 : SumWeightImag = 0;
1086 :
1087 : // Temporal working variables
1088 2249500 : Complex visibility = Complex(0,0);
1089 2249500 : vector<Double> realVisForMedian;
1090 2249500 : vector<Double> imagVisForMedian;
1091 : // Double realMedian = 0;
1092 : // Double imagMedian = 0;
1093 :
1094 140285004 : 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 138035504 : if (flags.getOriginalFlags(pol_k,chan_j,timestep_i)) continue;
1099 :
1100 138035484 : visibility = visibilities.correlationProduct(pol_k,chan_j,timestep_i);
1101 138035484 : realVisForMedian.push_back(visibility.real());
1102 138035484 : imagVisForMedian.push_back(visibility.imag());
1103 138035484 : SumWeightReal += 1;
1104 138035484 : SumWeightImag += 1;
1105 : }
1106 :
1107 2249500 : if (SumWeightReal)
1108 : {
1109 2249500 : AverageReal = median(realVisForMedian);
1110 2249500 : AverageImag = median(imagVisForMedian);
1111 :
1112 2249500 : noiseVsRef(realVisForMedian,AverageReal);
1113 2249500 : noiseVsRef(imagVisForMedian,AverageImag);
1114 :
1115 2249500 : StdReal = 1.4826 * median(realVisForMedian);
1116 2249500 : StdImag = 1.4826 * median(imagVisForMedian);
1117 : }
1118 :
1119 4499000 : return;
1120 2249500 : }
1121 :
1122 : bool
1123 7090 : FlagAgentRFlag::computeAntennaPairFlags(const vi::VisBuffer2 &visBuffer, VisMapper &visibilities,
1124 : FlagMapper &flags,Int /*antenna1*/,Int /*antenna2*/,vector<uInt> &/*rows*/)
1125 : {
1126 : // Set logger origin
1127 7090 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__,WHERE));
1128 :
1129 : // Get flag cube size
1130 : Int nPols,nChannels,nTimesteps;
1131 7090 : visibilities.shape(nPols, nChannels, nTimesteps);
1132 :
1133 : // Make field-spw pair
1134 7090 : auto field = visBuffer.fieldId()(0);
1135 7090 : auto spw = visBuffer.spectralWindows()(0);
1136 7090 : auto field_spw = std::make_pair(field,spw);
1137 :
1138 : // Check if frequency array has to be initialized
1139 7090 : Bool initFreq = false;
1140 :
1141 : // Get noise and scutofff levels
1142 7090 : Double noise = -1;
1143 10550 : if ( (field_spw_noise_map_p.find(field_spw) != field_spw_noise_map_p.end()) and
1144 3460 : field_spw_noise_map_p[field_spw] > 0)
1145 : {
1146 3448 : noise = field_spw_noise_map_p[field_spw];
1147 : }
1148 3642 : else if (noise_p > 0)
1149 : {
1150 24 : noise = noise_p;
1151 : }
1152 3618 : else if (field_spw_noise_histogram_sum_p.find(field_spw) == field_spw_noise_histogram_sum_p.end())
1153 : {
1154 174 : field_spw_noise_histogram_sum_p[field_spw] = vector<Double>(nChannels,0);
1155 174 : field_spw_noise_histogram_counts_p[field_spw] = vector<Double>(nChannels,0);
1156 174 : field_spw_noise_histogram_sum_squares_p[field_spw] = vector<Double>(nChannels,0);
1157 174 : field_spw_frequency_p[field_spw] = vector<Double>(nChannels,0);
1158 174 : if (doflag_p) prepass_p = true;
1159 174 : initFreq = true;
1160 : }
1161 : // Apply scale factor to the flagging threshold
1162 7090 : noise *= noiseScale_p;
1163 :
1164 : // Get cutoff level
1165 7090 : Double scutoff = -1;
1166 10502 : if ((field_spw_scutoff_map_p.find(field_spw) != field_spw_scutoff_map_p.end()) and
1167 3412 : (field_spw_scutoff_map_p[field_spw] > 0))
1168 : {
1169 3400 : scutoff = field_spw_scutoff_map_p[field_spw];
1170 : }
1171 3690 : else if (scutoff_p > 0)
1172 : {
1173 72 : scutoff = scutoff_p;
1174 : }
1175 3618 : else if (field_spw_scutoff_histogram_sum_p.find(field_spw) == field_spw_scutoff_histogram_sum_p.end())
1176 : {
1177 172 : field_spw_scutoff_histogram_sum_p[field_spw] = vector<Double>(nChannels,0);
1178 172 : field_spw_scutoff_histogram_counts_p[field_spw] = vector<Double>(nChannels,0);
1179 172 : field_spw_scutoff_histogram_sum_squares_p[field_spw] = vector<Double>(nChannels,0);
1180 :
1181 172 : 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 172 : if (doflag_p) prepass_p = true;
1188 : }
1189 : // Apply scale factor to the flagging threshold
1190 7090 : scutoff *= scutoffScale_p;
1191 :
1192 : // Initialize frequency array has to be initialized
1193 7090 : if (initFreq)
1194 : {
1195 176 : 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 8847 : for (uInt channel_idx=0; channel_idx < channelIndex_p.size(); ++channel_idx)
1198 : {
1199 8671 : field_spw_frequency_p[field_spw][channel_idx] = freqInHz[channelIndex_p[channel_idx]]/1E9;
1200 : }
1201 176 : field_spw_frequencies_p[field_spw] = freqInHz[channelIndex_p[0]]/1E9;
1202 176 : }
1203 :
1204 :
1205 : uInt effectiveNTimeSteps;
1206 7090 : if (nTimesteps > (Int) nTimeSteps_p)
1207 : {
1208 7076 : effectiveNTimeSteps = nTimeSteps_p;
1209 : }
1210 : else
1211 : {
1212 14 : effectiveNTimeSteps = nTimesteps;
1213 : }
1214 :
1215 7090 : 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 14166 : for (uInt timestep_i=0; timestep_i<effectiveNTimeStepsDelta; ++timestep_i)
1220 : {
1221 : // computeAntennaPairFlagsCore(field_spw,scutoff,0,effectiveNTimeSteps,timestep_i,visibilities,flags);
1222 7076 : computeAntennaPairFlagsCore(field_spw,noise,scutoff,-1,-2,timestep_i,visibilities,flags);
1223 : }
1224 :
1225 574234 : for (uInt timestep_i=effectiveNTimeStepsDelta; timestep_i<nTimesteps-effectiveNTimeStepsDelta; ++timestep_i)
1226 : {
1227 567144 : 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 14166 : 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 7076 : computeAntennaPairFlagsCore(field_spw,noise,scutoff,-1,-2,timestep_i,visibilities,flags);
1236 : }
1237 :
1238 7090 : return false;
1239 : }
1240 :
1241 : void
1242 146 : FlagAgentRFlag::passIntermediate(const vi::VisBuffer2 &visBuffer)
1243 : {
1244 : // Set logger origin
1245 146 : logger_p->origin(LogOrigin(agentName_p,__FUNCTION__,WHERE));
1246 :
1247 : // Make field-spw pair
1248 146 : auto field = visBuffer.fieldId()(0);
1249 146 : auto spw = visBuffer.spectralWindows()(0);
1250 146 : auto field_spw = std::make_pair(field,spw);
1251 :
1252 146 : if (field_spw_noise_map_p.find(field_spw) == field_spw_noise_map_p.end() && !noise_p)
1253 : {
1254 146 : Double noise = computeThreshold(field_spw_noise_histogram_sum_p[field_spw],
1255 146 : field_spw_noise_histogram_sum_squares_p[field_spw],
1256 146 : field_spw_noise_histogram_counts_p[field_spw]);
1257 :
1258 146 : field_spw_noise_map_p[field_spw] = noise;
1259 146 : *logger_p << LogIO::DEBUG1 << " field=" << field << " spw=" << spw << " noise=" << noise << LogIO::POST;
1260 : }
1261 :
1262 146 : if (field_spw_scutoff_map_p.find(field_spw) == field_spw_scutoff_map_p.end() && !scutoff_p)
1263 : {
1264 142 : Double scutoff = computeThreshold(field_spw_scutoff_histogram_sum_p[field_spw],
1265 142 : field_spw_scutoff_histogram_sum_squares_p[field_spw],
1266 142 : field_spw_scutoff_histogram_counts_p[field_spw]);
1267 :
1268 142 : field_spw_scutoff_map_p[field_spw] = scutoff;
1269 142 : *logger_p << LogIO::DEBUG1 << " field=" << field << " spw=" << spw << " scutoff=" << scutoff << LogIO::POST;
1270 : }
1271 :
1272 292 : return;
1273 : }
1274 :
1275 : void
1276 146 : FlagAgentRFlag::passFinal(const vi::VisBuffer2 &visBuffer)
1277 : {
1278 :
1279 : // Make field-spw pair
1280 146 : auto field = visBuffer.fieldId()(0);
1281 146 : auto spw = visBuffer.spectralWindows()(0);
1282 146 : auto field_spw = std::make_pair(field,spw);
1283 146 : if (user_field_spw_noise_map_p.find(field_spw) == user_field_spw_noise_map_p.end())
1284 : {
1285 146 : field_spw_noise_map_p.erase(field_spw);
1286 146 : field_spw_noise_histogram_sum_p.erase(field_spw);
1287 146 : field_spw_noise_histogram_sum_squares_p.erase(field_spw);
1288 146 : field_spw_noise_histogram_counts_p.erase(field_spw);
1289 : }
1290 :
1291 146 : if (user_field_spw_scutoff_map_p.find(field_spw) == user_field_spw_scutoff_map_p.end())
1292 : {
1293 146 : field_spw_scutoff_map_p.erase(field_spw);
1294 146 : field_spw_scutoff_histogram_sum_p.erase(field_spw);
1295 146 : field_spw_scutoff_histogram_sum_squares_p.erase(field_spw);
1296 146 : field_spw_scutoff_histogram_counts_p.erase(field_spw);
1297 : }
1298 :
1299 292 : return;
1300 : }
1301 :
1302 : } //# NAMESPACE CASA - END
1303 :
1304 :
|
