Line data Source code
1 : //# RFAFlagExaminer.cc: this defines RFAFlagExaminer
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: aips2-request@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 <flagging/Flagging/RFAFlagExaminer.h>
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 <map>
39 : #include <sstream>
40 : #include <cassert>
41 :
42 : namespace casa { //# NAMESPACE CASA - BEGIN
43 :
44 : const casacore::Bool dbg3 = false;
45 :
46 : // -----------------------------------------------------------------------
47 : // RFAFlagExaminer constructor
48 : // -----------------------------------------------------------------------
49 0 : RFAFlagExaminer::RFAFlagExaminer ( RFChunkStats &ch,const casacore::RecordInterface &parm ) :
50 0 : RFASelector(ch, parm)//,RFDataMapper(parm.asArrayString(RF_EXPR),parm.asString(RF_COLUMN))
51 : {
52 : if(dbg3) std::cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
53 : //desc_str = casacore::String("flagexaminer");
54 : if(dbg3) std::cout<<"FlagExaminer constructor "<<std::endl;
55 :
56 0 : totalflags = accumTotalFlags = 0;
57 0 : totalcount = accumTotalCount = 0;
58 0 : totalrowflags = accumTotalRowFlags = 0;
59 0 : totalrowcount = accumTotalRowCount = 0;
60 : //parseParm(parm);
61 :
62 0 : os = casacore::LogIO(casacore::LogOrigin("RFAFlagExaminer", "RFAFlagExaminer", WHERE));
63 :
64 0 : accumflags.clear();
65 0 : accumtotal.clear();
66 :
67 :
68 : // Handle in-row selections, the following is a
69 : // copy-paste from RFASelector2
70 :
71 : char s[256];
72 : // parse input arguments: channels
73 0 : if( parseRange(sel_chan,parm,RF_CHANS))
74 : {
75 0 : casacore::String sch;
76 0 : for( casacore::uInt i=0; i<sel_chan.ncolumn(); i++)
77 : {
78 0 : sprintf(s,"%d:%d",sel_chan(0,i),sel_chan(1,i));
79 0 : addString(sch,s,",");
80 : }
81 0 : addString(desc_str, casacore::String(RF_CHANS) + "=" +sch);
82 0 : sel_chan(sel_chan>=0) += -(casacore::Int)indexingBase();
83 : }
84 :
85 : // parse input arguments: correlations
86 0 : if( fieldType(parm,RF_CORR,casacore::TpString,casacore::TpArrayString))
87 : {
88 0 : casacore::String ss;
89 0 : casacore::Vector<casacore::String> scorr( parm.asArrayString(RF_CORR)) ;
90 0 : sel_corr.resize( scorr.nelements()) ;
91 0 : for( casacore::uInt i=0; i<scorr.nelements(); i++)
92 : {
93 0 : sel_corr(i) = casacore::Stokes::type( scorr(i)) ;
94 0 : if( sel_corr(i) == casacore::Stokes::Undefined)
95 0 : os<<"Illegal correlation "<<scorr(i)<<std::endl<<casacore::LogIO::EXCEPTION;
96 0 : addString(ss,scorr(i),",");
97 : }
98 0 : addString(desc_str,casacore::String(RF_CORR)+"="+ss);
99 : }
100 0 : }
101 :
102 :
103 0 : RFAFlagExaminer::~RFAFlagExaminer()
104 : {
105 : if(dbg3) std::cout << "FlagExaminer destructor " << std::endl;
106 0 : }
107 :
108 0 : casacore::Bool RFAFlagExaminer::newChunk(casacore::Int &maxmem)
109 : {
110 : /* For efficiency reasons, use arrays to collect
111 : histogram data for in-row selections
112 : */
113 0 : accumflags_channel = std::vector<casacore::uInt64>(chunk.num(CHAN), 0);
114 0 : accumtotal_channel = std::vector<casacore::uInt64>(chunk.num(CHAN), 0);
115 0 : accumflags_correlation = std::vector<casacore::uInt64>(chunk.num(CORR), 0);
116 0 : accumtotal_correlation = std::vector<casacore::uInt64>(chunk.num(CORR), 0);
117 :
118 0 : return RFASelector::newChunk(maxmem);
119 : }
120 :
121 :
122 0 : void RFAFlagExaminer::initialize()
123 : {
124 : if(dbg3) std::cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
125 :
126 0 : totalflags = 0;
127 0 : totalcount = 0;
128 0 : totalrowflags = 0;
129 0 : totalrowcount = 0;
130 0 : inTotalFlags =
131 0 : inTotalCount =
132 0 : inTotalRowCount =
133 0 : outTotalFlags =
134 0 : outTotalCount =
135 0 : outTotalRowCount =
136 0 : outTotalRowFlags = 0;
137 0 : }
138 :
139 : // Is not called if this is the only agent
140 0 : void RFAFlagExaminer::finalize()
141 : {
142 : //cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
143 :
144 0 : return;
145 : }
146 : // -----------------------------------------------------------------------
147 : // processRow
148 : // Raises/clears flags for a single row, depending on current selection
149 : // -----------------------------------------------------------------------
150 0 : void RFAFlagExaminer::processRow(casacore::uInt, casacore::uInt)
151 : {
152 : // called often... if(dbg3) std::cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
153 0 : return;
154 : }
155 :
156 0 : void RFAFlagExaminer::startFlag (bool verbose)
157 : {
158 : if(dbg3) std::cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
159 :
160 0 : totalflags = 0;
161 0 : totalcount = 0;
162 0 : totalrowflags = 0;
163 0 : totalrowcount = 0;
164 :
165 0 : inTotalFlags =
166 0 : inTotalCount =
167 0 : inTotalRowCount =
168 0 : outTotalFlags =
169 0 : outTotalCount =
170 0 : outTotalRowCount =
171 0 : outTotalRowFlags = 0;
172 :
173 0 : RFAFlagCubeBase::startFlag(verbose);
174 :
175 0 : return;
176 : }
177 :
178 0 : void RFAFlagExaminer::initializeIter (casacore::uInt)
179 : {
180 : if(dbg3) std::cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
181 :
182 0 : for(unsigned ii=0;
183 0 : ii<chunk.visBuf().flagRow().nelements();
184 : ii++)
185 0 : if (chunk.visBuf().flagRow()(ii)) {
186 0 : inTotalRowFlags++;
187 : }
188 :
189 0 : inTotalRowCount += chunk.visBuf().flagRow().nelements();
190 :
191 0 : for(casacore::Int ii=0;
192 0 : ii<chunk.visBuf().flag().shape()(0);
193 : ii++)
194 0 : for(casacore::Int jj=0;
195 0 : jj<chunk.visBuf().flag().shape()(1);
196 : jj++)
197 0 : if (chunk.visBuf().flag()(ii,jj)) inTotalFlags++;
198 0 : }
199 :
200 : // Is not called if this is the only agent
201 0 : void RFAFlagExaminer::finalizeIter (casacore::uInt)
202 : {
203 : //cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
204 :
205 0 : outTotalRowCount += chunk.visBuf().flagRow().nelements();
206 :
207 0 : for (unsigned ii = 0;
208 0 : ii < chunk.visBuf().flagRow().nelements();
209 : ii++)
210 0 : if (chunk.visBuf().flagRow()(ii)) {
211 0 : outTotalRowFlags++;
212 : }
213 :
214 0 : for (casacore::Int ii=0;
215 0 : ii<chunk.visBuf().flag().shape()(0);
216 : ii++) {
217 :
218 0 : outTotalCount += chunk.visBuf().flag().shape()(1);
219 :
220 0 : for (casacore::Int jj=0;
221 0 : jj<chunk.visBuf().flag().shape()(1);
222 : jj++) {
223 0 : if (chunk.visBuf().flag()(ii,jj)) outTotalFlags++;
224 : }
225 : }
226 0 : }
227 :
228 : //also need to call RFFlagCube::setMSFlags, which
229 : // updates some statistics void RFAFlagExaminer::endRows(casacore::uInt it)
230 :
231 : // it: time index
232 0 : void RFAFlagExaminer::iterFlag(casacore::uInt it)
233 : {
234 : if(dbg3) std::cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
235 :
236 : // Set the flags and count them up.
237 0 : RFASelector::iterFlag(it);
238 :
239 : // count if within specific timeslots
240 0 : const casacore::Vector<casacore::Double> ×( chunk.visBuf().time() );
241 0 : casacore::Double t0 = times(it);
242 :
243 0 : casacore::Bool within_time_slot = true;
244 :
245 0 : if (sel_time.ncolumn()) {
246 :
247 0 : if( anyEQ(sel_timerng.row(0) <= t0 &&
248 0 : sel_timerng.row(1) >= t0, true) )
249 0 : within_time_slot = true;
250 0 : else within_time_slot = false;
251 : }
252 :
253 0 : if (within_time_slot) {
254 :
255 : // More counting and fill up final display variables.
256 0 : const casacore::Vector<casacore::Int> &ifrs( chunk.ifrNums() );
257 0 : const casacore::Vector<casacore::Int> &feeds( chunk.feedNums() );
258 0 : const casacore::Vector<casacore::RigidVector<casacore::Double, 3> >&uvw( chunk.visBuf().uvw() );
259 :
260 0 : unsigned spw = chunk.visBuf().spectralWindow();
261 0 : unsigned field = chunk.visBuf().fieldId();
262 0 : const casacore::Vector<casacore::Int> &antenna1( chunk.visBuf().antenna1() );
263 0 : const casacore::Vector<casacore::Int> &antenna2( chunk.visBuf().antenna2() );
264 0 : const casacore::Vector<casacore::Int> &scan ( chunk.visBuf().scan() );
265 0 : const casacore::Vector<casacore::Int> &observation ( chunk.visBuf().observationId() );
266 0 : casacore::Int array ( chunk.visBuf().arrayId() );
267 :
268 0 : const casacore::Vector<casacore::String> &antenna_names( chunk.antNames()) ;
269 :
270 : // casacore::Vector<Vector<casacore::Double> > &uvw=NULL;//( chunk.visIter.uvw(uvw) );
271 : //chunk.visIter().uvw(uvw);
272 0 : casacore::Double uvdist=0.0;
273 :
274 : // loop over rows
275 0 : for (casacore::uInt i=0; i < ifrs.nelements(); i++) {
276 0 : casacore::Bool inrange=false;
277 :
278 0 : uvdist = sqrt( uvw(i)(0)*uvw(i)(0) + uvw(i)(1)*uvw(i)(1) );
279 :
280 0 : for( casacore::uInt j=0; j<sel_uvrange.ncolumn(); j++) {
281 0 : if( uvdist >= sel_uvrange(0, j) &&
282 0 : uvdist <= sel_uvrange(1, j) )
283 :
284 0 : inrange |= true;
285 : }
286 :
287 0 : if( (!sel_ifr.nelements() || sel_ifr(ifrs(i))) &&
288 0 : (!sel_feed.nelements() || sel_feed(feeds(i))) &&
289 0 : (!sel_uvrange.nelements() || inrange ) )
290 : {
291 : // Operate on the chosen row.
292 : // Collect counts.
293 :
294 : //cout << "selected row for " << ifrs(i) << "," << it << std::endl;
295 :
296 0 : if(chunk.nfIfrTime(ifrs(i),it) == chunk.num(CORR)*chunk.num(CHAN))
297 0 : totalrowflags++;
298 :
299 0 : totalrowcount++;
300 :
301 0 : casacore::uInt64 f = chunk.nfIfrTime(ifrs(i), it);
302 0 : casacore::uInt64 c = chunk.num(CORR) * chunk.num(CHAN);
303 :
304 : // need nfIfrTimeCorr
305 : // need nfIfrTimeChan
306 :
307 0 : totalflags += f;
308 0 : totalcount += c;
309 :
310 : /* Update histograms */
311 :
312 : /* histogram baseline */
313 : {
314 : string baseline =
315 0 : antenna_names(antenna1(i)) + "&&" +
316 0 : antenna_names(antenna2(i));
317 :
318 0 : accumflags["baseline"][baseline] += f;
319 0 : accumtotal["baseline"][baseline] += c;
320 : }
321 :
322 : /* histogram antenna */
323 : {
324 : /* Careful here, update the counts for both
325 : antenna1 and antenna2, unless they are the same.
326 : */
327 0 : accumflags["antenna"][antenna_names(antenna1(i))] += f;
328 0 : accumtotal["antenna"][antenna_names(antenna1(i))] += c;
329 0 : if (antenna1(i) != antenna2(i)) {
330 0 : accumflags["antenna"][antenna_names(antenna2(i))] += f;
331 0 : accumtotal["antenna"][antenna_names(antenna2(i))] += c;
332 : }
333 : }
334 :
335 : /* histogram spw */
336 : {
337 0 : std::stringstream spw_string;
338 0 : spw_string << spw;
339 0 : accumflags["spw"][spw_string.str()] += f;
340 0 : accumtotal["spw"][spw_string.str()] += c;
341 : }
342 :
343 : /* histogram fieldID */
344 : {
345 0 : std::stringstream fieldID_string;
346 0 : fieldID_string << field;
347 0 : accumflags["field"][fieldID_string.str()] += f;
348 0 : accumtotal["field"][fieldID_string.str()] += c;
349 : }
350 : /* histogram scan */
351 : {
352 0 : std::stringstream scan_string;
353 0 : scan_string << scan(i);
354 0 : accumflags["scan"][scan_string.str()] += f;
355 0 : accumtotal["scan"][scan_string.str()] += c;
356 : }
357 : /* histogram observation */
358 : {
359 0 : std::stringstream observation_string;
360 0 : observation_string << observation(i);
361 0 : accumflags["observation"][observation_string.str()] += f;
362 0 : accumtotal["observation"][observation_string.str()] += c;
363 : }
364 : /* histogram array */
365 : {
366 0 : std::stringstream array_string;
367 0 : array_string << array;
368 0 : accumflags["array"][array_string.str()] += f;
369 0 : accumtotal["array"][array_string.str()] += c;
370 : }
371 : }
372 : }
373 : }
374 :
375 0 : return;
376 : }
377 :
378 : /* Update histogram for channel and correlation
379 : This cannot happen in iterFlag(), which is called once per chunk,
380 : because the "flag" cursor needs to be updated once per row.
381 : */
382 0 : RFA::IterMode RFAFlagExaminer::iterRow(casacore::uInt irow)
383 : {
384 0 : unsigned ifr = chunk.ifrNum(irow);
385 :
386 0 : for( casacore::uInt ich = 0;
387 0 : ich < chunk.num(CHAN);
388 : ich++ ) {
389 :
390 0 : if( !flagchan.nelements() || flagchan(ich) ) {
391 :
392 0 : RFlagWord corrs = flag.getFlag(ich, ifr);
393 0 : unsigned n_flags = 0;
394 :
395 0 : for (casacore::uInt icorr = 0; icorr < chunk.num(CORR); icorr++) {
396 0 : if (corrs & 1) {
397 0 : accumflags_correlation[icorr] += 1;
398 0 : n_flags += 1;
399 : }
400 0 : accumtotal_correlation[icorr] += 1;
401 0 : corrs >>= 1;
402 : }
403 :
404 0 : accumflags_channel[ich] += n_flags;
405 0 : accumtotal_channel[ich] += chunk.num(CORR);
406 : }
407 : }
408 :
409 0 : return RFA::CONT;
410 : }
411 :
412 0 : void RFAFlagExaminer::endChunk ()
413 : {
414 0 : RFASelector::endChunk();
415 :
416 0 : for (unsigned ich = 0; ich < chunk.num(CHAN); ich++) {
417 0 : if (accumtotal_channel[ich] > 0) {
418 0 : std::stringstream ss;
419 0 : ss << chunk.visIter().spectralWindow() << ":" << ich;
420 0 : accumflags["channel"][ss.str()] = accumflags_channel[ich];
421 0 : accumtotal["channel"][ss.str()] = accumtotal_channel[ich];
422 : }
423 : }
424 :
425 0 : for (unsigned icorr = 0; icorr < chunk.num(CORR); icorr++) {
426 0 : if (accumtotal_correlation[icorr] > 0) {
427 0 : std::stringstream ss;
428 0 : ss << chunk.visIter().spectralWindow() << ":" << icorr;
429 0 : accumflags["correlation"][ss.str()] = accumflags_correlation[icorr];
430 0 : accumtotal["correlation"][ss.str()] = accumtotal_correlation[icorr];
431 : }
432 : }
433 0 : }
434 :
435 :
436 0 : void RFAFlagExaminer::endFlag ()
437 : {
438 : if(dbg3) std::cout << __FILE__ << ":" << __func__ << "():" << __LINE__ << std::endl;
439 :
440 : char s[1024];
441 0 : sprintf(s,"Chunk %d (field %s, fieldID %d, spw %d)",
442 0 : chunk.nchunk(),
443 0 : chunk.visIter().fieldName().chars(),
444 0 : chunk.visIter().fieldId(),
445 0 : chunk.visIter().spectralWindow());
446 0 : os << "---------------------------------------------------------------------" << casacore::LogIO::POST;
447 0 : os<<s<<casacore::LogIO::POST;
448 :
449 0 : sprintf(s, "%s, %d channel%s, %d time slots, %d baselines, %d rows\n",
450 0 : chunk.getCorrString().chars(),
451 0 : chunk.num(CHAN),
452 0 : chunk.num(CHAN) == 1 ? "" : "s",
453 0 : chunk.num(TIME),
454 0 : chunk.num(IFR),
455 0 : chunk.num(ROW));
456 0 : os << s << casacore::LogIO::POST;
457 :
458 0 : os << "\n\n\nData Selection to examine : " << desc_str ;
459 0 : if(flag_everything) os << " all " ;
460 0 : os << casacore::LogIO::POST;
461 :
462 0 : casacore::Double ffrac=0.0,rffrac=0.0;
463 0 : if(totalcount) ffrac = totalflags*100.0/totalcount;
464 0 : if(totalrowcount) rffrac = totalrowflags*100.0/totalrowcount;
465 :
466 : // casacore::LogSink cannot handle uInt64...
467 0 : std::stringstream ss;
468 0 : ss << totalrowflags << " out of " << totalrowcount <<
469 0 : " (" << rffrac << "%) rows are flagged.";
470 0 : os << ss.str() << casacore::LogIO::POST;
471 :
472 0 : ss.str("");
473 :
474 0 : ss << totalflags << " out of " << totalcount <<
475 0 : " (" << ffrac << "%) data points are flagged.\n\n";
476 :
477 0 : os << ss.str() << casacore::LogIO::POST;
478 :
479 0 : os << "---------------------------------------------------------------------" << casacore::LogIO::POST;
480 :
481 0 : accumTotalFlags += totalflags;
482 0 : accumTotalCount += totalcount;
483 0 : accumTotalRowFlags += totalrowflags;
484 0 : accumTotalRowCount += totalrowcount;
485 :
486 0 : return;
487 : }
488 :
489 : /*
490 : Return results of this run over all chunks:
491 : How many flags were set / unset
492 : */
493 0 : casacore::Record RFAFlagExaminer::getResult()
494 : {
495 0 : casacore::Record r;
496 :
497 0 : r.define("flagged", (casacore::Double) accumTotalFlags);
498 0 : r.define("total" , (casacore::Double) accumTotalCount);
499 :
500 :
501 0 : for (std::map<string, std::map<string, casacore::uInt64> >::iterator j = accumtotal.begin();
502 0 : j != accumtotal.end();
503 0 : j++) {
504 : /* Note here: loop over the keys of accumtotal, not accumflags,
505 : because accumflags may not have all channel keys */
506 :
507 0 : casacore::Record prop;
508 0 : for (std::map<string, casacore::uInt64>::const_iterator i = j->second.begin();
509 0 : i != j->second.end();
510 0 : i++) {
511 :
512 0 : casacore::Record t;
513 :
514 0 : t.define("flagged", (casacore::Double) accumflags[j->first][i->first]);
515 0 : t.define("total", (casacore::Double) i->second);
516 :
517 0 : prop.defineRecord(i->first, t);
518 : }
519 :
520 0 : r.defineRecord(j->first, prop);
521 : }
522 :
523 0 : return r;
524 : }
525 :
526 :
527 : } //# NAMESPACE CASA - END
528 :
|
