Line data Source code
1 : ///# Flagger.cc: this defines Flagger
2 : //# Copyright (C) 2000,2001,2002
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 : #include <casacore/casa/Arrays/Vector.h>
28 : #include <casacore/casa/Arrays/ArrayMath.h>
29 : #include <casacore/casa/Arrays/ArrayLogical.h>
30 : #include <casacore/ms/MeasurementSets/MSColumns.h>
31 : #include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
32 : #include <casacore/casa/BasicSL/Complex.h>
33 : #include <casacore/measures/Measures/Stokes.h>
34 : #include <casacore/casa/Utilities/Regex.h>
35 : #include <casacore/casa/OS/HostInfo.h>
36 : #include <flagging/Flagging/Flagger.h>
37 : #include <flagging/Flagging/RFAFlagExaminer.h>
38 : #include <flagging/Flagging/RFAMedianClip.h>
39 : #include <flagging/Flagging/RFASpectralRej.h>
40 : #include <flagging/Flagging/RFASelector.h>
41 : #include <flagging/Flagging/RFAUVBinner.h>
42 : #include <flagging/Flagging/RFATimeFreqCrop.h>
43 : #include <msvis/MSVis/VisibilityIterator.h>
44 : #include <msvis/MSVis/VisBuffer.h>
45 : #include <casacore/casa/System/ProgressMeter.h>
46 : #include <stdarg.h>
47 :
48 : #include <casacore/tables/Tables/Table.h>
49 : #include <casacore/tables/TaQL/TableParse.h>
50 : #include <casacore/tables/Tables/TableRecord.h>
51 : #include <casacore/tables/Tables/TableDesc.h>
52 : #include <casacore/tables/Tables/TableLock.h>
53 : #include <casacore/tables/Tables/SetupNewTab.h>
54 :
55 : #include <casacore/tables/TaQL/ExprNode.h>
56 : #include <msvis/MSVis/VisSet.h>
57 : #include <msvis/MSVis/VisSetUtil.h>
58 :
59 : #include <casacore/measures/Measures/Stokes.h>
60 : #include <casacore/casa/Quanta/UnitMap.h>
61 : #include <casacore/casa/Quanta/UnitVal.h>
62 : #include <casacore/casa/Quanta/MVAngle.h>
63 : #include <casacore/measures/Measures/MDirection.h>
64 : #include <casacore/measures/Measures/MPosition.h>
65 : #include <casacore/casa/Quanta/MVEpoch.h>
66 : #include <casacore/measures/Measures/MEpoch.h>
67 : #include <casacore/measures/Measures/MeasTable.h>
68 :
69 : #include <flagging/Flagging/RFANewMedianClip.h>
70 :
71 : #include <algorithm>
72 :
73 : #include <sstream>
74 :
75 : using namespace casacore;
76 : namespace casa {
77 :
78 : const bool Flagger::dbg = false;
79 :
80 : LogIO Flagger::os( LogOrigin("Flagger") );
81 : static char str[1024];
82 : // uInt debug_ifr=9999,debug_itime=9999;
83 :
84 :
85 : // -----------------------------------------------------------------------
86 : // Default Constructor
87 : // -----------------------------------------------------------------------
88 0 : Flagger::Flagger ():mssel_p(0), vs_p(0), scan_looping(false)
89 : {
90 0 : msselection_p = new MSSelection();
91 0 : spw_selection = false;
92 :
93 0 : agents_p = NULL;
94 0 : agentCount_p=0;
95 0 : opts_p = NULL;
96 :
97 0 : logSink_p = LogSink(LogMessage::NORMAL, false);
98 :
99 0 : nant = 0;
100 0 : nifr = 0;
101 0 : nfeed = 0;
102 0 : nfeedcorr = 0;
103 :
104 0 : setdata_p = false;
105 0 : selectdata_p = false;
106 : // setupAgentDefaults();
107 :
108 0 : quack_agent_exists = false;
109 0 : }
110 :
111 : // -----------------------------------------------------------------------
112 : // Constructor
113 : // constructs and attaches to MS
114 : // -----------------------------------------------------------------------
115 0 : Flagger::Flagger ( MeasurementSet &mset ) : mssel_p(0), vs_p(0), scan_looping(false)
116 : {
117 0 : msselection_p = new MSSelection();
118 0 : spw_selection = false;
119 :
120 0 : agents_p = NULL;
121 0 : agentCount_p=0;
122 0 : opts_p = NULL;
123 :
124 0 : logSink_p = LogSink(LogMessage::NORMAL, false);
125 :
126 0 : nant = 0;
127 0 : setdata_p = false;
128 0 : selectdata_p = false;
129 0 : attach(mset);
130 :
131 0 : quack_agent_exists = false;
132 0 : }
133 :
134 0 : Flagger::~Flagger ()
135 : {
136 : /*
137 : jmlarsen: Eh? The following code is probably a bug
138 : (it closes the MS when the length() is 0)
139 : What about calling detach()?
140 : */
141 :
142 0 : if ( !ms.tableName().length()) {
143 0 : os << "Flagger closing out " << ms.tableName() << LogIO::POST;
144 0 : ms.flush();
145 0 : ms.relinquishAutoLocks(true);
146 0 : ms.unlock();
147 :
148 : //originalms = NULL;
149 0 : originalms_p = NULL;
150 : }
151 :
152 0 : if (vs_p) delete vs_p;
153 0 : vs_p = 0;
154 :
155 : if (dbg) cout << "Flagger destructor :: about to clean mssel_p" << endl;
156 :
157 0 : if (mssel_p) delete mssel_p;
158 0 : mssel_p = 0;
159 :
160 : if (dbg) cout << "Flagger destructor :: cleaned mssel_p" << endl;
161 :
162 0 : if (msselection_p) delete msselection_p;
163 0 : msselection_p = 0;
164 :
165 0 : if (agents_p) delete agents_p;
166 0 : agents_p = NULL;
167 :
168 0 : if (opts_p) delete opts_p;
169 0 : opts_p = NULL;
170 0 : }
171 :
172 : // -----------------------------------------------------------------------
173 : // queryOptions
174 : // Returns record of available options and their default values
175 : // -----------------------------------------------------------------------
176 0 : const RecordInterface & Flagger::defaultOptions ()
177 : {
178 0 : static Record rec;
179 : // create record description on first entry
180 0 : if ( !rec.nfields() )
181 : {
182 0 : rec.defineRecord(RF_GLOBAL, Record());
183 0 : rec.define(RF_TRIAL, false);
184 0 : rec.define(RF_RESET, false);
185 :
186 0 : rec.setComment(RF_GLOBAL, "Record of global parameters applied to all methods");
187 0 : rec.setComment(RF_TRIAL, "T for trial run (no flags written out)");
188 0 : rec.setComment(RF_RESET, "T to reset existing flags before running");
189 : }
190 0 : return rec;
191 : }
192 :
193 : // -----------------------------------------------------------------------
194 : // Flagger::attach
195 : // attaches to MS
196 : // -----------------------------------------------------------------------
197 :
198 0 : bool Flagger::attach( MeasurementSet &mset, Bool setAgentDefaults )
199 : {
200 0 : nant = 0;
201 0 : setdata_p = false;
202 0 : selectdata_p = false;
203 0 : if (vs_p)
204 0 : delete vs_p;
205 0 : vs_p = 0;
206 0 : if (mssel_p)
207 0 : delete mssel_p;
208 0 : mssel_p = 0;
209 0 : if (setAgentDefaults)
210 0 : setupAgentDefaults();
211 0 : ms = mset;
212 :
213 0 : originalms = mset;
214 0 : originalms_p = &mset;
215 :
216 : // Matrix<Int> noselection;
217 : // Use default sort order - and no scratch cols....
218 : // vs_p = new VisSet(*originalms_p,noselection,0.0);
219 : // vs_p->resetVisIter(sort, 0.0);
220 :
221 0 : Matrix<Int> noselection;
222 0 : Block<Int> bsort(3);
223 0 : bsort[0] = MS::FIELD_ID;
224 0 : bsort[1] = MS::DATA_DESC_ID;
225 0 : bsort[2] = MS::TIME;
226 :
227 : // extract various interesting info from the MS
228 : // obtain number of distinct time slots
229 0 : MSColumns msc(ms);
230 0 : Vector<Double> time( msc.time().getColumn() );
231 :
232 0 : uInt nrows = time.nelements();
233 : unsigned ntime;
234 :
235 0 : if (nrows == 0) ntime = 0;
236 : else {
237 : Bool dum;
238 0 : Double *tp = time.getStorage(dum);
239 :
240 0 : std::sort(tp, tp + nrows);
241 :
242 0 : ntime = 1;
243 0 : for (unsigned i = 0; i < nrows-1; i++) {
244 0 : if (tp[i] < tp[i+1]) ntime += 1;
245 : }
246 : }
247 :
248 : // obtain central frequencies of spws.
249 0 : const MSSpectralWindow spwin( ms.spectralWindow() );
250 0 : ScalarColumn<Double> sfreqs(spwin, "REF_FREQUENCY");
251 0 : spwfreqs.resize();
252 0 : spwfreqs = sfreqs.getColumn();
253 0 : spwfreqs *= 1e+6;
254 :
255 : // obtain number of antennas and interferometers
256 0 : const MSAntenna msant( ms.antenna() );
257 0 : nant = msant.nrow();
258 0 : nifr = nant*(nant+1)/2; // cheap & dirty
259 0 : ScalarColumn<String> names(msant,"NAME");
260 0 : antnames.resize();
261 0 : antnames = names.getColumn();
262 0 : antnames.apply(stringUpper);
263 : // cerr<<"Antenna names: "<<antnames<<endl;
264 : // map ifrs to antennas
265 0 : ifr2ant1.resize(nifr);
266 0 : ifr2ant1.set(-1);
267 0 : ifr2ant2.resize(nifr);
268 0 : ifr2ant2.set(-1);
269 0 : for( uInt i1=0; i1<nant; i1++ ) {
270 0 : for( uInt i2=0; i2<=i1; i2++ ) {
271 0 : uInt ifr = ifrNumber(i1,i2);
272 0 : ifr2ant1(ifr) = i1;
273 0 : ifr2ant2(ifr) = i2;
274 : }
275 : }
276 :
277 : // get feed info
278 0 : const MSFeed msfeed( ms.feed() );
279 0 : nfeed = msfeed.nrow();
280 0 : nfeedcorr = nfeed*(nfeed+1)/2;
281 :
282 0 : sprintf(str,"attached MS %s: %d rows, %d times, %d baselines\n",ms.tableName().chars(),nrows,ntime,nifr);
283 : //os << "--------------------------------------------------" << LogIO::POST;
284 0 : os<<str<<LogIO::POST;
285 :
286 0 : return true;
287 0 : }
288 :
289 : // -----------------------------------------------------------------------
290 : // Flagger::detach
291 : // detaches from MS
292 : // -----------------------------------------------------------------------
293 0 : void Flagger::detach()
294 : {
295 0 : if ( ms.tableName().length() == 0) {
296 0 : os << "no measurement set was attached" << LogIO::POST;
297 : }
298 : else {
299 0 : os << "detaching from MS " << ms.tableName() << LogIO::POST;
300 0 : nant = 0;
301 0 : setdata_p = false;
302 0 : selectdata_p = false;
303 0 : ms.flush();
304 0 : ms.relinquishAutoLocks(true);
305 0 : ms.unlock();
306 0 : ms = MeasurementSet();
307 : }
308 0 : }
309 :
310 : /************************************ DATA SELECTION **************************************/
311 : /* return: true iff succesful
312 : */
313 0 : Bool Flagger::selectdata(Bool useoriginalms,
314 : String field, String spw, String array,
315 : String feed, String scan,
316 : String baseline, String uvrange, String time,
317 : String correlation, String intent, String observation)
318 : {
319 : if (dbg) cout << "selectdata: "
320 : << "useoriginalms=" << useoriginalms
321 : << " field=" << field << " spw=" << spw
322 : << " array=" << array << " feed=" << feed
323 : << " scan=" << scan << " baseline=" << baseline
324 : << " uvrange=" << uvrange << " time=" << time
325 : << " correlation=" << correlation << " scan intent="<< intent
326 : << " observation=" << observation << endl;
327 :
328 0 : LogIO os(LogOrigin("Flagger", "selectdata()", WHERE));
329 0 : if (ms.isNull()) {
330 : os << LogIO::SEVERE << "NO MeasurementSet attached"
331 0 : << LogIO::POST;
332 0 : return false;
333 : }
334 :
335 : /* If a diff sort-order is required, put it in here, and
336 : create the MSs with the new ordering */
337 :
338 0 : MeasurementSet *tms = NULL;
339 0 : if ( useoriginalms ) tms = &originalms;
340 : else
341 : {
342 0 : if (!mssel_p)
343 : {
344 0 : cout << "Flagger::selectdata -> mssel_p is NULL !!!" << endl;
345 0 : return false;
346 : }
347 0 : tms = mssel_p;
348 : }
349 :
350 : if (dbg) cout << "Setting selection strings" << endl;
351 :
352 :
353 : /* Row-Selection */
354 :
355 0 : if (spw.length() == 0 && uvrange.length() > 0) {
356 0 : spw = String("*");
357 : }
358 :
359 0 : const String dummyExpr = String("");
360 0 : if (msselection_p) {
361 0 : delete msselection_p;
362 0 : msselection_p = NULL;
363 0 : spw_selection = false;
364 :
365 : }
366 0 : msselection_p = new MSSelection(*tms,
367 0 : MSSelection::PARSE_NOW,
368 0 : (const String)time,
369 0 : (const String)baseline,
370 0 : (const String)field,
371 0 : (const String)spw,
372 0 : (const String)uvrange,
373 : dummyExpr, // taqlExpr
374 : dummyExpr, // corrExpr
375 0 : (const String)scan,
376 0 : (const String)array,
377 0 : (const String)intent,
378 0 : (const String)observation);
379 0 : spw_selection = ((spw != "" && spw != "*") || uvrange.length() > 0);
380 : // multiple spw agents are also needed for uvrange selections...
381 :
382 0 : selectdata_p = false;
383 : /* Print out selection info - before selecting ! */
384 : if (dbg)
385 : {
386 : cout.precision(16);
387 : cout << "Antenna 1 : " << msselection_p->getAntenna1List() << endl;
388 : cout << "Antenna 2 : " << msselection_p->getAntenna2List() << endl;
389 : Matrix<Int> baselinelist(msselection_p->getBaselineList());
390 : IPosition shp = baselinelist.shape();
391 : if (shp.product() < 20)
392 : {
393 : IPosition transposed = shp;
394 : transposed[0]=shp[1]; transposed[1]=shp[0];
395 : Matrix<Int> blist(transposed);
396 : for (Int i=0;i<shp[0];i++)
397 : for (Int j=0;j<shp[1];j++)
398 : blist(j,i) = baselinelist(i,j);
399 : cout << "Baselines : " << blist << endl;
400 : }
401 : cout << "Fields : " << msselection_p->getFieldList() << endl;
402 : cout << "Spw : " << msselection_p->getSpwList() << endl;
403 : cout << "SpwChans : " << msselection_p->getChanList() << endl;
404 : Matrix<Double> tlist(msselection_p->getTimeList());
405 : cout << "Time : " << tlist << endl;
406 : cout << "Time : " << endl;
407 : for(Int i=0;i<(tlist.shape())[1];i++)
408 : cout << "[" << MVTime(tlist(0,i)/C::day).string(MVTime::DMY,7) << " ~ " << MVTime(tlist(1,i)/C::day).string(MVTime::DMY,7) << "]" << endl;
409 : cout << "UVrange : " << msselection_p->getUVList() << endl;
410 : cout << "UV Units : " << msselection_p->getUVUnitsList() << endl;
411 : cout << "Scans : " << msselection_p->getScanList() << endl;
412 : cout << "Arrays : " << msselection_p->getSubArrayList() << endl;
413 : // cout << "Feeds : " << msselection_p->getFeedList() << endl;
414 : cout << "Scan intent : " << msselection_p->getStateObsModeList() << endl;
415 : cout << "Observation : " << msselection_p->getObservationList() << endl;
416 :
417 : }
418 :
419 : /* Correlations */
420 0 : correlations_p.resize(0);
421 0 : string tcorr[50];
422 0 : Regex delim("(,| )+");
423 0 : Int ncorr = split(correlation, tcorr, 50, delim);
424 0 : correlations_p.resize(ncorr);
425 0 : for(Int i=0;i<ncorr;i++) correlations_p[i] = upcase(String(tcorr[i]));
426 :
427 : if (dbg) cout << "Correlations : " << correlations_p << endl;
428 :
429 0 : selectdata_p = true;
430 0 : return true;
431 0 : }
432 :
433 0 : Bool Flagger::setdata(
434 : String field, String spw, String array,
435 : String feed, String scan,
436 : String baseline, String uvrange, String time,
437 : String correlation, String intent, String observation)
438 : {
439 : if (dbg) cout << "setdata: "
440 : << " field=" << field << " spw=" << spw
441 : << " array=" << array << " feed=" << feed
442 : << " scan=" << scan << " baseline=" << baseline
443 : << " uvrange=" << uvrange << " time=" << time
444 : << " correlation=" << correlation << " scan intent=" << intent
445 : << " observation=" << observation << endl;
446 :
447 0 : LogIO os(LogOrigin("Flagger", "setdata()", WHERE));
448 0 : setdata_p = false;
449 :
450 :
451 : /* check the MS */
452 0 : if (ms.isNull()) {
453 : os << LogIO::SEVERE << "NO MeasurementSet attached"
454 0 : << LogIO::POST;
455 0 : return false;
456 : }
457 :
458 : /* Parse selection parameters */
459 0 : if (!spw.length()) spw = String("*");
460 :
461 0 : if (!selectdata(true,field,spw,array,feed,scan, baseline,uvrange,time,correlation,intent,observation))
462 : {
463 : os << LogIO::SEVERE << "Selection failed !!"
464 0 : << LogIO::POST;
465 0 : return false;
466 : }
467 :
468 : /* Create selected reference MS */
469 0 : MeasurementSet mssel_p2(*originalms_p);
470 :
471 0 : msselection_p->getSelectedMS(mssel_p2, String(""));
472 :
473 : if (dbg) cout << "Original ms has nrows : " << originalms.nrow() << LogIO::POST;
474 : if (dbg) cout << "Selected ms has " << mssel_p2.nrow() << " rows." << LogIO::POST;
475 :
476 0 : if ( mssel_p2.nrow() ) {
477 0 : if (mssel_p) {
478 0 : delete mssel_p;
479 0 : mssel_p=NULL;
480 : }
481 0 : if (mssel_p2.nrow() == originalms_p->nrow()) {
482 0 : mssel_p = new MeasurementSet(*originalms_p);
483 : }
484 : else {
485 0 : mssel_p = new MeasurementSet(mssel_p2);
486 : }
487 : if (dbg)cout << "assigned new MS to mssel_p" << endl;
488 0 : ScalarColumn<String> fname( mssel_p->field(),"NAME" );
489 : if (dbg)cout << "fields : " << fname.getColumn() << endl;
490 :
491 : //mssel_p->rename("selectedms",Table::New);
492 : //mssel_p->flush();
493 0 : }
494 : else {
495 0 : os << LogIO::WARN << "Selected MS has zero rows" << LogIO::POST;
496 : //mssel_p = &originalms;
497 0 : if (mssel_p) {
498 0 : delete mssel_p;
499 0 : mssel_p=NULL;
500 : }
501 0 : return false;
502 : }
503 :
504 : /* Print out selection info - before selecting ! */
505 0 : if (mssel_p->nrow()!=ms.nrow()) {
506 : os << "By selection " << originalms.nrow() << " rows are reduced to "
507 0 : << mssel_p->nrow() << LogIO::POST;
508 : }
509 : else {
510 0 : os << "Selection did not drop any rows" << LogIO::NORMAL3;
511 : }
512 : /* Create a vis iter */
513 0 : Matrix<Int> noselection;
514 0 : Block<int> sort2(4);
515 : //sort2[0] = MS::SCAN_NUMBER;
516 : // Do scan priority only if quacking
517 0 : sort2[0]= MS::ARRAY_ID;
518 0 : sort2[1]= MS::FIELD_ID;
519 0 : sort2[2]= MS::DATA_DESC_ID;
520 0 : sort2[3] = MS::TIME;
521 : /* Set the chunk time interval to be the full time-range.
522 : - FlagExaminer counters depend on this.
523 : However, autoflag agents like tfcrop will take a performance hit
524 : since it will store flags for the entire chunk in memory
525 : */
526 0 : Double timeInterval = 7.0e9; //a few thousand years
527 : //Double timeInterval = 6000; // seconds : 100 minutes.
528 :
529 0 : if (vs_p) {
530 0 : delete vs_p; vs_p = NULL;
531 : }
532 0 : if (!vs_p) {
533 0 : if (!mssel_p)
534 0 : throw AipsError("No measurement set selection available");
535 : //vs_p = new VisSet(*mssel_p,sort,noselection,0.0);
536 :
537 0 : Bool addScratch = false;
538 0 : vs_p = new VisSet(*mssel_p, sort2, noselection,
539 0 : addScratch, timeInterval);
540 : // Use default sort order - and no scratch cols....
541 : //vs_p = new VisSet(*mssel_p,noselection,0.0);
542 : }
543 0 : vs_p->resetVisIter(sort2, timeInterval);
544 :
545 : // Optimize for iterating randomly through one MS
546 0 : vs_p->iter().slurp();
547 :
548 : /* Channel selection */ // Always select all chans
549 0 : selectDataChannel();
550 0 : ms = *mssel_p;
551 :
552 :
553 0 : setdata_p = true;
554 0 : return true;
555 0 : }
556 :
557 :
558 0 : Bool Flagger::selectDataChannel(){
559 :
560 0 : if (!vs_p || !msselection_p) return false;
561 : /* Set channel selection in visiter */
562 : /* Set channel selection per spectral window - from msselection_p->getChanList(); */
563 : // this is needed when "setdata" is used to select data for autoflag algorithms
564 : // This should not be done for manual flagging, because the channel indices for the
565 : // selected subset start from zero - and throw everything into confusion.
566 0 : Vector<Int> spwlist = msselection_p->getSpwList();
567 :
568 0 : if ( spwlist.nelements() && spw_selection ) {
569 0 : Matrix<Int> spwchan = msselection_p->getChanList();
570 0 : IPosition cshp = spwchan.shape();
571 0 : if ( (Int)spwlist.nelements() > (spwchan.shape())[0] )
572 0 : cout << "WARN : Using only the first channel range per spw" << endl;
573 0 : for( uInt i=0;i<spwlist.nelements();i++ )
574 : {
575 0 : Int j=0;
576 0 : for ( j=0;j<(spwchan.shape())[0];j++ )
577 0 : if ( spwchan(j,0) == spwlist[i] ) break;
578 0 : vs_p->iter().selectChannel(1, Int(spwchan(j,1)),
579 0 : Int(spwchan(j,2)-spwchan(j,1)+1),
580 0 : Int(spwchan(j,3)), spwlist[i]);
581 : }
582 0 : }
583 : else{
584 0 : return false;
585 : }
586 :
587 :
588 0 : return true;
589 0 : }
590 :
591 :
592 :
593 : /************************** Set Manual Flags ***************************/
594 :
595 0 : Bool Flagger::fillSelections(Record &rec)
596 : {
597 :
598 : //TableExprNode ten = msselection_p->toTableExprNode(originalms_p);
599 :
600 : /* Fill in record for selected values. */
601 : /* Field ID */
602 0 : Vector<Int> fieldlist = msselection_p->getFieldList();
603 0 : if (fieldlist.nelements())
604 : {
605 0 : RecordDesc flagDesc;
606 0 : flagDesc.addField(RF_FIELD, TpArrayInt);
607 0 : Record flagRec(flagDesc);
608 0 : flagRec.define(RF_FIELD, fieldlist);
609 0 : rec.mergeField(flagRec, RF_FIELD, RecordInterface::OverwriteDuplicates);
610 0 : }
611 : /* BASELINE */
612 0 : Matrix<Int> baselinelist = msselection_p->getBaselineList();
613 :
614 : /* Here, it may be necessary to convert negative indices to positive ones
615 : (see CAS-2021).
616 :
617 : For now, fail cleanly if getBaselineList returned negative indices.
618 : */
619 :
620 0 : if (baselinelist.nelements())
621 : {
622 0 : IPosition shp = baselinelist.shape();
623 : if (dbg)cout << "Original shape of baselinelist : " << shp << endl;
624 0 : IPosition transposed = shp;
625 0 : transposed[0] = shp[1]; transposed[1] = shp[0];
626 0 : Matrix<Int> blist(transposed);
627 :
628 0 : for(Int i=0; i < shp[0]; i++)
629 0 : for(Int j=0; j < shp[1]; j++) {
630 :
631 0 : if (baselinelist(i, j) < 0) {
632 0 : throw AipsError("Sorry, negated antenna selection (such as '!2') is not supported by flagger (CAS-2021)");
633 : }
634 0 : blist(j, i) = baselinelist(i, j);
635 : }
636 :
637 0 : RecordDesc flagDesc;
638 0 : flagDesc.addField(RF_BASELINE, TpArrayInt);
639 0 : Record flagRec(flagDesc);
640 0 : flagRec.define(RF_BASELINE, blist);
641 0 : rec.mergeField(flagRec, RF_BASELINE, RecordInterface::OverwriteDuplicates);
642 0 : }
643 : /* FEED */
644 : /*
645 : Matrix<Int> feedlist = msselection_p->getFeedList();
646 : if (feedlist.nelements())
647 : {
648 : IPosition shp = feedlist.shape();
649 : if (dbg)cout << "Original shape of feedlist : " << shp << endl;
650 : IPosition transposed = shp;
651 : transposed[0]=shp[1]; transposed[1]=shp[0];
652 : Matrix<Int> blist(transposed);
653 : for(Int i=0;i<shp[0];i++)
654 : for(Int j=0;j<shp[1];j++)
655 : blist(j,i) = feedlist(i,j);
656 : // need to add 1 because RFASelector expects 1-based indices.
657 :
658 : RecordDesc flagDesc;
659 : flagDesc.addField(RF_FEED, TpArrayInt);
660 : Record flagRec(flagDesc);
661 : flagRec.define(RF_FEED, blist);
662 : rec.mergeField(flagRec, RF_FEED, RecordInterface::OverwriteDuplicates);
663 : }
664 : */
665 : /* TIME */
666 0 : Matrix<Double> timelist = msselection_p->getTimeList();
667 0 : if (timelist.nelements())
668 : {
669 : /* Times need to be in MJD */
670 0 : for( Int i=0;i<(timelist.shape())[0];i++ )
671 0 : for( Int j=0;j<(timelist.shape())[1];j++ )
672 0 : timelist(i,j) /= (Double)(24*3600);
673 : //for( Int i=0;i<(timelist.shape())[0];i++ )
674 : // for( Int j=0;j<(timelist.shape())[1];j++ )
675 : // cout << "timelist(" << i << ", " << j << ")="
676 : // << timelist(i,j) << endl;
677 0 : RecordDesc flagDesc;
678 0 : flagDesc.addField(RF_TIMERANGE, TpArrayDouble);
679 0 : Record flagRec(flagDesc);
680 0 : flagRec.define(RF_TIMERANGE, timelist);
681 0 : rec.mergeField(flagRec, RF_TIMERANGE, RecordInterface::OverwriteDuplicates);
682 0 : }
683 : /* RF_CORR */
684 0 : if (correlations_p.nelements())
685 : {
686 0 : RecordDesc flagDesc;
687 0 : flagDesc.addField(RF_CORR, TpArrayString);
688 0 : Record flagRec(flagDesc);
689 0 : flagRec.define(RF_CORR, correlations_p);
690 0 : rec.mergeField(flagRec, RF_CORR, RecordInterface::OverwriteDuplicates);
691 0 : }
692 : /* Array ID */
693 0 : Vector<Int> arraylist = msselection_p->getSubArrayList();
694 0 : if (arraylist.nelements())
695 : {
696 0 : RecordDesc flagDesc;
697 0 : flagDesc.addField(RF_ARRAY, TpArrayInt);
698 0 : Record flagRec(flagDesc);
699 0 : flagRec.define(RF_ARRAY, arraylist);
700 0 : rec.mergeField(flagRec, RF_ARRAY, RecordInterface::OverwriteDuplicates);
701 0 : }
702 : /* Scan ID */
703 0 : Vector<Int> scanlist = msselection_p->getScanList();
704 0 : if (scanlist.nelements())
705 : {
706 0 : RecordDesc flagDesc;
707 0 : flagDesc.addField(RF_SCAN, TpArrayInt);
708 0 : Record flagRec(flagDesc);
709 0 : flagRec.define(RF_SCAN, scanlist);
710 0 : rec.mergeField(flagRec, RF_SCAN, RecordInterface::OverwriteDuplicates);
711 0 : }
712 : /* State ID (obsMode) */
713 0 : Vector<Int> stateobsmodelist = msselection_p->getStateObsModeList();
714 0 : if (stateobsmodelist.nelements())
715 : {
716 0 : RecordDesc flagDesc;
717 0 : flagDesc.addField(RF_INTENT, TpArrayInt);
718 0 : Record flagRec(flagDesc);
719 0 : flagRec.define(RF_INTENT, stateobsmodelist);
720 0 : rec.mergeField(flagRec, RF_INTENT, RecordInterface::OverwriteDuplicates);
721 0 : }
722 : /* Observation ID */
723 0 : Vector<Int> observationlist = msselection_p->getObservationList();
724 0 : if (observationlist.nelements())
725 : {
726 0 : RecordDesc flagDesc;
727 0 : flagDesc.addField(RF_OBSERVATION, TpArrayInt);
728 0 : Record flagRec(flagDesc);
729 0 : flagRec.define(RF_OBSERVATION, observationlist);
730 0 : rec.mergeField(flagRec, RF_OBSERVATION, RecordInterface::OverwriteDuplicates);
731 0 : }
732 :
733 0 : return true;
734 0 : }
735 :
736 :
737 :
738 : /*
739 : Sets up agents for mode = 'manualflag' and mode = 'summary'
740 : */
741 :
742 0 : Bool Flagger::setmanualflags(Bool autocorr,
743 : Bool unflag,
744 : String clipexpr,
745 : Vector<Double> cliprange,
746 : String clipcolumn,
747 : Bool outside,
748 : Bool channel_average,
749 : Double quackinterval,
750 : String quackmode,
751 : Bool quackincrement,
752 : String opmode,
753 : Double diameter,
754 : Double lowerlimit,
755 : Double upperlimit)
756 : {
757 : if (dbg)
758 : cout << "setmanualflags: "
759 : << "autocorr=" << autocorr
760 : << " unflag=" << unflag
761 : << " clipexpr=" << clipexpr
762 : << " cliprange=" << cliprange
763 : << " clipcolumn=" << clipcolumn
764 : << " outside=" << outside
765 : << " quackinterval=" << quackinterval
766 : << " quackmode=" << quackmode
767 : << " quackincrement=" << quackincrement
768 : << " opmode=" << opmode
769 : << " diameter=" << diameter
770 : << endl;
771 :
772 0 : LogIO os(LogOrigin("Flagger", "setmanualflags()", WHERE));
773 0 : if (ms.isNull()) {
774 : os << LogIO::SEVERE << "NO MeasurementSet attached"
775 0 : << LogIO::POST;
776 0 : return false;
777 : }
778 0 : if (!selectdata_p) {
779 : os << LogIO::SEVERE << "Please run selectdata with/without arguments before setmanualflags"
780 0 : << LogIO::POST;
781 0 : return false;
782 : }
783 :
784 : /* Fill in an agent record */
785 : /* This assumes that selectdata has already been called */
786 :
787 : /* Loop over SPW and chan ranges. */
788 :
789 0 : Vector<Int> spwlist = msselection_p->getSpwList();
790 :
791 0 : Matrix<Int> spwchan = msselection_p->getChanList();
792 0 : Matrix<Double> uvrangelist = msselection_p->getUVList();
793 0 : Vector<Bool> uvrangeunits = msselection_p->getUVUnitsList();
794 0 : IPosition cshp = spwchan.shape();
795 :
796 : /* If no selection depends on spw, ( i.e. no spw, chan, uvrange )
797 : then no need to make separate records for each spw. */
798 0 : bool separatespw = false;
799 : Int nrec;
800 0 : if (spwlist.nelements() && spw_selection) {
801 0 : separatespw = true;
802 0 : nrec = spwlist.nelements();
803 : }
804 : else {
805 0 : separatespw = false; nrec = 1;
806 : }
807 :
808 : if (dbg) {
809 : cout << "separatespw = " << separatespw << endl;
810 : cout << spwlist << endl;
811 : }
812 :
813 0 : scan_looping = false;
814 :
815 0 : for ( Int i=0; i < nrec; i++ ) {
816 0 : Record selrec;
817 0 : if (upcase(opmode).matches("FLAG") ||
818 0 : upcase(opmode).matches("SHADOW") ||
819 0 : upcase(opmode).matches("ELEVATION")) {
820 0 : selrec.define("id", String("select"));
821 :
822 0 : if (upcase(opmode).matches("SHADOW")) {
823 0 : selrec.define(RF_DIAMETER, Double(diameter));
824 : }
825 0 : else if (upcase(opmode).matches("ELEVATION")) {
826 0 : selrec.define(RF_LOWERLIMIT, Double(lowerlimit));
827 0 : selrec.define(RF_UPPERLIMIT, Double(upperlimit));
828 : }
829 : }
830 0 : else if (upcase(opmode).matches("SUMMARY")) {
831 0 : selrec.define("id", String("flagexaminer"));
832 : }
833 : else {
834 0 : throw AipsError("Unknown mode " + upcase(opmode));
835 : }
836 :
837 : /* Fill selections for all but spw, chan, corr */
838 0 : fillSelections(selrec);
839 :
840 0 : if (separatespw) {
841 : /* SPW ID */
842 : {
843 0 : RecordDesc flagDesc;
844 0 : flagDesc.addField(RF_SPWID, TpArrayInt);
845 0 : Record flagRec(flagDesc);
846 0 : Vector<Int> t_spw(1); t_spw[0] = spwlist[i];
847 0 : flagRec.define(RF_SPWID, t_spw);
848 0 : selrec.mergeField(flagRec, RF_SPWID, RecordInterface::OverwriteDuplicates);
849 0 : }
850 :
851 : /* reform chan ranges */
852 0 : Int ccount=0;
853 0 : for( Int j=0;j<cshp[0];j++ ) if ( spwlist[i] == spwchan(j,0) ) ccount++;
854 0 : Matrix<Int> chanlist(2,ccount); chanlist.set(0);
855 :
856 0 : ccount=0;
857 0 : for( Int j=0;j<cshp[0];j++ )
858 : {
859 0 : if ( spwlist[i] == spwchan(j,0) )
860 : {
861 0 : chanlist(0,ccount) = spwchan(j,1);
862 0 : chanlist(1,ccount) = spwchan(j,2);
863 0 : if ( spwchan(j,3) > 1 )
864 0 : os << LogIO::WARN << ".... ignoring chan 'step' for manual flags" << LogIO::POST;
865 0 : ccount++;
866 : }
867 : }
868 : /* RF_CHANS */
869 : {
870 0 : RecordDesc flagDesc;
871 0 : flagDesc.addField(RF_CHANS, TpArrayInt);
872 0 : Record flagRec(flagDesc);
873 0 : flagRec.define(RF_CHANS, chanlist);
874 0 : selrec.mergeField(flagRec, RF_CHANS, RecordInterface::OverwriteDuplicates);
875 0 : }
876 :
877 : /* UV-RANGE */
878 0 : if (uvrangelist.nelements())
879 : {
880 0 : Matrix<Double> templist(uvrangelist.shape());
881 : /* Convert to Metres... */
882 : /* or complain if units are not metres... */
883 : // current spw : spwlist[i];
884 :
885 0 : if ( (templist.shape())[1] != (Int)uvrangeunits.nelements() )
886 0 : cout << "UVRANGE units are wrong length ! " << endl;
887 0 : for( Int j=0;j<(templist.shape())[1];j++ )
888 : {
889 0 : Double unit=1.0;
890 0 : if ( ! uvrangeunits[j] ) unit = C::c/(spwfreqs[spwlist[i]]/1e+6);
891 0 : for( Int k=0;k<(templist.shape())[0];k++ )
892 0 : templist(k,j) = uvrangelist(k,j) * unit ;
893 : }
894 :
895 0 : RecordDesc flagDesc;
896 0 : flagDesc.addField(RF_UVRANGE, TpArrayDouble);
897 0 : Record flagRec(flagDesc);
898 0 : flagRec.define(RF_UVRANGE, templist);
899 0 : selrec.mergeField(flagRec, RF_UVRANGE, RecordInterface::OverwriteDuplicates);
900 : if (dbg) cout << "uv list (m) : " << templist << endl;
901 0 : }
902 0 : }
903 :
904 : // Operation related parameters.
905 0 : if (upcase(opmode).matches("SHADOW") ) {
906 0 : RecordDesc flagDesc;
907 0 : flagDesc.addField(RF_SHADOW, TpBool);
908 0 : Record flagRec(flagDesc);
909 0 : flagRec.define(RF_SHADOW, true);
910 0 : selrec.mergeField(flagRec, RF_SHADOW, RecordInterface::OverwriteDuplicates);
911 0 : }
912 :
913 0 : else if (upcase(opmode).matches("ELEVATION") ) {
914 0 : RecordDesc flagDesc;
915 0 : flagDesc.addField(RF_ELEVATION, TpBool);
916 0 : Record flagRec(flagDesc);
917 0 : flagRec.define(RF_ELEVATION, true);
918 0 : selrec.mergeField(flagRec, RF_ELEVATION, RecordInterface::OverwriteDuplicates);
919 0 : }
920 :
921 : /* Flag Autocorrelations too? */
922 0 : if (autocorr) {
923 0 : RecordDesc flagDesc;
924 0 : flagDesc.addField(RF_AUTOCORR, TpBool);
925 0 : Record flagRec(flagDesc);
926 0 : flagRec.define(RF_AUTOCORR, autocorr);
927 0 : selrec.mergeField(flagRec, RF_AUTOCORR, RecordInterface::OverwriteDuplicates);
928 0 : }
929 :
930 : /* Unflag! */
931 0 : if (unflag) {
932 0 : RecordDesc flagDesc;
933 0 : flagDesc.addField(RF_UNFLAG, TpBool);
934 0 : Record flagRec(flagDesc);
935 0 : flagRec.define(RF_UNFLAG, unflag);
936 0 : selrec.mergeField(flagRec, RF_UNFLAG, RecordInterface::OverwriteDuplicates);
937 0 : }
938 :
939 : /* Reset flags before applying new ones */
940 : // ( I think... )
941 : /*
942 : {
943 : RecordDesc flagDesc;
944 : flagDesc.addField(RF_RESET, TpBool);
945 : Record flagRec(flagDesc);
946 : flagRec.define(RF_RESET, true);
947 : selrec.mergeField(flagRec, RF_RESET, RecordInterface::OverwriteDuplicates);
948 : }
949 : */
950 :
951 : /* Clip/FlagRange */
952 : /* Jira Casa 212 : Check if "clipexpr" has multiple
953 : comma-separated expressions
954 : and loop here, creating multiple clipRecs.
955 : The RFASelector will handle it. */
956 0 : if (clipexpr.length() && cliprange.nelements() == 2 &&
957 0 : (cliprange[0] < cliprange[1] ||
958 0 : (cliprange[0] <= cliprange[1] && !outside) // i.e. exact matching
959 : )
960 : )
961 : {
962 0 : RecordDesc flagDesc;
963 0 : if ( outside )
964 0 : flagDesc.addField(RF_CLIP, TpRecord);
965 : else
966 0 : flagDesc.addField(RF_FLAGRANGE, TpRecord);
967 :
968 0 : Record flagRec(flagDesc);
969 :
970 0 : RecordDesc clipDesc;
971 0 : clipDesc.addField(RF_EXPR, TpString);
972 0 : clipDesc.addField(RF_MIN, TpDouble);
973 0 : clipDesc.addField(RF_MAX, TpDouble);
974 0 : clipDesc.addField(RF_CHANAVG, TpBool);
975 0 : Record clipRec(clipDesc);
976 0 : clipRec.define(RF_EXPR, clipexpr);
977 0 : clipRec.define(RF_MIN, cliprange[0]);
978 0 : clipRec.define(RF_MAX, cliprange[1]);
979 0 : clipRec.define(RF_CHANAVG, channel_average);
980 :
981 0 : if ( outside )
982 : {
983 0 : flagRec.defineRecord(RF_CLIP, clipRec);
984 0 : selrec.mergeField(flagRec, RF_CLIP, RecordInterface::OverwriteDuplicates);
985 : }
986 : else
987 : {
988 0 : flagRec.defineRecord(RF_FLAGRANGE, clipRec);
989 0 : selrec.mergeField(flagRec, RF_FLAGRANGE, RecordInterface::OverwriteDuplicates);
990 : }
991 :
992 : /* clip column */
993 0 : if (!clipcolumn.length()) clipcolumn=String("DATA");
994 0 : RecordDesc flagDesc2;
995 0 : flagDesc2.addField(RF_COLUMN, TpString);
996 0 : Record flagRec2(flagDesc2);
997 0 : flagRec2.define(RF_COLUMN, clipcolumn);
998 0 : selrec.mergeField(flagRec2, RF_COLUMN, RecordInterface::OverwriteDuplicates);
999 :
1000 0 : }
1001 :
1002 : /* Quack! */
1003 0 : if (quackinterval > 0.0) {
1004 : //Reset the Visiter to have SCAN before time
1005 0 : scan_looping = true;
1006 0 : Block<int> sort2(5);
1007 0 : sort2[0] = MS::ARRAY_ID;
1008 0 : sort2[1] = MS::FIELD_ID;
1009 0 : sort2[2] = MS::DATA_DESC_ID;
1010 0 : sort2[3] = MS::SCAN_NUMBER;
1011 0 : sort2[4] = MS::TIME;
1012 0 : Double timeInterval = 7.0e9; //a few thousand years
1013 :
1014 0 : vs_p->resetVisIter(sort2, timeInterval);
1015 0 : vs_p->iter().slurp();
1016 : //lets make sure the data channel selection is done
1017 0 : selectDataChannel();
1018 :
1019 0 : RecordDesc flagDesc;
1020 0 : flagDesc.addField(RF_QUACK, TpArrayDouble);
1021 0 : flagDesc.addField(RF_QUACKMODE, TpString);
1022 0 : flagDesc.addField(RF_QUACKINC, TpBool);
1023 :
1024 0 : Vector<Double> quackparams(2);
1025 0 : quackparams[0] = quackinterval;
1026 0 : quackparams[1] = 0.0;
1027 :
1028 0 : Record flagRec(flagDesc);
1029 0 : flagRec.define(RF_QUACK, quackparams);
1030 0 : flagRec.define(RF_QUACKMODE, quackmode);
1031 0 : flagRec.define(RF_QUACKINC, quackincrement);
1032 0 : selrec.mergeField(flagRec, RF_QUACK , RecordInterface::OverwriteDuplicates);
1033 0 : selrec.mergeField(flagRec, RF_QUACKMODE, RecordInterface::OverwriteDuplicates);
1034 0 : selrec.mergeField(flagRec, RF_QUACKINC , RecordInterface::OverwriteDuplicates);
1035 :
1036 0 : quack_agent_exists = true;
1037 0 : }
1038 : /* end if opmode = ... */
1039 :
1040 : /* Add this agent to the list */
1041 0 : addAgent(selrec);
1042 0 : }
1043 :
1044 0 : return true;
1045 0 : }
1046 :
1047 0 : Bool Flagger::setautoflagparams(String algorithm,Record ¶meters)
1048 : {
1049 0 : LogIO os(LogOrigin("Flagger", "setautoflagparams()", WHERE));
1050 0 : if (ms.isNull()) {
1051 : os << LogIO::SEVERE << "NO MeasurementSet attached"
1052 0 : << LogIO::POST;
1053 0 : return false;
1054 : }
1055 0 : if (!selectdata_p) {
1056 : os << LogIO::SEVERE << "Please run setdata with/without arguments before setautoflagparams"
1057 0 : << LogIO::POST;
1058 0 : return false;
1059 : }
1060 :
1061 : /* Create an agent record */
1062 0 : Record selrec;
1063 0 : selrec = getautoflagparams(algorithm);
1064 0 : selrec.merge(parameters,Record::OverwriteDuplicates);
1065 :
1066 : /* special case for "sprej". need to parse a param.*/
1067 0 : if (algorithm.matches("sprej") && selrec.isDefined("fitspwchan")) {
1068 : /* Get the "fitspwchan" string" */
1069 : /* Pass this through the msselection parser and getChanList */
1070 : /* Construct a list of regions records from this list */
1071 0 : String fitspwchan;
1072 0 : selrec.get(RecordFieldId("fitspwchan"),fitspwchan);
1073 :
1074 0 : if (fitspwchan.length())
1075 : {
1076 : /* Parse it */
1077 0 : const String dummy("");
1078 0 : MSSelection tmpmss(*mssel_p,MSSelection::PARSE_NOW,
1079 : dummy,dummy, dummy, fitspwchan, dummy,
1080 0 : dummy, dummy, dummy, dummy);
1081 0 : Matrix<Int> spwchanlist = tmpmss.getChanList();
1082 0 : Vector<Int> spwlist = tmpmss.getSpwList();
1083 :
1084 : /* Create region record template */
1085 0 : RecordDesc regdesc;
1086 0 : for(uInt i=0;i<spwlist.nelements();i++)
1087 0 : regdesc.addField(String(i),TpRecord);
1088 0 : Record regions(regdesc);
1089 :
1090 : /* reform chan ranges */
1091 0 : Int ccount=0;
1092 0 : IPosition cshp = spwchanlist.shape();
1093 0 : for(uInt i=0;i<spwlist.nelements();i++)
1094 : {
1095 0 : ccount=0;
1096 0 : for( Int j=0;j<cshp[0];j++ )
1097 0 : if ( spwlist[i] == spwchanlist(j,0) ) ccount++;
1098 0 : Matrix<Int> chanlist(2,ccount); chanlist.set(0);
1099 :
1100 0 : ccount=0;
1101 0 : for( Int j=0;j<cshp[0];j++ )
1102 : {
1103 0 : if ( spwlist[i] == spwchanlist(j,0) )
1104 : {
1105 0 : chanlist(0,ccount) = spwchanlist(j,1);
1106 0 : chanlist(1,ccount) = spwchanlist(j,2);
1107 0 : if ( spwchanlist(j,3) > 1 )
1108 0 : os << LogIO::WARN << ".... ignoring chan 'step' for 'sprej' fitting" << LogIO::POST;
1109 0 : ccount++;
1110 : }
1111 : }
1112 :
1113 0 : RecordDesc spwDesc;
1114 0 : spwDesc.addField(RF_SPWID, TpInt);
1115 0 : spwDesc.addField(RF_CHANS, TpArrayInt);
1116 0 : Record spwRec(spwDesc);
1117 0 : spwRec.define(RF_SPWID, spwlist[i]);
1118 0 : spwRec.define(RF_CHANS, chanlist);
1119 :
1120 : /* create a single region record */
1121 : /* add this to the list of region records */
1122 0 : regions.defineRecord(String(i),spwRec);
1123 0 : }
1124 :
1125 : /* Attach the list of region records */
1126 0 : selrec.defineRecord(RF_REGION, regions);
1127 0 : }
1128 0 : }
1129 :
1130 : /* Add this agent to the list */
1131 0 : addAgent(selrec);
1132 :
1133 0 : return true;
1134 0 : }
1135 :
1136 :
1137 0 : Record Flagger::getautoflagparams(String algorithm)
1138 : {
1139 0 : LogIO os(LogOrigin("Flagger", "getautoflagparams()", WHERE));
1140 :
1141 : // Use "RFATimeMedian::getDefaults()" !!!!!!!
1142 :
1143 0 : Record defrecord;
1144 0 : if ( agent_defaults.isDefined(algorithm) )
1145 : {
1146 0 : RecordFieldId rid(algorithm);
1147 0 : defrecord = agent_defaults.asRecord(rid);
1148 0 : defrecord.define("id",algorithm);
1149 :
1150 0 : if (defrecord.isDefined("expr")) defrecord.define("expr","ABS I");
1151 0 : }
1152 0 : return defrecord;
1153 :
1154 0 : }
1155 :
1156 : /* Clean up all selections */
1157 : // TODO Later add in the ability to clear specified
1158 : // agent types.
1159 : //
1160 : // subRecord = agents_p->getField(RecordFieldId(x));
1161 : // if ( subRecord.isDefined('id') && 'id' is "select" )
1162 : // agents_p->removeField(RecordFieldId(x));
1163 : //
1164 0 : Bool Flagger::clearflagselections(Int recordindex)
1165 : {
1166 0 : LogIO os(LogOrigin("Flagger", "clearflagselections()", WHERE));
1167 :
1168 0 : if ( agents_p && agents_p->nfields() )
1169 : {
1170 0 : if ( recordindex >= 0 )
1171 : {
1172 : if (dbg) cout << "Deleting only agent : " << recordindex << endl;
1173 0 : agents_p->removeField(RecordFieldId(recordindex));
1174 : }
1175 : else
1176 : {
1177 : if (dbg) cout << "Deleting all agents" << endl;
1178 0 : delete agents_p;
1179 0 : agents_p =0;
1180 0 : agentCount_p = 0;
1181 : }
1182 : }
1183 :
1184 : // printflagselections();
1185 :
1186 0 : return true;
1187 0 : }
1188 :
1189 0 : Bool Flagger::printflagselections()
1190 : {
1191 0 : LogIO os(LogOrigin("Flagger", "printflagselections()", WHERE));
1192 0 : if ( agents_p )
1193 : {
1194 0 : os << "Current list of agents : " << agents_p << LogIO::POST;
1195 0 : ostringstream out;
1196 0 : agents_p->print(out);
1197 0 : os << out.str() << LogIO::POST;
1198 0 : }
1199 0 : else os << " No current agents " << LogIO::POST;
1200 :
1201 0 : return true;
1202 0 : }
1203 :
1204 0 : Bool Flagger::addAgent(RecordInterface &newAgent)
1205 : {
1206 0 : if (!agents_p)
1207 : {
1208 0 : agentCount_p = 0;
1209 0 : agents_p = new Record;
1210 : if (dbg) cout << "creating new agent" << endl;
1211 : }
1212 :
1213 0 : ostringstream fieldName;
1214 0 : fieldName << agentCount_p++;
1215 0 : agents_p->defineRecord(fieldName.str(), newAgent);
1216 :
1217 : // printflagselections();
1218 :
1219 0 : return true;
1220 0 : }
1221 :
1222 :
1223 : // computes IFR index, given two antennas
1224 0 : uInt Flagger::ifrNumber ( Int ant1,Int ant2 ) const
1225 : {
1226 0 : if ( ant1<ant2 )
1227 0 : return ifrNumber(ant2,ant1);
1228 0 : return ant1*(ant1+1)/2 + ant2;
1229 : }
1230 : //TODO
1231 : // Here, fill in correct indices for the baseline ordering in the MS
1232 : // All agents will see this and will be fine.
1233 :
1234 : // computes vector of IFR indeces, given two antennas
1235 0 : Vector<Int> Flagger::ifrNumbers ( Vector<Int> ant1,Vector<Int> ant2 ) const
1236 : {
1237 0 : unsigned n = ant1.nelements();
1238 0 : Vector<Int> r(n);
1239 :
1240 0 : for (unsigned i = 0; i < n; i++) {
1241 0 : Int a1 = ant1(i);
1242 0 : Int a2 = ant2(i);
1243 0 : if (a1 > a2) {
1244 0 : r(i) = a1*(a1+1)/2 + a2;
1245 : }
1246 : else {
1247 0 : r(i) = a2*(a2+1)/2 + a1;
1248 : }
1249 : }
1250 0 : return r;
1251 0 : }
1252 :
1253 0 : void Flagger::ifrToAnt ( uInt &ant1,uInt &ant2,uInt ifr ) const
1254 : {
1255 0 : ant1 = ifr2ant1(ifr);
1256 0 : ant2 = ifr2ant2(ifr);
1257 0 : }
1258 :
1259 : // -----------------------------------------------------------------------
1260 : // Flagger::setupAgentDefaults
1261 : // Sets up record of available agents and their default parameters
1262 : // -----------------------------------------------------------------------
1263 0 : const RecordInterface & Flagger::setupAgentDefaults ()
1264 : {
1265 0 : agent_defaults = Record();
1266 0 : agent_defaults.defineRecord("timemed", RFATimeMedian::getDefaults());
1267 0 : agent_defaults.defineRecord("newtimemed", RFANewMedianClip::getDefaults());
1268 0 : agent_defaults.defineRecord("freqmed", RFAFreqMedian::getDefaults());
1269 0 : agent_defaults.defineRecord("sprej", RFASpectralRej::getDefaults());
1270 0 : agent_defaults.defineRecord("select", RFASelector::getDefaults());
1271 0 : agent_defaults.defineRecord("flagexaminer", RFAFlagExaminer::getDefaults());
1272 0 : agent_defaults.defineRecord("uvbin", RFAUVBinner::getDefaults());
1273 0 : agent_defaults.defineRecord("tfcrop", RFATimeFreqCrop::getDefaults());
1274 0 : return agent_defaults;
1275 : }
1276 :
1277 : // -----------------------------------------------------------------------
1278 : // Flagger::createAgent
1279 : // Creates flagging agent based on name
1280 : // -----------------------------------------------------------------------
1281 0 : std::shared_ptr<RFABase> Flagger::createAgent (const String &id,
1282 : RFChunkStats &chunk,
1283 : const RecordInterface &parms,
1284 : bool &only_selector)
1285 : {
1286 0 : if (id != "select" && id != "flagexaminer") {
1287 0 : only_selector = false;
1288 : }
1289 : // cerr << "Agent id: " << id << endl;
1290 0 : if ( id == "timemed" )
1291 0 : return std::shared_ptr<RFABase>(new RFATimeMedian(chunk, parms));
1292 0 : else if ( id == "newtimemed" )
1293 0 : return std::shared_ptr<RFABase>(new RFANewMedianClip(chunk, parms));
1294 0 : else if ( id == "freqmed" )
1295 0 : return std::shared_ptr<RFABase>(new RFAFreqMedian(chunk, parms));
1296 0 : else if ( id == "sprej" )
1297 0 : return std::shared_ptr<RFABase>(new RFASpectralRej(chunk, parms));
1298 0 : else if ( id == "select" )
1299 0 : return std::shared_ptr<RFABase>(new RFASelector(chunk, parms));
1300 0 : else if ( id == "flagexaminer" )
1301 0 : return std::shared_ptr<RFABase>(new RFAFlagExaminer(chunk, parms));
1302 0 : else if ( id == "uvbin" )
1303 0 : return std::shared_ptr<RFABase>(new RFAUVBinner(chunk, parms));
1304 0 : else if ( id == "tfcrop" )
1305 0 : return std::shared_ptr<RFABase>(new RFATimeFreqCrop(chunk, parms));
1306 : else
1307 0 : return std::shared_ptr<RFABase>();
1308 : }
1309 :
1310 :
1311 0 : void Flagger::summary( const RecordInterface &agents )
1312 : {
1313 : //os << "Autoflag summary will report results here" << LogIO::POST;
1314 0 : for(uInt i=0;i<agents.nfields(); i++){
1315 :
1316 0 : if (agents.dataType(i) != TpRecord){
1317 0 : os << "Unrecognized field: " << agents.name(i) << LogIO::EXCEPTION;
1318 : }
1319 0 : String agent_id(downcase(agents.name(i)));
1320 : // cerr << i << " " << agent_id << endl;
1321 0 : printAgentRecord(agent_id, i, agents.asRecord(i));
1322 0 : }
1323 0 : }
1324 :
1325 : /*
1326 : Well, for boolean values this function gives you 1 + 1 = 2,
1327 : where the AIPS++ built-in, never-should-have-been-written, sum()
1328 : surprisingly gives you 1 + 1 = 1.
1329 : */
1330 0 : int Flagger::my_aipspp_sum(const Array<Bool> &a)
1331 : {
1332 0 : return a.contiguousStorage() ?
1333 0 : std::accumulate(a.cbegin(), a.cend(), 0) :
1334 0 : std::accumulate(a.begin(), a.end(), 0);
1335 : }
1336 :
1337 0 : void Flagger::printAgentRecord(String &agent_id, uInt agentCount,
1338 : const RecordInterface &agent_rec){
1339 : // but if an id field is set in the sub-record, use that instead
1340 0 : if ( agent_rec.isDefined("id") && agent_rec.dataType("id") == TpString ){
1341 0 : agent_id = agent_rec.asString("id");
1342 : }
1343 0 : for(uInt i=0; i<agent_rec.nfields(); i++){
1344 0 : os << agent_id << "[" << agentCount+1 << "] : ";
1345 0 : String myName(agent_rec.name(i));
1346 0 : os << myName << ": ";
1347 0 : switch(agent_rec.type(i)){
1348 0 : case TpRecord :
1349 0 : printAgentRecord(myName, i, agent_rec.asRecord(i));
1350 0 : break;
1351 0 : case TpArrayBool :
1352 0 : os << agent_rec.asArrayBool(i);
1353 0 : break;
1354 0 : case TpArrayUChar :
1355 0 : os << agent_rec.asArrayuChar(i);
1356 0 : break;
1357 0 : case TpArrayShort:
1358 0 : os << agent_rec.asArrayShort(i);
1359 0 : break;
1360 0 : case TpArrayInt:
1361 0 : os << agent_rec.asArrayInt(i);
1362 0 : break;
1363 0 : case TpArrayUInt:
1364 0 : os << agent_rec.asArrayuInt(i);
1365 0 : break;
1366 0 : case TpArrayFloat:
1367 0 : os << agent_rec.asArrayFloat(i);
1368 0 : break;
1369 0 : case TpArrayDouble:
1370 0 : os << agent_rec.asArrayDouble(i);
1371 0 : break;
1372 0 : case TpArrayComplex:
1373 0 : os << agent_rec.asArrayComplex(i);
1374 0 : break;
1375 0 : case TpArrayDComplex:
1376 0 : os << agent_rec.asArrayDComplex(i);
1377 0 : break;
1378 0 : case TpArrayString:
1379 0 : os << agent_rec.asArrayString(i);
1380 0 : break;
1381 0 : case TpBool:
1382 0 : os << agent_rec.asBool(i);
1383 0 : break;
1384 0 : case TpUChar:
1385 0 : os << agent_rec.asuChar(i);
1386 0 : break;
1387 0 : case TpShort:
1388 0 : os << agent_rec.asShort(i);
1389 0 : break;
1390 0 : case TpInt:
1391 0 : os << agent_rec.asInt(i);
1392 0 : break;
1393 0 : case TpUInt:
1394 0 : os << agent_rec.asuInt(i);
1395 0 : break;
1396 0 : case TpFloat:
1397 0 : os << agent_rec.asFloat(i);
1398 0 : break;
1399 0 : case TpDouble:
1400 0 : os << agent_rec.asDouble(i);
1401 0 : break;
1402 0 : case TpComplex:
1403 0 : os << agent_rec.asComplex(i);
1404 0 : break;
1405 0 : case TpDComplex:
1406 0 : os << agent_rec.asDComplex(i);
1407 0 : break;
1408 0 : case TpString:
1409 0 : os << agent_rec.asString(i);
1410 0 : break;
1411 0 : default :
1412 0 : break;
1413 : }
1414 0 : os << endl << LogIO::POST;
1415 0 : }
1416 : //
1417 0 : }
1418 :
1419 : // -----------------------------------------------------------------------
1420 : // Flagger::run
1421 : // Performs the actual flagging
1422 : // -----------------------------------------------------------------------
1423 0 : Record Flagger::run (Bool trial, Bool reset)
1424 : {
1425 0 : if (!agents_p)
1426 : {
1427 0 : agentCount_p = 0;
1428 0 : agents_p = new Record;
1429 : if (dbg) cout << "creating new EMPTY agent and returning" << endl;
1430 0 : return Record();
1431 : }
1432 0 : Record agents = *agents_p;
1433 :
1434 : if (dbg) cout << agents << endl;
1435 :
1436 0 : if (!opts_p)
1437 : {
1438 0 : opts_p = new Record();
1439 : }
1440 0 : *opts_p = defaultOptions();
1441 0 : opts_p->define(RF_RESET,reset);
1442 0 : opts_p->define(RF_TRIAL,trial);
1443 0 : Record opt = *opts_p;
1444 :
1445 0 : if (!setdata_p) {
1446 : os << LogIO::SEVERE << "Please run setdata with/without arguments before any setmethod"
1447 0 : << LogIO::POST;
1448 0 : return Record();
1449 : }
1450 :
1451 :
1452 0 : Record result;
1453 :
1454 : //printflagselections();
1455 :
1456 0 : if ( !nant )
1457 0 : os<<"No Measurement Set has been attached\n"<<LogIO::EXCEPTION;
1458 :
1459 0 : RFABase::setIndexingBase(0);
1460 : // set debug level
1461 0 : Int debug_level = 0;
1462 0 : if ( opt.isDefined("debug") )
1463 0 : debug_level = opt.asInt("debug");
1464 :
1465 : // reset existing flags?
1466 0 : Bool reset_flags = isFieldSet(opt, RF_RESET);
1467 :
1468 : /* Don't use the progmeter if less than this
1469 : number of timestamps (for performance reasons;
1470 : just creating a ProgressMeter is relatively
1471 : expensive). */
1472 0 : const unsigned progmeter_limit = 1000;
1473 :
1474 : try { // all exceptions to be caught below
1475 :
1476 0 : bool didSomething=0;
1477 : // create iterator, visbuffer & chunk manager
1478 : // Block<Int> sortCol(1);
1479 : // sortCol[0] = MeasurementSet::SCAN_NUMBER;
1480 : //sortCol[0] = MeasurementSet::TIME;
1481 : // Setdata already made a data selection
1482 :
1483 0 : VisibilityIterator &vi(vs_p->iter());
1484 0 : vi.slurp();
1485 0 : VisBuffer vb(vi);
1486 :
1487 : RFChunkStats chunk(vi, vb,
1488 0 : *this);
1489 :
1490 : // setup global options for flagging agents
1491 0 : Record globopt(Record::Variable);
1492 0 : if ( opt.isDefined(RF_GLOBAL) )
1493 0 : globopt = opt.asRecord(RF_GLOBAL);
1494 :
1495 : // generate new array of agents by iterating through agents record
1496 0 : Record agcounts; // record of agent instance counts
1497 :
1498 0 : acc.resize(agents.nfields());
1499 :
1500 0 : for (uInt i=0; i<agents.nfields(); i++) {
1501 0 : acc[i] = std::shared_ptr<RFABase>();
1502 : }
1503 :
1504 0 : uInt nacc = 0;
1505 0 : bool only_selector = true; // only RFASelector agents?
1506 0 : for (uInt i=0; i<agents.nfields(); i++) {
1507 0 : if ( agents.dataType(i) != TpRecord )
1508 0 : os << "Unrecognized field '" << agents.name(i) << "' in agents\n" << LogIO::EXCEPTION;
1509 :
1510 0 : const RecordInterface & agent_rec( agents.asRecord(i) );
1511 :
1512 : // normally, the field name itself is the agent ID
1513 :
1514 0 : String agent_id( downcase(agents.name(i)) );
1515 :
1516 : // but if an id field is set in the sub-record, use that instead
1517 0 : if ( agent_rec.isDefined("id") && agent_rec.dataType("id") == TpString )
1518 : {
1519 0 : agent_id = agent_rec.asString("id");
1520 : }
1521 : // check that this is agent really exists
1522 0 : if ( !agent_defaults.isDefined(agent_id) )
1523 : {
1524 : //cerr << agent_defaults;
1525 : os << "Unknown flagging method '" <<
1526 0 : agents.name(i) << "'\n" << LogIO::EXCEPTION;
1527 : }
1528 :
1529 : // create parameter record by taking agent defaults, and merging in global
1530 : // and specified options
1531 0 : const RecordInterface & defparms(agent_defaults.asRecord(agent_id));
1532 0 : Record parms(defparms);
1533 0 : parms.merge(globopt,Record::OverwriteDuplicates);
1534 0 : parms.merge(agent_rec,Record::OverwriteDuplicates);
1535 :
1536 : // add the global reset argumnent
1537 0 : parms.define(RF_RESET,reset_flags);
1538 :
1539 : // see if this is a different instance of an already activated agent
1540 0 : if (agcounts.isDefined(agent_id)) {
1541 : // increment the instance counter
1542 0 : Int count = agcounts.asInt(agent_id)+1;
1543 0 : agcounts.define(agent_id,count);
1544 :
1545 : // modify the agent name to include an instance count
1546 : char s[1024];
1547 0 : sprintf(s,"%s#%d",defparms.asString(RF_NAME).chars(),count);
1548 0 : parms.define(RF_NAME,s);
1549 : }
1550 : else
1551 0 : agcounts.define(agent_id,1);
1552 : // create agent based on name
1553 : std::shared_ptr<RFABase> agent =
1554 : createAgent(agent_id,
1555 : chunk,
1556 : parms,
1557 0 : only_selector);
1558 0 : if ( agent.get() == NULL )
1559 0 : os<<"Unrecognized method name '"<<agents.name(i)<<"'\n"<<LogIO::EXCEPTION;
1560 0 : agent->init();
1561 0 : agent->setNAgent(agents.nfields());
1562 0 : String inp,st;
1563 : // agent->logSink()<<agent->getDesc()<<endl<<LogIO::POST;
1564 0 : acc[nacc++] = agent;
1565 0 : }
1566 :
1567 0 : for (unsigned i = 0; i < nacc; i++) {
1568 0 : acc[i]->setOnlySelector(only_selector);
1569 : }
1570 0 : acc.resize(nacc);
1571 :
1572 : // begin iterating over chunks
1573 0 : uInt nchunk=0;
1574 :
1575 0 : bool new_field_spw = true; /* Is the current chunk the first chunk for this (field,spw)? */
1576 :
1577 0 : Int64 inRowFlags=0, outRowFlags=0, totalRows=0, inDataFlags=0, outDataFlags=0, totalData=0;
1578 :
1579 0 : for (vi.originChunks();
1580 0 : vi.moreChunks();
1581 : ) { //vi.nextChunk(), nchunk++) {
1582 : //Start of loop over chunks
1583 0 : didSomething = false;
1584 0 : for (uInt i = 0; i < acc.size(); i++) {
1585 0 : acc[i]->initialize();
1586 : }
1587 :
1588 0 : chunk.newChunk(quack_agent_exists);
1589 :
1590 : // limit frequency of progmeter updates
1591 0 : Int pm_update_freq = chunk.num(TIME)/200;
1592 :
1593 : // How much memory do we have?
1594 0 : Int availmem = opt.isDefined("maxmem") ?
1595 0 : opt.asInt("maxmem") : HostInfo::memoryTotal()/1024;
1596 0 : if ( debug_level>0 )
1597 0 : dprintf(os,"%d MB memory available\n",availmem);
1598 :
1599 : // call newChunk() for all accumulators; determine which ones are active
1600 0 : Vector<Int> iter_mode(acc.size(),RFA::DATA);
1601 0 : Vector<Bool> active(acc.size());
1602 :
1603 0 : for (uInt i = 0; i < acc.size(); i++)
1604 : {
1605 : Int maxmem;
1606 0 : maxmem = availmem;
1607 0 : if ( ! (active(i) = acc[i]->newChunk(maxmem)) ) // refused this chunk?
1608 : {
1609 0 : iter_mode(i) = RFA::STOP; // skip over it
1610 : }
1611 : else
1612 : { // active, so reserve its memory
1613 0 : if ( debug_level>0 )
1614 0 : dprintf(os,"%s reserving %d MB of memory, %d left in pool\n",
1615 0 : acc[i]->name().chars(),availmem-maxmem,maxmem);
1616 0 : availmem = maxmem>0 ? maxmem : 0;
1617 : }
1618 : }
1619 : if (dbg) cout << "Active for this chunk: " << my_aipspp_sum(active) << endl;
1620 :
1621 : // initially active agents
1622 0 : Vector<Bool> active_init = active;
1623 : // start executing passes
1624 : char subtitle[1024];
1625 0 : sprintf(subtitle,"Flagging %s chunk %d: ",ms.tableName().chars(),nchunk+1);
1626 0 : String title(subtitle);
1627 : //cout << "--------title=" << title << endl;
1628 :
1629 0 : if ( !sum(active) )
1630 : {
1631 : //os<<LogIO::WARN<<"Unable to process this chunk with any active method.\n"<<LogIO::POST;
1632 :
1633 0 : goto end_of_loop;
1634 : /* Oh no, an evil goto statement...
1635 : This used to be a
1636 : continue;
1637 : which (in this context) is just as evil. But goto is used because some
1638 : looping code (the visIter update) needs to be always executed at the
1639 : end of this loop.
1640 :
1641 : The good solution would be to
1642 : - make the remainder of this loop
1643 : work also in the special case with sum(active) == 0
1644 : - simplify this overly long loop
1645 : */
1646 : }
1647 :
1648 0 : for( uInt npass=0; anyNE(iter_mode,(Int)RFA::STOP); npass++ ) // repeat passes while someone is active
1649 : {
1650 0 : uInt itime=0;
1651 0 : chunk.newPass(npass);
1652 : // count up who wants a data pass and who wants a dry pass
1653 0 : Int ndata = sum(iter_mode==(Int)RFA::DATA);
1654 0 : Int ndry = sum(iter_mode==(Int)RFA::DRY);
1655 0 : Int nactive = ndata+ndry;
1656 0 : if ( !nactive ) // no-one? break out then
1657 0 : break;
1658 : // didSomething++;
1659 : // Decide when to schedule a full data iteration, and when do dry runs only.
1660 : // There's probably room for optimizations here, but let's keep it simple
1661 : // for now: since data iterations are more expensive, hold them off as long
1662 : // as someone is requesting a dry run.
1663 0 : Bool data_pass = !ndry;
1664 : // Doing a full data iteration
1665 0 : if ( data_pass )
1666 : {
1667 :
1668 0 : sprintf(subtitle,"pass %d (data)",npass+1);
1669 :
1670 0 : std::unique_ptr<ProgressMeter> progmeter;
1671 0 : if (chunk.num(TIME) > progmeter_limit) {
1672 0 : progmeter = std::unique_ptr<ProgressMeter>(new ProgressMeter(1.0,static_cast<Double>(chunk.num(TIME)+0.001),title+subtitle,"","","",true,pm_update_freq));
1673 : }
1674 :
1675 : // start pass for all active agents
1676 : //cout << "-----------subtitle=" << subtitle << endl;
1677 0 : for( uInt ival = 0; ival<acc.size(); ival++ )
1678 0 : if ( active(ival) ) {
1679 0 : if ( iter_mode(ival) == RFA::DATA )
1680 0 : acc[ival]->startData(new_field_spw);
1681 0 : else if ( iter_mode(ival) == RFA::DRY )
1682 0 : acc[ival]->startDry(new_field_spw);
1683 : }
1684 : // iterate over visbuffers
1685 0 : for( vi.origin(); vi.more() && nactive; vi++,itime++ ) {
1686 :
1687 0 : if (progmeter.get() != NULL) {
1688 0 : progmeter->update(itime);
1689 : }
1690 0 : chunk.newTime();
1691 0 : Bool anyActive = false;
1692 0 : anyActive=false;
1693 0 : for( uInt i = 0; i<acc.size(); i++ )
1694 : {
1695 : //if ((acc[i]->getID() != "FlagExaminer") &&
1696 0 : if (active_init(i)) {
1697 0 : anyActive=true;
1698 : }
1699 : }
1700 :
1701 0 : for(uInt i=0;i<acc.size();i++) {
1702 0 : if (anyActive) acc[i]->initializeIter(itime);
1703 : }
1704 :
1705 0 : inRowFlags += my_aipspp_sum(vb.flagRow());
1706 0 : totalRows += vb.flagRow().nelements();
1707 0 : totalData += vb.flagCube().shape().product();
1708 :
1709 0 : inDataFlags += my_aipspp_sum(vb.flagCube());
1710 :
1711 : // now, call individual VisBuffer iterators
1712 0 : for( uInt ival = 0; ival<acc.size(); ival++ )
1713 0 : if ( active(ival) ) {
1714 : // call iterTime/iterDry as appropriate
1715 0 : RFA::IterMode res = RFA::STOP;
1716 0 : if ( iter_mode(ival) == RFA::DATA ) {
1717 0 : res = acc[ival]->iterTime(itime);
1718 : }
1719 0 : else if ( iter_mode(ival) == RFA::DRY )
1720 0 : res = acc[ival]->iterDry(itime);
1721 : // change requested? Deactivate agent
1722 0 : if ( ! ( res == RFA::CONT || res == iter_mode(ival) ) )
1723 : {
1724 0 : active(ival) = false;
1725 0 : nactive--;
1726 0 : iter_mode(ival)==RFA::DATA ? ndata-- : ndry--;
1727 0 : iter_mode(ival) = res;
1728 0 : if ( nactive <= 0 )
1729 0 : break;
1730 : }
1731 : }
1732 :
1733 : // also iterate over rows for data passes
1734 0 : for( Int ir=0; ir<vb.nRow() && ndata; ir++ ) {
1735 0 : for( uInt ival = 0; ival<acc.size(); ival++ )
1736 0 : if ( iter_mode(ival) == RFA::DATA )
1737 : {
1738 0 : RFA::IterMode res = acc[ival]->iterRow(ir);
1739 0 : if ( ! ( res == RFA::CONT || res == RFA::DATA ) )
1740 : {
1741 0 : ndata--; nactive--;
1742 0 : iter_mode(ival) = res;
1743 0 : active(ival) = false;
1744 0 : if ( ndata <= 0 )
1745 0 : break;
1746 : }
1747 : }
1748 : }
1749 :
1750 0 : for( uInt ival = 0; ival<acc.size(); ival++ ) {
1751 0 : if ( active(ival) ) {
1752 0 : acc[ival]->endRows(itime);
1753 : }
1754 : }
1755 : } /* for vi... */
1756 :
1757 : // end pass for all agents
1758 0 : for( uInt ival = 0; ival<acc.size(); ival++ )
1759 : {
1760 0 : if ( active(ival) ) {
1761 0 : if ( iter_mode(ival) == RFA::DATA )
1762 0 : iter_mode(ival) = acc[ival]->endData();
1763 0 : else if ( iter_mode(ival) == RFA::DRY )
1764 0 : iter_mode(ival) = acc[ival]->endDry();
1765 : }
1766 : }
1767 0 : }
1768 : else // dry pass only
1769 : {
1770 0 : sprintf(subtitle,"pass %d (dry)",npass+1);
1771 : //cout << "-----------subtitle=" << subtitle << endl;
1772 :
1773 0 : std::unique_ptr<ProgressMeter> progmeter;
1774 0 : if (chunk.num(TIME) > progmeter_limit) {
1775 0 : progmeter = std::unique_ptr<ProgressMeter>(new ProgressMeter (1.0,static_cast<Double>(chunk.num(TIME)+0.001),title+subtitle,"","","",true,pm_update_freq));
1776 : }
1777 : // start pass for all active agents
1778 0 : for( uInt ival = 0; ival<acc.size(); ival++ )
1779 0 : if ( iter_mode(ival) == RFA::DRY )
1780 0 : acc[ival]->startDry(new_field_spw);
1781 0 : for( uInt itime=0; itime<chunk.num(TIME) && ndry; itime++ )
1782 : {
1783 0 : if (progmeter.get() != NULL) {
1784 0 : progmeter->update(itime);
1785 : }
1786 : // now, call individual VisBuffer iterators
1787 0 : for( uInt ival = 0; ival<acc.size(); ival++ )
1788 0 : if ( iter_mode(ival) == RFA::DRY )
1789 : {
1790 : // call iterTime/iterDry as appropriate
1791 0 : RFA::IterMode res = acc[ival]->iterDry(itime);
1792 : // change requested? Deactivate agent
1793 0 : if ( ! ( res == RFA::CONT || res == RFA::DRY ) )
1794 : {
1795 0 : iter_mode(ival) = res;
1796 0 : active(ival) = false;
1797 0 : if ( --ndry <= 0 )
1798 0 : break;
1799 : }
1800 : }
1801 : }
1802 : // end pass for all agents
1803 0 : for( uInt ival = 0; ival<acc.size(); ival++ )
1804 0 : if ( iter_mode(ival) == RFA::DRY )
1805 0 : iter_mode(ival) = acc[ival]->endDry();
1806 0 : } // end of dry pass
1807 : } // end loop over passes
1808 :
1809 : //cout << opt << endl;
1810 : //cout << "any active = " << active_init << endl;
1811 :
1812 0 : if ( !isFieldSet(opt, RF_TRIAL) && anyNE(active_init, false) )
1813 : {
1814 0 : sprintf(subtitle,"pass (flag)");
1815 : //cout << "-----------subtitle=" << subtitle << endl;
1816 :
1817 0 : std::unique_ptr<ProgressMeter> progmeter;
1818 0 : if (chunk.num(TIME) > progmeter_limit) {
1819 0 : progmeter = std::unique_ptr<ProgressMeter>(new ProgressMeter(1.0,static_cast<Double>(chunk.num(TIME)+0.001),title+"storing flags","","","",true,pm_update_freq));
1820 : }
1821 0 : for (uInt i = 0; i<acc.size(); i++)
1822 0 : if (active_init(i))
1823 0 : acc[i]->startFlag(new_field_spw);
1824 0 : uInt itime=0;
1825 0 : for( vi.origin(); vi.more(); vi++,itime++ ) {
1826 0 : if (progmeter.get() != NULL) {
1827 0 : progmeter->update(itime);
1828 : }
1829 :
1830 0 : chunk.newTime();
1831 : // inRowFlags += sum(chunk.nrfIfr());
1832 :
1833 0 : Bool anyActive = false;
1834 0 : anyActive=false;
1835 0 : for( uInt i = 0; i<acc.size(); i++ ) {
1836 : // cout << i << " " << acc[i]->getID() << " " << active_init(i) << endl;
1837 0 : if ((acc[i]->getID() != "FlagExaminer") &&
1838 0 : active_init(i))
1839 0 : anyActive=true;
1840 : }
1841 :
1842 : //cout << "anyActive" << anyActive << endl;
1843 :
1844 0 : didSomething = (anyActive==true);
1845 0 : for( uInt i = 0; i<acc.size(); i++ ) {
1846 0 : if ( active_init(i) ) {
1847 : //if (acc[i]->getID() != "FlagExaminer" )
1848 0 : acc[i]->iterFlag(itime);
1849 : }
1850 0 : if (anyActive) acc[i]->finalizeIter(itime);
1851 : }
1852 :
1853 : // outRowFlags += sum(chunk.nrfIfr());
1854 :
1855 0 : outRowFlags += my_aipspp_sum(vb.flagRow());
1856 0 : outDataFlags += my_aipspp_sum(vb.flagCube());
1857 :
1858 : } // for (vi ... ) loop over time
1859 0 : if (didSomething) {
1860 0 : for (uInt i = 0; i < acc.size(); i++) {
1861 0 : if (acc[i].get()) acc[i]->finalize();
1862 : }
1863 : }
1864 0 : for (uInt i = 0; i < acc.size(); i++) {
1865 :
1866 0 : if ( active_init(i) ) {
1867 :
1868 0 : acc[i]->endFlag();
1869 : // summary mode prints here
1870 :
1871 : // cerr << "Agent = " << acc[i]->getID() << endl;
1872 : // cerr << "Agent's name = " << acc[i]->name() << endl;
1873 : }
1874 : }
1875 : {
1876 0 : logSink_p.clearLocally();
1877 0 : LogIO oss(LogOrigin("Flagger", "run()"), logSink_p);
1878 0 : os=oss;
1879 0 : }
1880 0 : } /* end if not trial run and some agent is active */
1881 :
1882 0 : end_of_loop:
1883 :
1884 : // call endChunk on all agents
1885 0 : for( uInt i = 0; i<acc.size(); i++ )
1886 0 : acc[i]->endChunk();
1887 :
1888 0 : int field_id = chunk.visBuf().fieldId();
1889 0 : int spw_id = chunk.visBuf().spectralWindow();
1890 0 : int scan_number = chunk.visBuf().scan()(0);
1891 :
1892 :
1893 : /* Is this the end of the current (field, spw)? */
1894 0 : new_field_spw = ( !vi.moreChunks() ||
1895 0 : !scan_looping ||
1896 0 : ! (chunk.visBuf().fieldId() == field_id &&
1897 0 : chunk.visBuf().spectralWindow() == spw_id &&
1898 0 : chunk.visBuf().scan()(0) > scan_number )
1899 : );
1900 :
1901 0 : if (didSomething && new_field_spw) {
1902 :
1903 0 : LogIO osss(LogOrigin("Flagger", "run"),logSink_p);
1904 :
1905 0 : stringstream ss;
1906 0 : ss << "Chunk = " << chunk.nchunk()
1907 0 : << ", Field = " << field_id << " (" << chunk.visIter().fieldName() << ")"
1908 0 : << ", Spw Id = " << spw_id << ", Corrs = [" << chunk.getCorrString() << "]"
1909 0 : << ", nTime = " << vi.nSubInterval() << ", Total rows = " << totalRows
1910 0 : << ", Total data = " << totalData << endl;
1911 : ss << "Input: "
1912 0 : << " Rows flagged = " << inRowFlags << " "
1913 0 : << "(" << 100.0*inRowFlags/totalRows << " %)."
1914 0 : << " Data flagged = " << inDataFlags << " "
1915 0 : << "(" << 100.0*inDataFlags/totalData << " %)."
1916 0 : << endl;
1917 : ss << "This run: "
1918 0 : << " Rows flagged = " << outRowFlags - inRowFlags << " "
1919 0 : << "(" << 100.0*(outRowFlags-inRowFlags)/totalRows << " %)."
1920 0 : << " Data flagged = " << outDataFlags - inDataFlags << " "
1921 0 : << "(" << 100.0*(outDataFlags-inDataFlags)/totalData << " %)."
1922 0 : << endl;
1923 0 : ss << "------------------------------------------------------------------------------------" << endl;
1924 :
1925 0 : osss << ss.str() << LogIO::POST;
1926 :
1927 : /* Reset counters */
1928 0 : inRowFlags = outRowFlags = totalRows = inDataFlags = outDataFlags = totalData = 0;
1929 0 : }
1930 :
1931 : // CAS-2798: moved vi.nextChunk() to here to avoid printing the wrong information above
1932 0 : vi.nextChunk(); nchunk++;
1933 :
1934 0 : } // end loop over chunks
1935 :
1936 : // get results for all agents
1937 : // (gets just last agent; doesn't work for multiple agents,
1938 : // but only used in mode summary)
1939 0 : for (uInt i = 0; i < acc.size(); i++) {
1940 0 : result = acc[i]->getResult();
1941 : }
1942 :
1943 : if (dbg) {
1944 : cout << "Total number of data chunks : " << nchunk << endl;
1945 : }
1946 :
1947 0 : }
1948 0 : catch(AipsError &x)
1949 : {
1950 0 : throw;
1951 0 : }
1952 0 : catch(std::exception &e)
1953 : {
1954 0 : throw AipsError(e.what());
1955 0 : }
1956 :
1957 0 : ms.flush();
1958 : //os<<"Flagging complete\n"<<LogIO::POST;
1959 :
1960 : /* Clear the current flag selections */
1961 0 : clearflagselections(0);
1962 :
1963 0 : return result;
1964 0 : }
1965 :
1966 : // -----------------------------------------------------------------------
1967 : // printSummaryReport
1968 : // Generates a summary flagging report for current chunk
1969 : // -----------------------------------------------------------------------
1970 0 : void Flagger::printSummaryReport (RFChunkStats &chunk)
1971 : {
1972 : if (dbg) cout << "Flagger:: printSummaryReport" << endl;
1973 : // generate a short text report in the first pane
1974 : char s[1024];
1975 0 : sprintf(s,"MS '%s'\nchunk %d (field %s, spw %d)",ms.tableName().chars(),
1976 0 : chunk.nchunk(),chunk.visIter().fieldName().chars(),chunk.visIter().spectralWindow());
1977 0 : os << "---------------------------------------------------------------------" << LogIO::POST;
1978 0 : os<<s<<LogIO::POST;
1979 :
1980 : // print overall flagging stats
1981 0 : uInt n=0,n0;
1982 :
1983 0 : sprintf(s,"%s, %d channels, %d time slots, %d baselines, %d rows\n",
1984 0 : chunk.getCorrString().chars(),chunk.num(CHAN),chunk.num(TIME),
1985 : chunk.num(IFR),chunk.num(ROW));
1986 0 : os<<s<<LogIO::POST;
1987 :
1988 : // % of rows flagged
1989 0 : n = sum(chunk.nrfIfr());
1990 0 : n0 = chunk.num(ROW);
1991 0 : sprintf(s,"%d (%0.2f%%) rows are flagged (all baselines/times/chans/corrs in this chunk).",n,n*100.0/n0);
1992 0 : os<<s<<LogIO::POST;
1993 :
1994 : // % of data points flagged
1995 0 : n = sum(chunk.nfIfrTime());
1996 0 : n0 = chunk.num(ROW)*chunk.num(CHAN)*chunk.num(CORR);
1997 0 : sprintf(s,"%d of %d (%0.2f%%) data points are flagged (all baselines/times/chans/corrs in this chunk).",n,n0,n*100.0/n0);
1998 0 : os<<s<<LogIO::POST;
1999 : //os << "---------------------------------------------------------------------" << LogIO::POST;
2000 :
2001 : // % flagged per baseline (ifr)
2002 : // % flagged per timestep (itime)
2003 : // // there is info about (ifr,itime)
2004 : // % flagged per antenna (ifr) -> decompose into a1,a2
2005 : // % flagged per channel/corr -> (ich,ifr)
2006 : // % flagged per ( field, spw, array, scan ) are chunks - i think !!!
2007 :
2008 :
2009 : // print per-agent flagging summary
2010 : /*
2011 : for( uInt i=0; i<acc.size(); i++ )
2012 : {
2013 : String name(acc[i]->name() + "["+i+"]"+": ");
2014 : String stats( acc[i]->isActive() ? acc[i]->getStats() : String("can't process this chunk") );
2015 : os<<name+stats<<LogIO::POST;
2016 : }
2017 : */
2018 : if (dbg) cout << "end of.... Flagger:: printSummaryReport" << endl;
2019 0 : }
2020 :
2021 : // -----------------------------------------------------------------------
2022 : // printAgentReport
2023 : // Generates per-agent reports for current chunk of data
2024 : // Meant to be called before doing endChunk() on all the flagging
2025 : // agents.
2026 : // -----------------------------------------------------------------------
2027 0 : void Flagger::printAgentReports( )
2028 : {
2029 : // call each agent to produce summary plots
2030 0 : for( uInt i=0; i<acc.size(); i++ )
2031 0 : acc[i]->printFlaggingReport();
2032 0 : }
2033 :
2034 : /* FLAG VERSION SUPPORT */
2035 0 : Bool Flagger::saveFlagVersion(String versionname, String comment, String merge )
2036 : {
2037 : try
2038 : {
2039 0 : FlagVersion fv(originalms.tableName(),"FLAG","FLAG_ROW");
2040 0 : fv.saveFlagVersion(versionname, comment, merge);
2041 0 : }
2042 0 : catch (AipsError x)
2043 : {
2044 0 : os << LogIO::SEVERE << "Could not save Flag Version : " << x.getMesg() << LogIO::POST;
2045 0 : throw;
2046 0 : }
2047 0 : return true;
2048 : }
2049 0 : Bool Flagger::restoreFlagVersion(Vector<String> versionname, String merge )
2050 : {
2051 : try
2052 : {
2053 0 : FlagVersion fv(originalms.tableName(),"FLAG","FLAG_ROW");
2054 0 : for(Int j=0;j<(Int)versionname.nelements();j++)
2055 0 : fv.restoreFlagVersion(versionname[j], merge);
2056 0 : }
2057 0 : catch (AipsError x)
2058 : {
2059 0 : os << LogIO::SEVERE << "Could not restore Flag Version : " << x.getMesg() << LogIO::POST;
2060 0 : return false;
2061 0 : }
2062 0 : return true;
2063 : }
2064 0 : Bool Flagger::deleteFlagVersion(Vector<String> versionname)
2065 : {
2066 : try
2067 : {
2068 0 : FlagVersion fv(originalms.tableName(),"FLAG","FLAG_ROW");
2069 0 : for(Int j=0;j<(Int)versionname.nelements();j++)
2070 0 : fv.deleteFlagVersion(versionname[j]);
2071 0 : }
2072 0 : catch (AipsError x)
2073 : {
2074 0 : os << LogIO::SEVERE << "Could not delete Flag Version : " << x.getMesg() << LogIO::POST;
2075 0 : return false;
2076 0 : }
2077 0 : return true;
2078 : }
2079 0 : Bool Flagger::getFlagVersionList(Vector<String> &verlist)
2080 : {
2081 : try
2082 : {
2083 0 : verlist.resize(0);
2084 : Int num;
2085 0 : FlagVersion fv(originalms.tableName(),"FLAG","FLAG_ROW");
2086 0 : Vector<String> vlist = fv.getVersionList();
2087 :
2088 0 : num = verlist.nelements();
2089 0 : verlist.resize( num + vlist.nelements() + 1, true );
2090 0 : verlist[num] = String("\nMS : ") + originalms.tableName() + String("\n");
2091 0 : for(Int j=0;j<(Int)vlist.nelements();j++)
2092 0 : verlist[num+j+1] = vlist[j];
2093 0 : }
2094 0 : catch (AipsError x)
2095 : {
2096 0 : os << LogIO::SEVERE << "Could not get Flag Version List : " << x.getMesg() << LogIO::POST;
2097 0 : return false;
2098 0 : }
2099 0 : return true;
2100 : }
2101 :
2102 :
2103 :
2104 :
2105 0 : void Flagger::reform_baselinelist(Matrix<Int> &baselinelist, unsigned nant)
2106 : {
2107 0 : for (unsigned i = 0; (int)i < baselinelist.shape()(1); i++) {
2108 0 : int a1 = baselinelist(0, i);
2109 0 : if (a1 < 0) {
2110 0 : for (unsigned a = 0; a < nant; a++) {
2111 0 : if (a != (unsigned) (-a1)) {
2112 0 : IPosition longer = baselinelist.shape();
2113 0 : longer(1) += 1;
2114 0 : baselinelist.resize(longer);
2115 0 : baselinelist(0, longer(1)-1) = a;
2116 0 : baselinelist(1, longer(1)-1) = baselinelist(1, i);
2117 0 : }
2118 : }
2119 : }
2120 : }
2121 0 : }
2122 :
2123 : // -----------------------------------------------------------------------
2124 : // dprintf
2125 : // Function for printfing stuff to a debug stream
2126 : // -----------------------------------------------------------------------
2127 0 : int dprintf( LogIO &os,const char *format, ...)
2128 : {
2129 : char str[1024];
2130 : va_list ap;
2131 0 : va_start(ap, format);
2132 0 : int ret = vsnprintf(str, 1024, format, ap);
2133 0 : va_end(ap);
2134 0 : os << LogIO::DEBUGGING << str << LogIO::POST;
2135 0 : return ret;
2136 : }
2137 :
2138 : } //#end casa namespace
|
