Line data Source code
1 :
2 : //# RFASelector.cc: this defines RFASelector
3 : //# Copyright (C) 2000,2001,2002
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: casa-feedback@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : #include <casacore/casa/Exceptions/Error.h>
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <casacore/casa/Arrays/ArrayLogical.h>
31 : #include <casacore/casa/Arrays/MaskedArray.h>
32 : #include <casacore/casa/Arrays/MaskArrMath.h>
33 : #include <casacore/casa/Quanta/Quantum.h>
34 : #include <casacore/casa/Quanta/MVTime.h>
35 : #include <casacore/casa/Logging/LogIO.h>
36 : #include <msvis/MSVis/VisibilityIterator.h>
37 : #include <msvis/MSVis/VisBuffer.h>
38 : #include <flagging/Flagging/RFASelector.h>
39 :
40 : #include <iomanip>
41 : #include <ostream>
42 : #include <cassert>
43 :
44 : using namespace casacore;
45 : namespace casa { //# NAMESPACE CASA - BEGIN
46 :
47 : Bool dbg2 = false;
48 : Bool verbose2 = false;
49 :
50 : // -----------------------------------------------------------------------
51 : // reformRange
52 : // Reforms an array of 2N elements into a [2,N] matrix
53 : // -----------------------------------------------------------------------
54 :
55 0 : template<> Array<Int> fieldToArray<Int>( const RecordInterface &parm,const String &id )
56 0 : { return parm.toArrayInt(id); }
57 0 : template<> Array<Double> fieldToArray<Double>( const RecordInterface &parm,const String &id )
58 0 : { return parm.toArrayDouble(id); }
59 0 : template<> Array<String> fieldToArray<String>( const RecordInterface &parm,const String &id )
60 0 : { return parm.toArrayString(id); }
61 :
62 :
63 : template Bool RFASelector::reformRange<Int>( Matrix<Int>&,const Array<Int>& );
64 : template Bool RFASelector::reformRange<Double>( Matrix<Double>&,const Array<Double>& );
65 : template Bool RFASelector::reformRange<String>( Matrix<String>&,const Array<String>& );
66 : template Bool RFASelector::parseRange<Int>( Matrix<Int>&,const RecordInterface&,const String&);
67 : template Bool RFASelector::parseRange<Double>( Matrix<Double>&,const RecordInterface&,const String&);
68 : template Bool RFASelector::parseRange<String>( Matrix<String>&,const RecordInterface&,const String&);
69 :
70 : // -----------------------------------------------------------------------
71 : // RFA_selector::find
72 : // Helper templated method to find an object in an array
73 : // -----------------------------------------------------------------------
74 :
75 : template Bool RFASelector::find<uInt>(uInt&,const uInt&,const Vector<uInt>&);
76 : template Bool RFASelector::find<Int>(uInt&,const Int&,const Vector<Int>&);
77 : template Bool RFASelector::find<String>(uInt&,const String&,const Vector<String>&);
78 :
79 : // parseTimes
80 : // Converts a field that is an array of Int/Double or Strings into
81 : // an array of Doubles. Numeric values are converted as is. Strings
82 : // are fed through MVTime::read to convert into MJDs (or MJSs, if
83 : // secs is true).
84 : // -----------------------------------------------------------------------
85 0 : Bool RFASelector::parseTimes ( Array<Double> ×,const RecordInterface &parm,const String &id,Bool secs )
86 : {
87 0 : LogIO os(LogOrigin("RFASelector", "parseTimes()", WHERE));
88 0 : if( !isFieldSet(parm,id) )
89 0 : return false;
90 0 : if( fieldType(parm,id,TpString,TpArrayString) ) // String date/times
91 : {
92 0 : Array<String> tt( parm.asArrayString(id) );
93 0 : times.resize(tt.shape());
94 : Bool deltt,deltimes;
95 0 : const String *ptt = tt.getStorage(deltt);
96 0 : Double *ptimes = times.getStorage(deltimes);
97 0 : Int scale = secs ? 24*3600 : 1;
98 0 : for( uInt i=0; i<tt.nelements(); i++ )
99 : {
100 0 : Quantity q;
101 0 : if( !MVTime::read(q,ptt[i]) ) {
102 0 : os<<"bad "<<id<<" specified: "<<ptt[i]<<endl<<LogIO::EXCEPTION;
103 : }
104 0 : ptimes[i] = scale*(Double)MVTime(q);
105 0 : }
106 0 : tt.freeStorage(ptt,deltt);
107 0 : times.putStorage(ptimes,deltimes);
108 0 : }
109 0 : else if( isField(parm,id,isReal) ) // if not strings, try numeric MJDs
110 : {
111 0 : times = parm.toArrayDouble(id);
112 : }
113 : else // else can't parse
114 0 : return false;
115 0 : return true;
116 0 : }
117 :
118 :
119 : // -----------------------------------------------------------------------
120 : // addString
121 : // Helper method to build up description strings
122 : // -----------------------------------------------------------------------
123 0 : void RFASelector::addString( String &str,const String &s1,const char *sep )
124 : {
125 0 : if( str.length() )
126 0 : str += sep;
127 0 : str += s1;
128 0 : }
129 :
130 : // -----------------------------------------------------------------------
131 : // parseMinMax
132 : // Helper function to parse a range specification
133 : // -----------------------------------------------------------------------
134 0 : Bool RFASelector::parseMinMax( Float &vmin,Float &vmax,const RecordInterface &spec,uInt f0 )
135 : {
136 :
137 0 : vmin = -C::flt_max; vmax=C::flt_max;
138 : // Option 1: fields named min/max exist... so use them
139 0 : Bool named = false;
140 0 : if( spec.isDefined(RF_MIN) )
141 0 : { vmin = spec.asFloat(RF_MIN); named = true; }
142 0 : if( spec.isDefined(RF_MAX) )
143 0 : { vmax = spec.asFloat(RF_MAX); named = true; }
144 0 : if( named )
145 0 : return true;
146 : // Else look at first available field, if a 2-element array, assume
147 : // [min,max] has been specified
148 0 : if (spec.nfields() < f0) {
149 0 : return false;
150 : }
151 0 : if( spec.shape(f0).nelements()==1 && spec.shape(f0)(0) == 2 )
152 : {
153 0 : Vector<Double> mm = spec.toArrayDouble(f0);
154 0 : vmin=mm(0); vmax=mm(1);
155 0 : return true;
156 0 : }
157 : // Else assume next two record fields are {min,max}
158 0 : if( spec.nfields()-f0 > 2 )
159 : {
160 0 : vmin = spec.asFloat(f0);
161 0 : vmax = spec.asFloat(f0+1);
162 : }
163 : else
164 0 : vmax = spec.asFloat(f0);
165 0 : return true;
166 : }
167 :
168 0 : Bool RFASelector::fortestingonly_parseMinMax( Float &vmin,Float &vmax,const RecordInterface &spec,uInt f0 )
169 : {
170 0 : return parseMinMax( vmin, vmax, spec, f0 );
171 : }
172 :
173 :
174 : // -----------------------------------------------------------------------
175 : // normalize
176 : // Helper function to shift a cyclic value (i.e. angle) into
177 : // the interval [base,base+cycle) by adding/subtarcting an integer
178 : // number of cycles
179 : // -----------------------------------------------------------------------
180 0 : static Double normalize( Double value,Double base,Double cycle )
181 : {
182 0 : if( value < base )
183 0 : value += (floor((base-value)/cycle)+1)*cycle;
184 0 : else if( value >= base+cycle )
185 0 : value -= floor((value-base)/cycle)*cycle;
186 0 : return value;
187 : }
188 :
189 :
190 : // -----------------------------------------------------------------------
191 : // addClipInfo
192 : // -----------------------------------------------------------------------
193 0 : void RFASelector::addClipInfo( const Vector<String> &expr,
194 : Float vmin,Float vmax,
195 : Bool clip,
196 : Bool channel_average)
197 : {
198 : // create mapper and clipinfo block
199 0 : RFDataMapper *mapper = new RFDataMapper(expr,sel_column);
200 :
201 0 : ClipInfo clipinfo = { mapper,
202 : vmin,vmax,
203 : channel_average,
204 0 : clip,0.0 };
205 :
206 : // if dealing with cyclic values, normalize min/max accordingly
207 0 : Double cycle = mapper->getValueCycle(); // e.g. 360 for angles
208 :
209 0 : if (cycle > 0) {
210 0 : Double base = mapper->getValueBase(); // e.g. -180 for angles
211 :
212 : // normalize min/max angle into [base,base+cycle)
213 0 : clipinfo.vmin = normalize(clipinfo.vmin, base, cycle);
214 0 : clipinfo.vmax = normalize(clipinfo.vmax, base, cycle);
215 :
216 : // if order is reversed, then we're spanning a cycle boundary (i.e. 355->5)
217 0 : if( clipinfo.vmin > clipinfo.vmax ) {
218 0 : clipinfo.vmax += cycle; // ...so add a cycle
219 : }
220 :
221 : // use vmin as offset
222 0 : clipinfo.offset = clipinfo.vmin;
223 0 : clipinfo.vmin = 0;
224 0 : clipinfo.vmax -= clipinfo.offset;
225 : }
226 : // add block to appropriate list
227 0 : Block<ClipInfo> & block( mapper->type()==RFDataMapper::MAPROW ? sel_clip_row : sel_clip );
228 0 : uInt ncl = block.nelements();
229 0 : block.resize(ncl+1,false,true);
230 0 : block[ncl] = clipinfo;
231 0 : }
232 :
233 : // -----------------------------------------------------------------------
234 : // parseClipField
235 : // Helper function to parse a clip specification
236 : // -----------------------------------------------------------------------
237 0 : void RFASelector::parseClipField( const RecordInterface &spec,Bool clip )
238 : {
239 : //cerr << "spec: " << spec << endl;
240 : //cerr << "clip: " << clip << endl;
241 : //cerr << "spec.name(0): " << spec.name(0)<< endl;
242 : //cerr << "RF_EXPR: " << RF_EXPR << endl;
243 :
244 : try {
245 :
246 0 : bool ca = spec.asBool(RF_CHANAVG);
247 :
248 : // Syntax one - we have a record of {expr,min,max} or {expr,[min,max]}
249 : // or {expr,max}
250 0 : if( spec.name(0)== RF_EXPR )
251 : {
252 0 : Vector<String> expr = spec.asArrayString(0);
253 :
254 : Float vmin,vmax;
255 0 : if( !parseMinMax(vmin,vmax,spec,1))
256 0 : throw(AipsError(""));
257 0 : addClipInfo(expr, vmin, vmax, clip, ca);
258 0 : }
259 : // Syntax two: we have a record of { expr1=[min,max],expr2=[min,max],.. }
260 : else
261 : {
262 0 : for( uInt i=0; i<spec.nfields(); i++ )
263 : {
264 0 : Vector<String> expr(1,spec.name(i));
265 0 : Float vmin=-C::flt_max,vmax=C::flt_max;
266 0 : if( spec.dataType(i) == TpRecord )
267 : {
268 0 : uInt f0=0;
269 0 : if( spec.asRecord(i).name(0) == RF_EXPR )
270 : {
271 0 : expr = spec.asRecord(i).asArrayString(0);
272 0 : f0++;
273 : }
274 0 : if( !parseMinMax(vmin,vmax,spec.asRecord(i),f0))
275 0 : throw(AipsError(""));
276 : }
277 : else
278 : {
279 0 : if( isArray( spec.type(i)))
280 : {
281 0 : Vector<Double> vec = spec.toArrayDouble(i);
282 0 : if( vec.nelements() == 1 )
283 0 : vmax=vec(0);
284 0 : else if( vec.nelements() == 2 )
285 0 : { vmin=vec(0); vmax=vec(1); }
286 : else
287 0 : throw(AipsError(""));
288 0 : }
289 : else
290 : {
291 0 : vmax = spec.asFloat(i);
292 : }
293 : }
294 0 : addClipInfo(expr,vmin,vmax,clip, ca);
295 0 : }
296 : }
297 : }
298 0 : catch( AipsError x ) {
299 0 : os<<"Illegal \""<<(clip?RF_CLIP:RF_FLAGRANGE) <<
300 : "\" record: " << x.what() << "\n" <<
301 0 : LogIO::EXCEPTION;
302 0 : }
303 0 : }
304 :
305 0 : void RFASelector::fortestingonly_parseClipField( const RecordInterface &spec,Bool clip )
306 : {
307 0 : parseClipField(spec, clip);
308 0 : }
309 :
310 0 : void RFASelector::addClipInfoDesc ( const Block<ClipInfo> &clip)
311 : {
312 0 : for( uInt i=0; i<clip.nelements(); i++ )
313 : {
314 0 : String ss;
315 0 : char s1[32]="",s2[32]="";
316 :
317 0 : if (clip[i].channel_average) {
318 0 : ss = "average=1 ";
319 : }
320 : else {
321 0 : ss = "average=0 ";
322 : }
323 :
324 0 : if( clip[i].vmin != -C::flt_max )
325 0 : sprintf(s1,"%g",clip[i].vmin + clip[i].offset);
326 :
327 0 : if( clip[i].vmax != C::flt_max )
328 0 : sprintf(s2,"%g",clip[i].vmax + clip[i].offset);
329 :
330 0 : if( clip[i].clip )
331 : {
332 0 : ss += clip[i].mapper->description();
333 0 : if( s1[0] )
334 : {
335 0 : ss += String("<") +s1;
336 0 : if( s2[0] ) ss += ",";
337 : }
338 0 : if( s2[0] )
339 0 : ss += String(">")+s2;
340 : }
341 : else
342 : {
343 0 : if( s1[0] )
344 0 : ss += String(s1)+"<=";
345 0 : ss += clip[i].mapper->description();
346 0 : if( s2[0] )
347 0 : ss += String("<=")+s2;
348 : }
349 :
350 0 : addString(desc_str,ss);
351 0 : }
352 0 : }
353 :
354 : // -----------------------------------------------------------------------
355 : // RFASelector constructor
356 : // -----------------------------------------------------------------------
357 0 : RFASelector::RFASelector ( RFChunkStats &ch,const RecordInterface &parm) :
358 0 : RFAFlagCubeBase(ch,parm)
359 : {
360 0 : if ( chunk.measSet().isNull() ) {
361 0 : throw AipsError("Received chunk referring to NULL MS!");
362 : }
363 :
364 : char s[256];
365 :
366 0 : if( fieldType(parm,RF_FREQS,TpArrayString)) // frequency range[s], as measures
367 : {
368 0 : Matrix<String> sfq;
369 0 : parseRange(sfq,parm,RF_FREQS);
370 0 : sel_freq.resize( sfq.shape()) ;
371 : // parse array of String frequency quanta
372 0 : for( uInt i=0; i<sfq.nrow(); i++)
373 0 : for( uInt j=0; i<sfq.ncolumn(); j++)
374 : {
375 : Float q; char unit[32];
376 0 : if( sscanf(sfq(i,j).chars(),"%f%s",&q,unit)<2)
377 0 : os<<"Illegal \""<<RF_FREQS<<"\" array\n"<<LogIO::EXCEPTION;
378 0 : Quantity qq(q,unit);
379 0 : sel_freq(i,j) = qq.getValue("Hz");
380 0 : }
381 0 : }
382 : else // freq. specified as MHz
383 : {
384 0 : parseRange(sel_freq,parm,RF_FREQS);
385 0 : sel_freq *= 1e+6;
386 : }
387 0 : if( sel_freq.nelements())
388 : {
389 0 : String fq;
390 0 : for( uInt i=0; i<sel_freq.ncolumn(); i++)
391 : {
392 0 : sprintf(s,"%.2f-%.2f",sel_freq(0,i)*1e-6,sel_freq(1,i)*1e-6);
393 0 : addString(fq,s,",");
394 : }
395 0 : addString(desc_str,String(RF_FREQS)+"="+fq+"MHz");
396 0 : }
397 :
398 : // parse input arguments: channels
399 0 : if( parseRange(sel_chan,parm,RF_CHANS))
400 : {
401 0 : String sch;
402 0 : for( uInt i=0; i<sel_chan.ncolumn(); i++)
403 : {
404 0 : sprintf(s,"%d:%d",sel_chan(0,i),sel_chan(1,i));
405 0 : addString(sch,s,",");
406 : }
407 0 : addString(desc_str,String(RF_CHANS)+"="+sch);
408 0 : sel_chan(sel_chan>=0) += -(Int)indexingBase();
409 0 : }
410 :
411 : // parse input arguments: correlations
412 0 : if( fieldType(parm,RF_CORR,TpString,TpArrayString))
413 : {
414 0 : String ss;
415 0 : Vector<String> scorr( parm.asArrayString(RF_CORR)) ;
416 0 : sel_corr.resize( scorr.nelements()) ;
417 0 : for( uInt i=0; i<scorr.nelements(); i++)
418 : {
419 0 : sel_corr(i) = Stokes::type( scorr(i)) ;
420 0 : if( sel_corr(i) == Stokes::Undefined)
421 0 : os<<"Illegal correlation "<<scorr(i)<<endl<<LogIO::EXCEPTION;
422 0 : addString(ss,scorr(i),",");
423 : }
424 0 : addString(desc_str,String(RF_CORR)+"="+ss);
425 0 : }
426 : // read clip column
427 0 : if( fieldType(parm,RF_COLUMN,TpString))
428 : {
429 0 : parm.get(RF_COLUMN,sel_column);
430 0 : addString(desc_str, String(RF_COLUMN)+"="+sel_column);
431 : }
432 :
433 : // parse input arguments: Spw ID(s)
434 0 : if( fieldType(parm,RF_SPWID,TpInt,TpArrayInt))
435 : {
436 0 : parm.get(RF_SPWID,sel_spwid);
437 0 : String ss;
438 0 : for( uInt i=0; i<sel_spwid.nelements(); i++)
439 0 : addString(ss,String::toString(sel_spwid(i)),",");
440 0 : addString(desc_str,String(RF_SPWID)+"="+ss);
441 0 : sel_spwid -= (Int)indexingBase();
442 0 : }
443 :
444 : // parse input arguments: Field names or ID(s)
445 0 : if( fieldType(parm,RF_FIELD,TpString,TpArrayString))
446 : {
447 0 : parm.get(RF_FIELD,sel_fieldnames);
448 0 : sel_fieldnames.apply(stringUpper);
449 0 : String ss;
450 0 : for( uInt i=0; i<sel_fieldnames.nelements(); i++)
451 0 : addString(ss,sel_fieldnames(i),",");
452 0 : addString(desc_str,String(RF_FIELD)+"="+ss);
453 0 : }
454 0 : else if( fieldType(parm,RF_FIELD,TpInt,TpArrayInt))
455 : {
456 0 : parm.get(RF_FIELD,sel_fieldid);
457 0 : String ss;
458 0 : for( uInt i=0; i<sel_fieldid.nelements(); i++)
459 0 : addString(ss,String::toString(sel_fieldid(i)),",");
460 0 : addString(desc_str,String(RF_FIELD)+"="+ss);
461 0 : sel_fieldid -= (Int)indexingBase();
462 0 : }
463 : // parse input arguments: Scan Number(s)
464 0 : if( fieldType(parm,RF_SCAN,TpInt,TpArrayInt))
465 : {
466 0 : parm.get(RF_SCAN,sel_scannumber);
467 0 : String ss;
468 0 : for( uInt i=0; i<sel_scannumber.nelements(); i++)
469 0 : addString(ss,String::toString(sel_scannumber(i)),",");
470 0 : addString(desc_str,String(RF_SCAN)+"="+ss);
471 0 : sel_scannumber -= (Int)indexingBase();
472 0 : }
473 : // parse input arguments: Scan intent
474 0 : if ( fieldType(parm,RF_INTENT,TpInt,TpArrayInt))
475 : {
476 0 : parm.get(RF_INTENT,sel_stateid);
477 0 : String ss;
478 0 : for( uInt i=0; i<sel_stateid.nelements(); i++)
479 0 : addString(ss,String::toString(sel_stateid(i)),",");
480 0 : addString(desc_str,String(RF_INTENT)+"="+ss);
481 0 : sel_stateid -= (Int)indexingBase();
482 0 : }
483 : // parse input arguments: Array ID(s)
484 0 : if( fieldType(parm,RF_ARRAY,TpInt,TpArrayInt))
485 : {
486 0 : parm.get(RF_ARRAY,sel_arrayid);
487 0 : String ss;
488 0 : for( uInt i=0; i<sel_arrayid.nelements(); i++)
489 0 : addString(ss,String::toString(sel_arrayid(i)),",");
490 0 : addString(desc_str,String(RF_ARRAY)+"="+ss);
491 0 : sel_arrayid -= (Int)indexingBase();
492 0 : }
493 : // parse input arguments: Observation ID(s)
494 0 : if( fieldType(parm,RF_OBSERVATION,TpInt,TpArrayInt))
495 : {
496 0 : parm.get(RF_OBSERVATION,sel_observation);
497 0 : String ss;
498 0 : for( uInt i=0; i<sel_observation.nelements(); i++)
499 0 : addString(ss,String::toString(sel_observation(i)),",");
500 0 : addString(desc_str,String(RF_OBSERVATION)+"="+ss);
501 0 : sel_observation -= (Int)indexingBase();
502 0 : }
503 : // parse input: specific time ranges
504 0 : Array<Double> rng;
505 0 : Matrix<Double> timerng;
506 0 : if( parseTimes(rng,parm,RF_TIMERANGE))
507 : {
508 0 : if( !reformRange(timerng,rng))
509 0 : os << "Illegal \"" << RF_TIMERANGE << "\" array\n" << LogIO::EXCEPTION;
510 0 : sel_timerng = timerng*(Double)(24*3600);
511 :
512 0 : String s(String(RF_TIMERANGE) + "("); // + String::toString(timerng.ncolumn()));
513 :
514 : Bool del;
515 0 : Double *ptimes = rng.getStorage(del);
516 0 : for (unsigned i = 0; i < rng.nelements(); i++) {
517 0 : s += String::toString(ptimes[i]);
518 0 : if (i < rng.nelements()-1) {
519 0 : s += ' ';
520 : }
521 : }
522 0 : rng.putStorage(ptimes, del);
523 0 : s += ")";
524 :
525 0 : addString(desc_str, s);
526 0 : }
527 :
528 : // parse input: specific UV ranges
529 0 : Array<Double> uvrng;
530 0 : Matrix<Double> uvrange;
531 0 : if( parseTimes(uvrng,parm,RF_UVRANGE))
532 : {
533 0 : if( !reformRange(uvrange,uvrng))
534 0 : os<<"Illegal \""<<RF_UVRANGE<<"\" array\n"<<LogIO::EXCEPTION;
535 0 : sel_uvrange = uvrange;
536 0 : addString(desc_str,String(RF_UVRANGE)+"("+String::toString(uvrange.ncolumn())+")");
537 : }
538 :
539 : // parse input arguments: ANT specified by string ID
540 0 : LogicalVector sel_ant(num(ANT),false);
541 0 : if( fieldType(parm,RF_ANT,TpString,TpArrayString))
542 : {
543 0 : Vector<String> sant( parm.asArrayString(RF_ANT)) ;
544 0 : sant.apply(stringUpper);
545 0 : const Vector<String> &names( chunk.antNames()) ;
546 0 : for( uInt i=0; i<sant.nelements(); i++)
547 : {
548 : uInt iant;
549 0 : if( !find( iant,sant(i),names))
550 0 : os<<"Illegal antenna ID "<<sant(i)<<endl<<LogIO::EXCEPTION;
551 0 : sel_ant(iant)=true;
552 : }
553 0 : }
554 : // else ANT specified directly by indexes
555 0 : else if( fieldType(parm,RF_ANT,TpInt,TpArrayInt))
556 : {
557 0 : Vector<Int> sant = parm.asArrayInt(RF_ANT);
558 0 : for( uInt i=0; i<sant.nelements(); i++)
559 0 : sel_ant( sant(i) - (Int)indexingBase()) = true;
560 0 : }
561 0 : if( sum(sel_ant))
562 : {
563 0 : String sant;
564 0 : for( uInt i=0; i<num(ANT); i++)
565 0 : if( sel_ant(i))
566 0 : addString(sant, chunk.antNames()(i),",");
567 0 : addString(desc_str, String(RF_ANT)+"="+sant);
568 0 : }
569 :
570 : // parse input: baselines as "X-Y"
571 0 : sel_ifr = LogicalVector(num(IFR),false);
572 0 : String ifrdesc;
573 0 : const Vector<String> &names( chunk.antNames()) ;
574 0 : if( fieldType(parm, RF_BASELINE, TpString, TpArrayString))
575 : {
576 0 : Vector<String> ss(parm.asArrayString(RF_BASELINE));
577 0 : ss.apply(stringUpper);
578 0 : for( uInt i=0; i<ss.nelements(); i++)
579 : {
580 : uInt ant1,ant2;
581 0 : String ants[2];
582 0 : Int res = split(ss(i),ants,2,"-");
583 0 : Bool wild1 = (ants[0]=="*" || ants[0]=="") ; // is it a wildcard instead of ID?
584 0 : Bool wild2 = (ants[1]=="*" || ants[1]=="") ;
585 0 : if( res<2 || ( wild1 && wild2))
586 0 : os<<"Illegal baseline specification "<<ss(i)<<endl<<LogIO::EXCEPTION;
587 0 : Bool val1 = wild1 || find(ant1,ants[0],names);
588 0 : Bool val2 = wild2 || find(ant2,ants[1],names);
589 : // if both antenna IDs are valid, use them
590 0 : if( val1 && val2)
591 : {
592 0 : if( wild1)
593 : {
594 0 : addString(ifrdesc,ants[1]+"-*",",");
595 0 : for( uInt a=0; a<num(ANT); a++)
596 0 : sel_ifr( chunk.antToIfr(a,ant2)) = true;
597 : }
598 0 : else if( wild2)
599 : {
600 0 : addString(ifrdesc,ants[0]+"-*",",");
601 0 : for( uInt a=0; a<num(ANT); a++)
602 0 : sel_ifr( chunk.antToIfr(ant1,a)) = true;
603 : }
604 : else
605 : {
606 0 : addString(ifrdesc,ants[0]+"-"+ants[1],",");
607 0 : sel_ifr( chunk.antToIfr(ant1,ant2)) = true;
608 : }
609 0 : }
610 : else // try to interpret them as numbers instead
611 : {
612 0 : if( sscanf(ss(i).chars(),"%d-%d",&ant1,&ant2)<2 ||
613 0 : ant1>=num(ANT) || ant2>=num(ANT))
614 0 : os<<"Illegal baseline specification "<<ss(i)<<endl<<LogIO::EXCEPTION;
615 0 : sel_ifr( chunk.antToIfr(ant1-(Int)indexingBase(),ant2-(Int)indexingBase())) = true;
616 0 : addString(ifrdesc,ss(i),",");
617 : }
618 0 : }
619 0 : }
620 : // parse input: baselines as [[x1,y1],[x2,y2],... etc.
621 0 : else if( fieldType(parm,RF_BASELINE,TpInt,TpArrayInt))
622 : {
623 0 : Matrix<Int> ant;
624 0 : if( parseRange(ant,parm,RF_BASELINE))
625 : {
626 0 : ant -= (Int)indexingBase();
627 0 : for( uInt i=0; i<ant.ncolumn(); i++)
628 : {
629 0 : if( ant(0,i)==-1)
630 : {
631 0 : if( ant(1,i)==-1)
632 0 : os<<"Illegal baseline specification [-1,-1]"<<LogIO::EXCEPTION<<endl;
633 0 : for( uInt a=0; a<num(ANT); a++)
634 0 : sel_ifr( chunk.antToIfr(a,ant(1,i))) = true;
635 0 : addString(ifrdesc,names(ant(1,i))+"-*",",");
636 : }
637 0 : else if( ant(1,i)==-1)
638 : {
639 0 : for( uInt a=0; a<num(ANT); a++)
640 0 : sel_ifr( chunk.antToIfr(ant(0,i),a)) = true;
641 0 : addString(ifrdesc,names(ant(0,i))+"-*",",");
642 : }
643 : else
644 : {
645 0 : unsigned indx = chunk.antToIfr(ant(0,i),ant(1,i));
646 :
647 0 : assert( indx < sel_ifr.nelements() );
648 :
649 0 : sel_ifr(indx) = true;
650 0 : addString(ifrdesc,names(ant(0,i))+"-"+names(ant(1,i)),",");
651 : }
652 : }
653 : }
654 0 : }
655 :
656 0 : if( sum(sel_ifr))
657 : {
658 0 : String ss;
659 0 : addString(desc_str,String(RF_BASELINE)+"="+ifrdesc);
660 0 : }
661 : else // no IFRs were specified
662 : {
663 0 : if( sum(sel_ant)) // antennas specified? flag only their baselines
664 : {
665 0 : for( uInt a1=0; a1<num(ANT); a1++)
666 0 : if( sel_ant(a1))
667 0 : for( uInt a2=0; a2<num(ANT); a2++)
668 0 : sel_ifr(chunk.antToIfr(a1,a2)) = true;
669 : }
670 : else // no antennas either? flag everything
671 0 : sel_ifr.resize();
672 : }
673 :
674 : // parse input: feeds as [[x1,y1],[x2,y2],... etc.
675 0 : if( fieldType(parm,RF_FEED,TpInt,TpArrayInt))
676 : {
677 0 : sel_feed = LogicalVector(num(FEEDCORR),false);
678 0 : String feeddesc;
679 0 : Matrix<Int> feed;
680 0 : if( parseRange(feed,parm,RF_FEED))
681 : {
682 0 : feed -= (Int)indexingBase();
683 0 : for( uInt i=0; i<feed.ncolumn(); i++)
684 : {
685 0 : if( feed(0,i)==-1)
686 : {
687 0 : if( feed(1,i)==-1)
688 0 : os<<"Illegal feed specification [-1,-1]"<<LogIO::EXCEPTION<<endl;
689 0 : for( uInt a=0; a<num(FEED); a++)
690 0 : sel_feed( chunk.antToIfr(a,feed(1,i))) = true;
691 0 : addString(feeddesc,names(feed(1,i))+"-*",",");
692 : }
693 0 : else if( feed(1,i)==-1)
694 : {
695 0 : for( uInt a=0; a<num(FEED); a++)
696 0 : sel_feed( chunk.antToIfr(feed(0,i),a)) = true;
697 0 : addString(feeddesc,names(feed(0,i))+"-*",",");
698 : }
699 : else
700 : {
701 0 : sel_feed(chunk.antToIfr(feed(0,i),feed(1,i))) = true;
702 0 : addString(feeddesc,names(feed(0,i))+"-"+names(feed(1,i)),",");
703 : }
704 : }
705 : }
706 0 : }
707 :
708 : // now, all selection-related arguments are accounted for.
709 : // set flag if some subset has been selected
710 0 : Bool have_subset = ( desc_str.length());
711 :
712 0 : if (have_subset)
713 0 : desc_str+=";";
714 : // unflag specified?
715 0 : unflag = (fieldType(parm, RF_UNFLAG,TpBool) && parm.asBool(RF_UNFLAG));
716 : //addString(desc_str,unflag?RF_UNFLAG:"flag");
717 :
718 0 : ac = new MSAntennaColumns(chunk.measSet().antenna());
719 0 : diameters = ac->dishDiameter().getColumn();
720 :
721 0 : shadow = fieldType(parm, RF_SHADOW, TpBool) && parm.asBool(RF_SHADOW);
722 0 : if (shadow) {
723 0 : diameter = parm.asDouble(RF_DIAMETER);
724 : }
725 :
726 0 : elevation = fieldType(parm, RF_ELEVATION, TpBool) && parm.asBool(RF_ELEVATION);
727 0 : if (elevation) {
728 0 : lowerlimit = parm.asDouble(RF_LOWERLIMIT);
729 0 : upperlimit = parm.asDouble(RF_UPPERLIMIT);
730 : }
731 :
732 : // now, scan arguments for what to flag within the selection
733 : // specific times (specified by center times)
734 0 : Vector<Double> ctimes;
735 0 : Double timedelta = 10;
736 0 : if (parseTimes(ctimes,parm,RF_CENTERTIME))
737 : {
738 0 : ctimes *= (Double)(24*3600);
739 0 : String ss (String(RF_CENTERTIME)+"("+String::toString(ctimes.nelements())+")");
740 0 : Vector<Double> dt;
741 0 : if (parseTimes(dt,parm,RF_TIMEDELTA,true))
742 : {
743 0 : timedelta = dt(0);
744 0 : sprintf(s,",dtime=%.1fs",timedelta);
745 0 : ss += s;
746 : }
747 0 : addString(desc_str,ss);
748 0 : uInt n = ctimes.nelements();
749 0 : sel_time.resize(2,n);
750 0 : sel_time.row(0) = ctimes - timedelta;
751 0 : sel_time.row(1) = ctimes + timedelta;
752 0 : }
753 :
754 : // flag autocorrelations too?
755 0 : sel_autocorr = (fieldType(parm,RF_AUTOCORR,TpBool) && parm.asBool(RF_AUTOCORR));
756 0 : if (sel_autocorr)
757 0 : addString(desc_str,RF_AUTOCORR);
758 :
759 : // parse input: quack mode (for VLA)
760 0 : if (isFieldSet(parm, RF_QUACK)) {
761 :
762 : //quack_si = 30.0; // scan interval
763 0 : quack_dt = 10.0; // dt to flag at start of scan
764 :
765 : // are specific values given?
766 0 : Vector<Double> qt;
767 0 : if (parseTimes(qt, parm, RF_QUACK, true)) {
768 :
769 0 : if (qt.nelements() > 3) {
770 0 : os << RF_QUACK << " must be specified as T, <scaninterval> or [scaninterval,dt]"
771 0 : << endl << LogIO::EXCEPTION;
772 : }
773 :
774 0 : quack_si = qt(0);
775 0 : if (qt.nelements() > 2) {
776 0 : quack_dt = qt(1);
777 : }
778 :
779 0 : assert(parm.isDefined(RF_QUACKMODE));
780 0 : assert(parm.isDefined(RF_QUACKINC));
781 :
782 0 : quack_mode = parm.asString(RF_QUACKMODE);
783 0 : quack_increment = parm.asBool(RF_QUACKINC);
784 :
785 : }
786 0 : sprintf(s, "%s=%ds", // %s=%s; %s=%s",
787 0 : RF_QUACK, (Int)quack_si);
788 : //RF_QUACKMODE, quack_mode,
789 : //RF_QUACKINC, quack_increment ? "true" : "false");
790 0 : addString(desc_str, s);
791 : // quack_si /= (24*3600);
792 : // quack_dt /= (24*3600);
793 0 : }
794 : else {
795 0 : quack_si = 0;
796 : }
797 :
798 : // flag a specific range or clip outside of range?
799 :
800 0 : if (isValidRecord(parm,RF_CLIP)) {
801 0 : parseClipField(parm.asRecord(RF_CLIP),true);
802 : }
803 :
804 0 : if (isValidRecord(parm,RF_FLAGRANGE))
805 0 : parseClipField(parm.asRecord(RF_FLAGRANGE),false);
806 :
807 : // add to description strings, if something was parsed
808 0 : if (sel_clip.nelements())
809 : {
810 0 : addClipInfoDesc(sel_clip);
811 0 : sel_clip_active.resize(sel_clip.nelements());
812 : }
813 0 : if (sel_clip_row.nelements())
814 0 : addClipInfoDesc(sel_clip_row);
815 :
816 : // if nothing has been specified to flag, flag everything within selection
817 0 : flag_everything =
818 0 : (quack_si == 0 &&
819 0 : !sel_time.nelements() &&
820 0 : !sel_autocorr &&
821 0 : !sel_clip.nelements() &&
822 0 : !sel_clip_row.nelements() &&
823 0 : !shadow &&
824 0 : !elevation &&
825 0 : !sel_stateid.nelements());
826 : //!elevation);
827 : /*
828 : if (flag_everything)
829 : {
830 : if (!have_subset && !unflag)
831 : os<<"FLAG ALL requested, but no MS subset specified.\n"
832 : "Refusing to flag the whole measurement set!\n"<<LogIO::EXCEPTION;
833 : // addString(desc_str,"all");
834 : }
835 : */
836 : // cout<<"Selector: "<<desc_str<<endl;
837 :
838 0 : }
839 :
840 :
841 0 : RFASelector::~RFASelector ()
842 : {
843 0 : delete ac;
844 :
845 0 : for( uInt i=0; i<sel_clip.nelements(); i++ )
846 0 : if( sel_clip[i].mapper )
847 0 : delete sel_clip[i].mapper;
848 0 : }
849 :
850 0 : void RFASelector::startData(bool verbose)
851 : {
852 0 : RFAFlagCubeBase::startData(verbose);
853 :
854 0 : String flagstring = unflag?String("unflag"):String("flag");
855 :
856 0 : os << LogIO::DEBUGGING << "Data flagged/unflagged : " << desc_str << " " << flagstring;
857 0 : if (flag_everything) os << " all" ;
858 0 : os << LogIO::POST;
859 :
860 0 : Bool have_subset = ( desc_str.length() );
861 :
862 0 : if( flag_everything && !shadow)
863 : {
864 : /* jmlarsen: This does not seem useful nor necessary */
865 :
866 : if (false) if( !have_subset && !unflag)
867 : os<<"FLAG ALL requested, but no MS subset specified.\n"
868 : "Refusing to flag the whole measurement set!\n"<<LogIO::EXCEPTION;
869 : }
870 :
871 0 : return;
872 0 : }
873 :
874 : // -----------------------------------------------------------------------
875 : // newChunk
876 : // At each new chunk, figure out what goes where
877 : // -----------------------------------------------------------------------
878 0 : Bool RFASelector::newChunk (Int &maxmem)
879 : {
880 : // check correlations and figure out the active correlations mask
881 0 : Vector<Int> corrtype;
882 0 : chunk.visIter().corrType(corrtype);
883 0 : corrmask = 0;
884 0 : if( sel_corr.nelements() )
885 : {
886 0 : corrmask = chunk.getCorrMask(sel_corr);
887 0 : if( !corrmask )
888 : {
889 0 : if(verbose2) os<<"No matching correlations in this chunk\n"<<LogIO::POST;
890 0 : return active=false;
891 : }
892 : }
893 : else // no correlations specified so flag everything
894 0 : corrmask = chunk.fullCorrMask();
895 :
896 : // check field IDs and spectral window IDs
897 : uInt dum;
898 0 : if( sel_spwid.nelements() && !find(dum,chunk.visBuf().spectralWindow(),sel_spwid) )
899 : {
900 0 : if(verbose2) os<<"Spectral window does not match in this chunk\n"<<LogIO::POST;
901 0 : return active=false;
902 : }
903 0 : if( sel_fieldid.nelements() && !find(dum,chunk.visIter().fieldId(),sel_fieldid) )
904 : {
905 0 : if(verbose2) os<<"Field ID does not match in this chunk\n"<<LogIO::POST;
906 0 : return active=false;
907 : }
908 0 : if( sel_fieldnames.nelements() && !find(dum,chunk.visIter().fieldName(),sel_fieldnames) )
909 : {
910 0 : if(verbose2) os<<"Field name does not match in this chunk\n"<<LogIO::POST;
911 0 : return active=false;
912 : }
913 0 : if( sel_arrayid.nelements() && !find(dum,chunk.visIter().arrayId(),sel_arrayid) )
914 : {
915 0 : if(verbose2) os<<"Array ID does not match in this chunk\n"<<LogIO::POST;
916 0 : return active=false;
917 : }
918 0 : if( sel_observation.nelements() && !find(dum,chunk.visBuf().observationId()[0],sel_observation) )
919 : {
920 0 : if(verbose2) os<<"Observation ID does not match in this chunk\n"<<LogIO::POST;
921 0 : return active=false;
922 : }
923 : //Vector<Int> tempstateid(0);
924 : //tempstateid = chunk.visIter().stateId(tempstateid);
925 : //if( tempstateid.nelements() && sel_stateid.nelements() && !find(dum,tempstateid[0],sel_stateid) )
926 : //{
927 : // if(verbose2) os<<"State ID does not match in this chunk\n"<<LogIO::POST;
928 : // return active=false;
929 : //}
930 : //
931 : /*
932 : Vector<Int> temp(0);
933 : temp = chunk.visIter().scan(temp);
934 : if( temp.nelements() && sel_scannumber.nelements() && !find(dum,temp[0],sel_scannumber) )
935 : {
936 : os<<"Scan Number does not match the FIRST scan number in this chunk\n"<<LogIO::POST;
937 : return active=false;
938 : }
939 : */
940 :
941 : // Get the times at the beginning and end of this scan.
942 :
943 0 : const Vector<Int> &scans(chunk.visBuf().scan());
944 :
945 0 : Int s0 = scans(0);
946 :
947 0 : if (!allEQ(scans, s0) && quack_si > 0) {
948 0 : os << "RFASelector: VisBuffer has given us different scans (in this chunk)."
949 0 : << LogIO::EXCEPTION;
950 : }
951 :
952 :
953 : //cout << "scan range is = " <<
954 : //MVTime( scan_start/C::day).string( MVTime::DMY,7) << " - " <<
955 : //MVTime( scan_end /C::day).string( MVTime::DMY,7) << endl;
956 :
957 :
958 : // figure out active channels (i.e. within specified segments)
959 0 : flagchan.resize();
960 0 : if( sel_freq.ncolumn() || sel_chan.ncolumn() )
961 : {
962 0 : flagchan = LogicalVector(num(CHAN),false);
963 0 : const Vector<Double> & fq( chunk.frequency() );
964 0 : for( uInt i=0; i<sel_freq.ncolumn(); i++ )
965 0 : flagchan = flagchan || ( fq >= sel_freq(0,i) && fq <= sel_freq(1,i) );
966 0 : Vector<Int> ch( num(CHAN) );
967 0 : indgen(ch);
968 0 : Matrix<Int> schan = sel_chan;
969 0 : schan( sel_chan<0 ) += (Int)num(CHAN);
970 0 : for( uInt i=0; i<sel_chan.ncolumn(); i++ )
971 : {
972 0 : flagchan = flagchan || ( ch >= schan(0,i) && ch <= schan(1,i) );
973 : }
974 0 : if( !sum(flagchan) )
975 : {
976 0 : if(verbose2) os<<"No matching frequencies/channels in this chunk\n"<<LogIO::POST;
977 0 : return active=false;
978 : }
979 0 : if( allEQ(flagchan,true) )
980 0 : flagchan.resize(); // null array = all true
981 0 : }
982 : // init all clipping mappers, and check their correlation masks
983 0 : if( sel_clip.nelements() )
984 : {
985 : // see which mappers are active for this chunk, and accumulate their
986 : // masks in clip_corrmask
987 0 : RFlagWord clip_corrmask=0;
988 0 : for( uInt i=0; i<sel_clip.nelements(); i++ )
989 : {
990 0 : RFlagWord mask = sel_clip[i].mapper->corrMask(chunk.visIter());
991 0 : sel_clip_active(i) = (mask!=0);
992 0 : clip_corrmask |= mask;
993 : }
994 0 : sum_sel_clip_active = sum(sel_clip_active);
995 : // If no explicit correlations were selected by the user,
996 : // then use clip_corrmask as the overall mask
997 0 : if( !sel_corr.nelements() )
998 : {
999 0 : corrmask = clip_corrmask;
1000 0 : if( !corrmask )
1001 : {
1002 0 : if (verbose2) {
1003 0 : os << "No matching correlations in this chunk\n" << LogIO::POST;
1004 : }
1005 0 : return active=false;
1006 : }
1007 : }
1008 : }
1009 : // reserve a minimum of memory for ourselves
1010 0 : maxmem -= 1;
1011 : // init flagging cube and off we go...
1012 0 : RFAFlagCubeBase::newChunk(maxmem);
1013 : // see if full row is being flagged, i.e. no subset of channels was selected,
1014 : // and no explicit correlations (or all correlations).
1015 0 : select_fullrow = (!flagchan.nelements() &&
1016 0 : (!sel_corr.nelements() || corrmask==chunk.fullCorrMask()) );
1017 0 : return active=true;
1018 0 : }
1019 :
1020 : // -----------------------------------------------------------------------
1021 : // processRow
1022 : // Raises/clears flags for a single row, depending on current selection
1023 : // -----------------------------------------------------------------------
1024 0 : void RFASelector::processRow(uInt ifr,uInt it)
1025 : {
1026 0 : if(dbg2) cout << ifr << ",";
1027 : // does the selection include whole rows?
1028 0 : if (select_fullrow) {
1029 :
1030 0 : unflag ? flag.clearRowFlag(ifr,it) : flag.setRowFlag(ifr,it);
1031 : // apply this to the entire row...
1032 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
1033 0 : unflag ? flag.clearFlag(ich,ifr) : flag.setFlag(ich,ifr);
1034 :
1035 : }
1036 : else {
1037 : // else apply data flags to selection
1038 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
1039 0 : if( !flagchan.nelements() || flagchan(ich) )
1040 0 : unflag ? flag.clearFlag(ich,ifr) : flag.setFlag(ich,ifr);
1041 : }
1042 0 : }
1043 :
1044 : // -----------------------------------------------------------------------
1045 : // Processes 1 time slot in the MS
1046 : // There is one MS row per baseline
1047 : //
1048 : // -----------------------------------------------------------------------
1049 0 : RFA::IterMode RFASelector::iterTime (uInt it)
1050 : {
1051 0 : RFAFlagCubeBase::iterTime(it);
1052 :
1053 : // setup each correlation clip mapper
1054 0 : for (uInt i=0; i<sel_clip.nelements(); i++)
1055 0 : if (sel_clip_active(i))
1056 0 : sel_clip[i].mapper->setVisBuffer(chunk.visBuf());
1057 :
1058 : // extract time
1059 0 : const Vector<Double> ×( chunk.visBuf().time() );
1060 0 : const Vector<Double> &dtimes( chunk.visBuf().timeInterval() );
1061 0 : Double t0 = times(0);
1062 0 : Double dt0 = dtimes(0);
1063 0 : if( !allEQ(times,t0) )
1064 0 : os << "RFASelector: VisBuffer has given us different times." << LogIO::EXCEPTION;
1065 :
1066 : // extract scan number
1067 0 : const Vector<Int> &scans( chunk.visBuf().scan() );
1068 :
1069 0 : Int s0 = scans(0);
1070 :
1071 0 : if (!allEQ(scans,s0) && quack_si > 0) {
1072 0 : os << "RFASelector: crash&burn, VisBuffer's given us different scans."
1073 0 : << LogIO::EXCEPTION;
1074 : }
1075 :
1076 : // is current scan number within selected list of scan numbers?
1077 0 : bool scanselect = true;
1078 0 : if (sel_scannumber.nelements()) {
1079 0 : bool sel = false;
1080 0 : for (uInt i = 0;
1081 0 : i < sel_scannumber.nelements();
1082 : i++) {
1083 0 : if( sel_scannumber(i) == s0 ) {
1084 0 : sel = true;
1085 : }
1086 : }
1087 :
1088 0 : if( ! sel ) scanselect = false;
1089 : //if( ! sel ) return RFA::CONT;
1090 : }
1091 :
1092 : // flag if within specific timeslots
1093 0 : bool timeselect = true;
1094 0 : if (sel_timerng.ncolumn()) {
1095 0 : timeselect = false;
1096 0 : if( anyEQ(sel_timerng.row(0) <= t0 && sel_timerng.row(1) >= t0,
1097 0 : true) ) {
1098 0 : timeselect = true;
1099 : }
1100 : }
1101 :
1102 0 : if ( ! (timeselect && scanselect) ) {
1103 0 : return RFA::CONT;
1104 : }
1105 :
1106 0 : const Vector<Int> &ifrs( chunk.ifrNums() );
1107 0 : const Vector<Int> &feeds( chunk.feedNums() );
1108 0 : const Vector<casacore::RigidVector<casacore::Double, 3> >&uvw( chunk.visBuf().uvw() );
1109 : // Vector<Vector<Double> > &uvw=NULL;//( chunk.visIter.uvw(uvw) );
1110 : //chunk.visIter().uvw(uvw);
1111 0 : Double uvdist=0.0;
1112 :
1113 0 : if( ifrs.nelements() != feeds.nelements() ||
1114 0 : ifrs.nelements() != uvw.nelements() )
1115 0 : cout << "RFASelector::iterTime -> nelements mismatch " << endl;
1116 :
1117 : // do we flag the whole selection?
1118 0 : bool flagall = flag_everything;
1119 :
1120 0 : if (!flagall) {
1121 :
1122 0 : if (elevation) {
1123 0 : const Vector<MDirection> &azimuth_elevation = chunk.visIter().azel(t0);
1124 :
1125 0 : for (uInt i = 0; i < ifrs.nelements(); i++) {
1126 :
1127 : unsigned a1, a2;
1128 0 : chunk.ifrToAnt(a1, a2, chunk.ifrNum(i));
1129 :
1130 0 : Bool inrange = false;
1131 0 : uvdist = sqrt( uvw(i)(0)*uvw(i)(0) + uvw(i)(1)*uvw(i)(1) );
1132 0 : for (uInt j = 0; j < sel_uvrange.ncolumn(); j++)
1133 0 : if (uvdist >= sel_uvrange(0, j) && uvdist <= sel_uvrange(1, j))
1134 0 : inrange |= true;
1135 :
1136 0 : if( (!sel_ifr.nelements() || sel_ifr(ifrs(i))) &&
1137 0 : (!sel_feed.nelements() || sel_feed(feeds(i))) &&
1138 0 : (!sel_uvrange.nelements() || inrange ) ) {
1139 :
1140 : // double antenna_elevation = azimuth_elevation[i].getAngle("deg").getValue()[1];
1141 0 : double antenna1_elevation = azimuth_elevation[a1].getAngle("deg").getValue()[1];
1142 0 : double antenna2_elevation = azimuth_elevation[a2].getAngle("deg").getValue()[1];
1143 :
1144 0 : if ( antenna1_elevation < lowerlimit ||
1145 0 : antenna2_elevation < lowerlimit ||
1146 0 : antenna1_elevation > upperlimit ||
1147 0 : antenna2_elevation > upperlimit ) {
1148 0 : processRow(ifrs(i), it);
1149 : }
1150 : }
1151 : }
1152 :
1153 0 : }
1154 :
1155 0 : if (shadow) {
1156 :
1157 : /*
1158 : 1st loop: Figure out which antennas are shadowed
1159 : 2nd loop: Flag all data (incl self correlations)
1160 : involving shadowed antennas
1161 : */
1162 :
1163 0 : std::vector<bool> shadowed(diameters.nelements(), false);
1164 :
1165 0 : for (uInt i = 0; i < ifrs.nelements(); i++) {
1166 :
1167 : unsigned a1, a2;
1168 0 : chunk.ifrToAnt(a1, a2, chunk.ifrNum(i));
1169 :
1170 0 : if (a1 != a2) { /* Antennas don't shadow themselves. */
1171 : double d1, d2;
1172 0 : if (diameter < 0) {
1173 0 : d1 = diameters(a1);
1174 0 : d2 = diameters(a2);
1175 : }
1176 : else {
1177 0 : d1 = diameter;
1178 0 : d2 = diameter;
1179 : }
1180 :
1181 : Double uvdist2 =
1182 0 : uvw(i)(0) * uvw(i)(0) +
1183 0 : uvw(i)(1) * uvw(i)(1);
1184 :
1185 : /* The relevant threshold distance for shadowing is
1186 : (d1+d2)/2 */
1187 0 : if (uvdist2 < (d1+d2)*(d1+d2)/4.0) {
1188 :
1189 : if (0) cerr << "antenna is shadowed " << a1 << "-" << a2 << ": " <<
1190 : "(u, v, w) = (" <<
1191 : uvw(i)(0) << ", " <<
1192 : uvw(i)(1) << ", " <<
1193 : uvw(i)(2) << ")" << endl;
1194 :
1195 0 : if (uvw(i)(2) > 0) {
1196 0 : shadowed[a1] = true;
1197 : }
1198 : else {
1199 0 : shadowed[a2] = true;
1200 : }
1201 : }
1202 : }
1203 : }
1204 :
1205 : if (0) {
1206 : std::copy(shadowed.begin(),
1207 : shadowed.end(),
1208 : std::ostream_iterator<bool>(std::cout, "] [" ));
1209 : cout << endl;
1210 : }
1211 :
1212 0 : for (uInt i = 0; i < ifrs.nelements(); i++) {
1213 :
1214 : unsigned a1, a2;
1215 0 : chunk.ifrToAnt(a1, a2, chunk.ifrNum(i));
1216 :
1217 0 : Bool inrange = false;
1218 0 : uvdist = sqrt( uvw(i)(0)*uvw(i)(0) + uvw(i)(1)*uvw(i)(1) );
1219 0 : for (uInt j = 0; j < sel_uvrange.ncolumn(); j++)
1220 0 : if (uvdist >= sel_uvrange(0, j) && uvdist <= sel_uvrange(1, j))
1221 0 : inrange |= true;
1222 :
1223 0 : if( (!sel_ifr.nelements() || sel_ifr(ifrs(i))) &&
1224 0 : (!sel_feed.nelements() || sel_feed(feeds(i))) &&
1225 0 : (!sel_uvrange.nelements() || inrange ) ) {
1226 :
1227 0 : if (shadowed[a1] || shadowed[a2]) {
1228 0 : processRow(ifrs(i),it);
1229 : }
1230 : }
1231 : }
1232 0 : } /* end if shadow */
1233 :
1234 : // flag autocorrelations
1235 0 : if (sel_autocorr) {
1236 0 : for (uInt i=0; i < ifrs.nelements(); i++) {
1237 0 : Bool inrange=false;
1238 0 : uvdist = sqrt( uvw(i)(0)*uvw(i)(0) + uvw(i)(1)*uvw(i)(1) );
1239 0 : for( uInt j=0; j<sel_uvrange.ncolumn(); j++)
1240 0 : if( uvdist >= sel_uvrange(0,j) && uvdist <= sel_uvrange(1,j) ) inrange |= true;
1241 : // if( inrange ) cout << "x selected : " << i << " : " << uvdist << endl;
1242 0 : if ((!sel_ifr.nelements() || sel_ifr(ifrs(i))) &&
1243 0 : (!sel_feed.nelements() || sel_feed(feeds(i))) &&
1244 0 : (!sel_uvrange.nelements() || inrange ))
1245 : {
1246 : uInt a1,a2;
1247 0 : chunk.ifrToAnt(a1,a2,ifrs(i));
1248 0 : if( a1==a2 )
1249 0 : processRow(ifrs(i),it);
1250 : }
1251 : }
1252 : }
1253 :
1254 : // flag if quacked
1255 0 : if (quack_si > 0) {
1256 :
1257 : double scan_start;
1258 : double scan_end;
1259 :
1260 0 : if (quack_increment) {
1261 0 : scan_start = chunk.get_scan_start_unflagged(s0);
1262 0 : scan_end = chunk.get_scan_end_unflagged(s0);
1263 : // returns negative if there's no unflagged data
1264 : }
1265 : else {
1266 0 : scan_start = chunk.get_scan_start(s0);
1267 0 : scan_end = chunk.get_scan_end(s0);
1268 : }
1269 : //cout << "Start time for scan : " << MVTime( scan_start/C::day).string( MVTime::DMY,7) ;
1270 : //cout << " ::: iterTime for " << MVTime( t0/C::day).string( MVTime::DMY,7) << " and scan " << s0;
1271 :
1272 0 : if (quack_mode == "beg") {
1273 0 : if (scan_start > 0 &&
1274 0 : t0 <= (scan_start + quack_si)) flagall = true;
1275 : }
1276 0 : else if (quack_mode == "endb") {
1277 0 : if (scan_end > 0 &&
1278 0 : t0 >= (scan_end - quack_si)) flagall = true;
1279 : }
1280 0 : else if (quack_mode == "end") {
1281 0 : if (scan_end > 0 &&
1282 0 : t0 < (scan_end - quack_si)) flagall = true;
1283 : }
1284 0 : else if (quack_mode == "tail") {
1285 0 : if (scan_start > 0 &&
1286 0 : t0 > (scan_start + quack_si)) flagall = true;
1287 : }
1288 : else {
1289 0 : throw(AipsError("Illegal quack mode '" + quack_mode + "'"));
1290 : }
1291 :
1292 : }
1293 :
1294 : // flag for specific row-based clipping expressions
1295 0 : if (sel_clip_row.nelements()) {
1296 :
1297 0 : for( uInt i=0; i < ifrs.nelements(); i++ ) {// loop over rows
1298 0 : for( uInt j=0; j < sel_clip_row.nelements(); j++ ) {
1299 :
1300 0 : Float vmin = sel_clip_row[j].vmin;
1301 0 : Float vmax = sel_clip_row[j].vmax;
1302 0 : Float val = sel_clip_row[j].mapper->mapValue(i) - sel_clip_row[j].offset;
1303 0 : if( (sel_clip_row[j].clip && ( val<vmin || val>vmax ) ) ||
1304 0 : (!sel_clip_row[j].clip && val>=vmin && val<=vmax ) )
1305 0 : processRow(ifrs(i),it);
1306 : }
1307 : }
1308 : }
1309 :
1310 0 : bool clip_based_on_tsys = false;
1311 0 : if (clip_based_on_tsys) {
1312 : /*
1313 : for each row:
1314 :
1315 : find antenna1 and antenna2
1316 : lookup temperature for
1317 : (antenna1,spwid,time) and
1318 : (antenna2,spwid,time) in the SYSCAL subtable
1319 :
1320 : */
1321 0 : cout << "Get sysCal table" << endl;
1322 : // note, this is an optional subtable
1323 :
1324 0 : const MSSysCal syscal(chunk.measSet().sysCal());
1325 :
1326 0 : ScalarColumn<uInt> sc_antenna_id(syscal, "ANTENNA_ID");
1327 0 : ScalarColumn<uInt> sc_spwid (syscal, "SPECTRAL_WINDOW_ID");
1328 0 : ScalarColumn<Double> sc_time (syscal, "TIME");
1329 0 : ArrayColumn<Float> sc_tsys (syscal, "TSYS");
1330 :
1331 0 : unsigned spwid = chunk.visBuf().spectralWindow();
1332 :
1333 0 : for (uInt i = 0; i < ifrs.nelements(); i++) {
1334 :
1335 : unsigned a1, a2;
1336 0 : chunk.ifrToAnt(a1, a2, chunk.ifrNum(i));
1337 :
1338 : //for each SYSCAL table row
1339 : // if antenna1 matches or antenna2 matches and
1340 : // spwid matches and time matches
1341 : // if (tsys out of allowed range)
1342 : // flag
1343 :
1344 0 : unsigned nrows_syscal = sc_antenna_id.getColumn().nelements();
1345 0 : cout << nrows_syscal << " rows in SYSCAL table" << endl;
1346 0 : for (unsigned row = 0; row < nrows_syscal; row++) {
1347 0 : cout << "a1, a2 = " << a1 << ", " << a2 << " ?= " << sc_antenna_id(row) << endl;
1348 0 : cout << "time = " << t0 << " ?= " << sc_time(row) << endl;
1349 0 : cout << "spwid = " << spwid << " ?= " << sc_spwid(row) << endl;
1350 0 : cout << "tsys = " << sc_tsys(row) << endl;
1351 :
1352 0 : if ((a1 == sc_antenna_id(row) || a2 == sc_antenna_id(row)) &&
1353 0 : spwid == sc_spwid(row) &&
1354 0 : t0-dt0/2 < sc_time(row) && sc_time(row) < t0+dt0/2) {
1355 0 : cout << " MATCH!" << endl;
1356 : }
1357 : }
1358 : }
1359 0 : }
1360 : // scan intent based flagging (check match per antenna/baseline)
1361 0 : if(sel_stateid.nelements()) {
1362 0 : const Vector<Int> &stateids( chunk.visBuf().stateId() );
1363 : //for (uInt i = 0; i < ifrs.nelements(); i++) {
1364 : // for (uInt j = 0; j < sel_stateid.nelements(); j++) {
1365 : // if( sel_stateid(j) == stateid0 )
1366 0 : for (uInt i = 0; i < stateids.nelements(); i++) {
1367 :
1368 0 : Bool inrange=false;
1369 0 : uvdist = sqrt( uvw(i)(0)*uvw(i)(0) + uvw(i)(1)*uvw(i)(1) );
1370 0 : for( uInt j=0; j<sel_uvrange.ncolumn(); j++)
1371 0 : if( uvdist >= sel_uvrange(0,j) && uvdist <= sel_uvrange(1,j) ) inrange |= true;
1372 :
1373 0 : for (uInt j = 0; j < sel_stateid.nelements(); j++) {
1374 0 : if( stateids(i) == sel_stateid(j) &&
1375 0 : ifrs.nelements()==stateids.nelements() &&
1376 0 : (!sel_ifr.nelements()||sel_ifr(ifrs(i))) &&
1377 0 : (!sel_feed.nelements() || sel_feed(feeds(i))) &&
1378 0 : (!sel_uvrange.nelements() || inrange ) )
1379 0 : processRow(ifrs(i),it);
1380 : }
1381 : }
1382 : }
1383 : }
1384 :
1385 : // flag whole selection, if still needed
1386 0 : if (flagall)
1387 : {
1388 0 : for (uInt i=0; i<ifrs.nelements(); i++) // loop over rows
1389 : {
1390 0 : Bool inrange=false;
1391 0 : uvdist = sqrt( uvw(i)(0)*uvw(i)(0) + uvw(i)(1)*uvw(i)(1) );
1392 0 : for( uInt j=0; j<sel_uvrange.ncolumn(); j++)
1393 0 : if( uvdist >= sel_uvrange(0,j) && uvdist <= sel_uvrange(1,j) ) inrange |= true;
1394 0 : if( (!sel_ifr.nelements() || sel_ifr(ifrs(i))) &&
1395 0 : (!sel_feed.nelements() || sel_feed(feeds(i))) &&
1396 0 : (!sel_uvrange.nelements() || inrange ) )
1397 0 : processRow(ifrs(i),it);
1398 : }
1399 : }
1400 :
1401 0 : return RFA::CONT;
1402 : }
1403 :
1404 : // -----------------------------------------------------------------------
1405 : // iterRow
1406 : // -----------------------------------------------------------------------
1407 0 : RFA::IterMode RFASelector::iterRow (uInt ir)
1408 : {
1409 0 : uInt ifr = chunk.ifrNum(ir);
1410 :
1411 0 : if( sel_clip.nelements() && sum_sel_clip_active )
1412 : {
1413 : // apply data flags
1414 0 : for( uInt j=0; j<sel_clip.nelements(); j++ ) {
1415 :
1416 0 : if (sel_clip[j].channel_average) {
1417 : // Compute average
1418 0 : Float average = 0;
1419 0 : unsigned n_unflagged = 0;
1420 :
1421 0 : for (uInt ich=0; ich < num(CHAN); ich++ ) {
1422 :
1423 : // if active channel, and not flagged...
1424 0 : if ((!flagchan.nelements() || flagchan(ich))
1425 0 : &&
1426 0 : !flag.preFlagged(ich, ifr)) {
1427 :
1428 :
1429 0 : average += sel_clip[j].mapper->mapValue(ich, ir);
1430 0 : n_unflagged += 1;
1431 : }
1432 : }
1433 :
1434 :
1435 : // Decide whether to flag row (or active channels), based on average
1436 : //cout << "number of unflagged channels = " << n_unflagged << endl;
1437 0 : if (n_unflagged > 0) {
1438 0 : average = average / n_unflagged;
1439 :
1440 : //cout << " average = " << average << endl;
1441 :
1442 0 : Float vmin = sel_clip[j].vmin;
1443 0 : Float vmax = sel_clip[j].vmax;
1444 0 : Float val = average;
1445 :
1446 0 : if( ( sel_clip[j].clip && ( val < vmin || val > vmax ) ) ||
1447 0 : (!sel_clip[j].clip && vmin <= val && val <= vmax ) )
1448 :
1449 : /* No channel selections, flag all channels.
1450 :
1451 : jmlarsen: Unfortunately, clearRowFlag and setRowFlag
1452 : don't seem to work, therefore don't do like this.
1453 : (This could probably be brought to work by
1454 : fixing RFFlagCube::setMSFlags() to use the flagrow.)
1455 :
1456 : cout << " flag entire row " << endl;
1457 : unflag ?
1458 : flag.clearRowFlag(ifr, it) :
1459 : flag.setRowFlag(ifr, it);
1460 : */
1461 :
1462 0 : for (uInt ich = 0; ich < num(CHAN); ich++) {
1463 0 : if (!flagchan.nelements() || flagchan(ich)) {
1464 0 : unflag ?
1465 0 : flag.clearFlag(ich, ifr) :
1466 0 : flag.setFlag(ich, ifr);
1467 : }
1468 : }
1469 : }
1470 : else {
1471 : /* If all channels were already flagged, don't flag/unflag anything. */
1472 : }
1473 : }
1474 : else {
1475 : /* No averaging */
1476 0 : for( uInt ich=0; ich<num(CHAN); ich++ ) {
1477 0 : if( !flagchan.nelements() || flagchan(ich) ) {
1478 0 : Float vmin = sel_clip[j].vmin;
1479 0 : Float vmax = sel_clip[j].vmax;
1480 0 : Float val = sel_clip[j].mapper->mapValue(ich,ir);
1481 :
1482 : // jagonzal: Added ISO isnan check to catch extremely large values (CAS-3355)
1483 0 : if( ( sel_clip[j].clip && (val<vmin || val>vmax || std::isnan(val)) ) ||
1484 0 : (!sel_clip[j].clip && val>=vmin && val<=vmax ) )
1485 0 : unflag ? flag.clearFlag(ich,ifr) : flag.setFlag(ich,ifr);
1486 : }
1487 : }
1488 : }
1489 : }
1490 : }
1491 :
1492 0 : return RFA::CONT;
1493 : }
1494 :
1495 : /* Flush flags after each time stamp, if in kiss mode */
1496 0 : void RFASelector::endRows(uInt it)
1497 : {
1498 0 : if (only_selector) {
1499 0 : flag.advance(it);
1500 0 : flag.setMSFlags(it);
1501 : }
1502 0 : }
1503 :
1504 0 : void RFASelector::iterFlag(uInt it)
1505 : {
1506 0 : if (!only_selector) {
1507 0 : RFAFlagCubeBase::iterFlag(it);
1508 : }
1509 0 : }
1510 :
1511 0 : String RFASelector::getDesc ()
1512 : {
1513 0 : return desc_str+" "+RFAFlagCubeBase::getDesc();
1514 : }
1515 :
1516 0 : const RecordInterface & RFASelector::getDefaults ()
1517 : {
1518 0 : static Record rec;
1519 : // create record description on first entry
1520 0 : if( !rec.nfields() )
1521 : {
1522 0 : rec = RFAFlagCubeBase::getDefaults();
1523 0 : rec.removeField(RF_FIGNORE); // fignore is meaningless
1524 0 : rec.define(RF_NAME,"Selector");
1525 0 : rec.define(RF_SPWID,false);
1526 0 : rec.define(RF_FIELD,false);
1527 0 : rec.define(RF_FREQS,false);
1528 0 : rec.define(RF_CHANS,false);
1529 0 : rec.define(RF_CORR,false);
1530 0 : rec.define(RF_ANT,false);
1531 0 : rec.define(RF_BASELINE,false);
1532 0 : rec.define(RF_TIMERANGE,false);
1533 0 : rec.define(RF_AUTOCORR,false);
1534 0 : rec.define(RF_CENTERTIME,false);
1535 0 : rec.define(RF_TIMEDELTA,10.0);
1536 0 : rec.define(RF_QUACK,false);
1537 0 : rec.define(RF_CLIP,false);
1538 0 : rec.define(RF_CHANAVG, false);
1539 0 : rec.define(RF_FLAGRANGE,false);
1540 0 : rec.define(RF_UNFLAG,false);
1541 0 : rec.define(RF_SHADOW,false);
1542 0 : rec.define(RF_ELEVATION, false);
1543 0 : rec.define(RF_SCAN,false);
1544 0 : rec.define(RF_INTENT,false);
1545 0 : rec.define(RF_ARRAY,false);
1546 0 : rec.define(RF_FEED,false);
1547 0 : rec.define(RF_UVRANGE,false);
1548 0 : rec.define(RF_COLUMN,false);
1549 0 : rec.define(RF_DIAMETER, false);
1550 0 : rec.define(RF_LOWERLIMIT, false);
1551 0 : rec.define(RF_UPPERLIMIT, false);
1552 :
1553 0 : rec.setComment(RF_SPWID,"Restrict flagging to specific spectral windows (integers)");
1554 0 : rec.setComment(RF_FIELD,"Restrict flagging to specific field IDs or field names (integers/strings)");
1555 0 : rec.setComment(RF_FREQS,"Restrict flagging to specific frequency ranges (2,N array of doubles:MHz)");
1556 0 : rec.setComment(RF_CHANS,"Restrict flagging to specific channels (2,N array of integers)");
1557 0 : rec.setComment(RF_CORR,"Restrict flagging to specific correlations (array of strings)");
1558 0 : rec.setComment(RF_ANT,"Restrict flagging to specific antennas (array of strings/integers)");
1559 0 : rec.setComment(RF_BASELINE,"Restrict flagging to specific baselines (array of strings, e.g., 'A1-A2','A1-*', or 2,N array of integers [[A1,A2],[B1,B2],...])");
1560 0 : rec.setComment(RF_TIMERANGE,"Restrict flagging to specific time ranges (2,N array of strings or MJDs");
1561 0 : rec.setComment(RF_AUTOCORR,"Flag autocorrelations (F/T)");
1562 0 : rec.setComment(RF_CENTERTIME,"Flag specific timeslots (array of strings or MJDs)");
1563 0 : rec.setComment(RF_TIMEDELTA,String("Time delta for ")+RF_CENTERTIME+", in seconds");
1564 0 : rec.setComment(RF_QUACK,"Use [SI,DT] for VLA quack-flagging");
1565 0 : rec.setComment(RF_CLIP,"Flag outside a specific range of values");
1566 0 : rec.setComment(RF_CHANAVG, "Flag based on channel average");
1567 0 : rec.setComment(RF_FLAGRANGE,"Flag inside a specific range of values");
1568 0 : rec.setComment(RF_UNFLAG,"If T, specified flags are CLEARED");
1569 0 : rec.setComment(RF_SHADOW, "If T, flag shadowed antennas");
1570 0 : rec.setComment(RF_ELEVATION, "If T, flag based on elevation");
1571 0 : rec.setComment(RF_SCAN,"Restrict flagging to specific scans (integers)");
1572 0 : rec.setComment(RF_INTENT,"Restrict flagging to specific scan intent -corresponding state IDs (integers)");
1573 0 : rec.setComment(RF_ARRAY,"Restrict flagging to specific array ids (integers)");
1574 0 : rec.setComment(RF_FEED,"Restrict flagging to specific feeds (2,N array of integers)");
1575 0 : rec.setComment(RF_UVRANGE,"Restrict flagging to specific uv-distance ranges in meters (2,N array of doubles )");
1576 0 : rec.setComment(RF_COLUMN,"Data column to clip on.");
1577 0 : rec.setComment(RF_DIAMETER, "Effective diameter to use. If negative, the true antenna diameters are used");
1578 0 : rec.setComment(RF_LOWERLIMIT, "Limiting elevation");
1579 0 : rec.setComment(RF_UPPERLIMIT, "Limiting elevation");
1580 : }
1581 0 : return rec;
1582 : }
1583 :
1584 : } //# NAMESPACE CASA - END
1585 :
|
