Line data Source code
1 : //# RFATimeFreqCrop.cc: this defines RFATimeFreqCrop
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 :
28 : #include <flagging/Flagging/RFATimeFreqCrop.h>
29 :
30 : namespace casa { //# NAMESPACE CASA - BEGIN
31 :
32 : //#define baselinecnt ( (NumAnt)*((NumAnt)+1)/2 - ((NumAnt)-ant1[bs])*((NumAnt)-ant1[bs]+1)/2 + (ant2[bs] - ant1[bs]) )
33 : #define SELF(ant) ( (NumAnt)*((NumAnt)+1)/2 - ((NumAnt)-ant)*((NumAnt)-ant+1)/2 )
34 : #define MIN(a,b) ((a)<=(b) ? (a) : (b))
35 :
36 : //#define PLOT // to activate the mean and clean bandpass plots
37 : //#define UPLOT // to activate bandpass plots of each fit iteration
38 :
39 : //#define DOPLOT // to activate ds9 and bandpass-fit plots
40 :
41 : /* Constructor for 'RFATimeFreqCrop' */
42 0 : RFATimeFreqCrop :: RFATimeFreqCrop( RFChunkStats &ch,const casacore::RecordInterface &parm ):
43 : RFAFlagCubeBase(ch,parm) ,
44 : RFDataMapper(parm.asArrayString(RF_EXPR),parm.asString(RF_COLUMN)),
45 0 : vi(ch.visIter()),
46 0 : vb(ch.visBuf())
47 : {
48 0 : ANT_TOL = parm.asDouble("ant_cutoff");
49 0 : BASELN_TOL = parm.asDouble("baseline_cutoff");
50 0 : T_TOL = parm.asDouble("time_amp_cutoff");
51 0 : F_TOL = parm.asDouble("freq_amp_cutoff");
52 0 : FlagLevel = parm.asInt("flag_level");
53 0 : CorrChoice = parm.asInt("auto_cross");
54 0 : NumTime = parm.asInt("num_time");
55 0 : ShowPlots = parm.asBool("showplots");
56 0 : FreqLineFit = parm.asBool("freqlinefit");
57 0 : MaxNPieces = parm.asInt("maxnpieces");
58 0 : DryRun = parm.asBool("dryrun");
59 0 : Expr = parm.asArrayString(RF_EXPR);
60 0 : Column = parm.asString(RF_COLUMN);
61 0 : IgnorePreflags = parm.asBool(RF_FIGNORE);
62 : // cout << "Flagging on " << parm.asArrayString(RF_EXPR) << " for column : " << parm.asString(RF_COLUMN) << endl;
63 0 : StopAndExit=false;
64 : /*
65 : if(ShowPlots)
66 : {
67 : cmd = "xpaset -p ds9 plot new name flagger \n";
68 : system(cmd.data());
69 : }
70 : */
71 :
72 0 : }
73 :
74 :
75 : /* Sets default values to parameters */
76 0 : const casacore::RecordInterface & RFATimeFreqCrop::getDefaults ()
77 : {
78 0 : static casacore::Record rec;
79 0 : if( !rec.nfields() )
80 : {
81 0 : rec = RFAFlagCubeBase::getDefaults();
82 0 : rec.define(RF_NAME,"RFATimeFreqCrop");
83 0 : rec.define(RF_EXPR,"ABS I");
84 0 : rec.define(RF_COLUMN,"DATA");
85 0 : rec.define("ant_cutoff",0);
86 0 : rec.define("baseline_cutoff",0);
87 0 : rec.define("time_amp_cutoff",3);
88 0 : rec.define("freq_amp_cutoff",3);
89 0 : rec.define("flag_level",1);
90 0 : rec.define("auto_cross",1);
91 0 : rec.define("num_time",50);
92 0 : rec.define("column","DATA");
93 0 : rec.define("expr","ABS I");
94 0 : rec.define("fignore",false);
95 0 : rec.define("showplots",false);
96 0 : rec.define("freqlinefit",false);
97 0 : rec.define("maxnpieces",6);
98 0 : rec.define("dryrun",false);
99 : // rec.setcomment("ant_cutoff","Total autocorrelation amplitude threshold for a functional antenna");
100 : // rec.setcomment("time_amp_cutoff","Multiple/fraction of standard deviation, to set the threshold, while flagging across time");
101 : // rec.setcomment("freq_amp_cutoff","Multiple/fraction of standard deviation, to set the threshold, while flagging across frequency");
102 : // rec.setcomment("flag_level","Levels of Flagging");
103 : }
104 0 : return rec;
105 : }
106 :
107 :
108 :
109 : /* Called at the beginning of each chunk of data */
110 0 : casacore::Bool RFATimeFreqCrop :: newChunk (casacore::Int &i)
111 : {
112 0 : casacore::LogIO os(casacore::LogOrigin("tfcrop","newChunk","WHERE"));
113 :
114 0 : if(StopAndExit)
115 : {
116 : // cout << "newChunk :: NOT Working with data chunk : " << chunk.msName() << endl;
117 0 : return false;
118 : }
119 :
120 :
121 : // cout << "newChunk :: Working with data chunk : " << chunk.msName() << endl;
122 : // cout << "TimeSteps = " << num(TIME) << ", Baselines = " << num(IFR) << ", Chans = " << num(CHAN) << ", Polns = " << num(POLZN) << ", Ants = " << num(ANT) << endl;
123 : // cout << "Parameters : " << " Antenna_tol=" << ANT_TOL << ", Baseline_tol=" << BASELN_TOL << ", Time_tol=" << T_TOL << "sigma, Freq_tol=" << F_TOL << "sigma, FlagLevel=" << FlagLevel << ", Flag_corr=" << casacore::String(CorrChoice?"Cross":"Auto") << endl;
124 :
125 : /* Initialize NumT - number of timestamps to work with in one go */
126 0 : if(NumTime==0) NumTime=50;
127 0 : if(num(TIME) >= NumTime) NumT = NumTime;
128 0 : else NumT = num(TIME);
129 :
130 : /* Assume that all baselines are present in this chunk */
131 : // TODO : check this.
132 0 : NumAnt = num(ANT);
133 0 : NumB = ((NumAnt)*((NumAnt)-1)/2 + (NumAnt));
134 :
135 : /* Number of polarizations */
136 0 : NumP = num(POLZN);
137 : //UUU : FORCE it to be 1 for now.....
138 0 : NumP=1;
139 :
140 : /* Number of channels */
141 0 : NumC = num(CHAN);
142 :
143 : /* Polarizations/correlations */
144 :
145 0 : corrmask = RFDataMapper::corrMask(chunk.visIter());
146 0 : casacore::Vector<casacore::Int> corrlist(num(POLZN));
147 0 : for(casacore::uInt corr=0; corr<num(POLZN); corr++)
148 0 : corrlist[corr] = (casacore::Int) ( (corrmask >> corr) & 1 );
149 :
150 : /* Check that the above makes sense */
151 0 : if(NumC<1 || NumP<1 || NumB <1 || NumAnt<1 || NumT<1)
152 : {
153 0 : cout << "Invalid chunk shapes" << endl;
154 0 : return false;
155 : }
156 :
157 : // os << "Chunk=" << chunk.nchunk() << ", Field=" << chunk.visIter().fieldName() << ", FieldID=" << chunk.visIter().fieldId() << ", Spw=" << chunk.visIter().spectralWindow() << ", nTime=" << num(TIME) << " (" << NumT << " at a time), nBaseline=" << num(IFR) << ", nChan=" << num(CHAN) << ", nCorrs=" << num(POLZN) << " [" << chunk.getCorrString() << "]" << casacore::LogIO::POST;
158 : //cout << "Chunk=" << chunk.nchunk() << ", Field=" << chunk.visIter().fieldName() << ", FieldID=" << chunk.visIter().fieldId() << ", Spw=" << chunk.visIter().spectralWindow() << ", nTime=" << num(TIME) << " (" << NumT << " at a time), nBaseline=" << num(IFR) << ", nChan=" << num(CHAN) << ", nCorrs=" << num(POLZN) << " [" << chunk.getCorrString() << "]" << endl; // ". Flagging on " << Expr << endl; //" -> correlations : " << corrlist << endl;
159 : // cout << "Working with " << NumC << " x " << NumT << " subsets (nchan x ntime), "<< Expr << " on the " << Column << " column, and applying flags to correlations : " << corrlist << endl;
160 :
161 :
162 : /* UUU : What is this ? */
163 : // RFAFlagCubeBase::newChunk(i-=1);
164 :
165 : /* Allocate memory for one set of NumT timesteps. */
166 0 : AllocateMemory();
167 :
168 0 : return RFAFlagCubeBase::newChunk(i);
169 0 : }
170 :
171 :
172 : /* Called at the beginning of each PASS */
173 0 : void RFATimeFreqCrop :: startData (bool verbose)
174 : {
175 : // cout << "StartData - reset time-counter" << endl;
176 :
177 0 : RFAFlagCubeBase::startData(verbose);
178 :
179 0 : iterTimecnt=0; // running count of visbuffers gone by.
180 0 : timecnt=0;
181 :
182 : /// (chunk.visIter()).setRowBlocking(0);
183 :
184 0 : }
185 :
186 0 : void RFATimeFreqCrop :: AllocateMemory()
187 : {
188 :
189 : /* casacore::Cube to hold visibility amplitudes : POLZN x CHAN x (IFR*TIME) */
190 0 : visc.resize(NumP,NumC,NumB*NumT);
191 0 : visc=0;
192 :
193 : /* casacore::Cube to hold visibility flags : POLZN x CHAN x (IFR*TIME) */
194 0 : flagc.resize(NumP,NumC,NumB*NumT);
195 0 : flagc=true;
196 :
197 : /* casacore::Vector to hold Row Flags : (IFR*TIME) */
198 0 : rowflags.resize(NumB*NumT);
199 0 : rowflags=true;
200 :
201 : /* casacore::Vector to hold baseline flags - to prevent unnecessary computation */
202 0 : baselineflags.resize(NumB);
203 0 : baselineflags=false;
204 :
205 : /* casacore::Cube to hold MEAN bandpasses : POLZN x IFR x CHAN */
206 : /* casacore::Cube to hold CLEAN bandpasses : POLZN x IFR x CHAN */
207 0 : if(CorrChoice == 0)
208 : {
209 0 : meanBP.resize(NumP,NumAnt,NumC);
210 0 : cleanBP.resize(NumP,NumAnt,NumC);
211 : }
212 : else
213 : {
214 0 : meanBP.resize(NumP,NumB,NumC);
215 0 : cleanBP.resize(NumP,NumB,NumC);
216 : }
217 :
218 0 : matpos = meanBP.shape();
219 0 : meanBP=0.0;
220 0 : cleanBP=0.0;
221 :
222 : // cout << " BP Shape = " << matpos << endl;
223 :
224 : /* casacore::Cube to hold flags for the entire Chunk (channel subset, but all times) */
225 0 : chunkflags.resize(NumP,NumC,NumB*num(TIME));
226 0 : chunkflags=true;
227 :
228 :
229 : /* Temporary workspace vectors */
230 0 : tempBP.resize(NumC);tempTS.resize(NumT);
231 0 : flagBP.resize(NumC);flagTS.resize(NumT);
232 0 : fitBP.resize(NumC);fitTS.resize(NumT);
233 :
234 0 : tempBP=0;tempTS=0;flagBP=false;flagTS=false;fitBP=0;fitTS=0;
235 :
236 :
237 0 : }
238 :
239 :
240 :
241 : /* Called once for every TIMESTAMP - for each VisBuf */
242 0 : RFA::IterMode RFATimeFreqCrop :: iterTime (casacore::uInt itime)
243 : {
244 :
245 : // cout << "iterTime :: " << itime << endl;
246 : // RFAFlagCubeBase::iterTime(itime);
247 0 : RFDataMapper::setVisBuffer(vb);
248 0 : flag.advance(itime,true);
249 :
250 : // vv = &vb.visCube(); // extract a viscube - one timestamp - one VisBuf
251 : // vi.flag(ff); // extract the corresponding flags
252 :
253 0 : if(ant1.shape() != (vb.antenna1()).shape())
254 0 : ant1.resize((vb.antenna1()).shape());
255 0 : ant1 = vb.antenna1();
256 :
257 0 : if(ant2.shape() != (vb.antenna2()).shape())
258 0 : ant2.resize((vb.antenna2()).shape());
259 0 : ant2 = vb.antenna2();
260 : //const casacore::Vector<casacore::Int> &ifrs( chunk.ifrNums() );
261 :
262 0 : casacore::uInt nBaselinesInData = ant1.nelements();
263 : // cout << "ant1 nelements : " << nBaselinesInData << " timecnt : " << timecnt << " itertimecnt : " << iterTimecnt << endl;
264 :
265 : /* Polarizations/correlations */
266 0 : corrmask = RFDataMapper::corrMask(chunk.visIter());
267 0 : casacore::Vector<casacore::Int> corrlist(num(POLZN));
268 0 : for(casacore::uInt corr=0; corr<num(POLZN); corr++)
269 0 : corrlist[corr] = (casacore::Int) ( (corrmask >> corr) & 1 );
270 :
271 : // if(nBaselinesInData != num(IFR)) cout << "nbaselines is not consistent !" << endl;
272 :
273 : // Read in the data AND any existing flags into the flagCube - accumulate
274 : // timecnt : casacore::Time counter for each NumT subset
275 : // itime : The row counter for each chunk
276 : casacore::Bool tfl;
277 0 : for(casacore::uInt pl=0;pl<NumP;pl++)
278 : {
279 0 : for(casacore::uInt bs=0;bs<nBaselinesInData;bs++)
280 : {
281 0 : casacore::Int countflags=0, countpnts=0;
282 0 : casacore::uInt baselinecnt = BaselineIndex(bs,ant1[bs],ant2[bs]);
283 0 : AlwaysAssert( baselinecnt<NumB, casacore::AipsError );
284 : // Read in rowflag
285 0 : rowflags( (timecnt*NumB)+baselinecnt ) = flag.getRowFlag( chunk.ifrNum(bs), itime);
286 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
287 : {
288 : // read the data. mapvalue evaluates 'expr'.
289 0 : visc(pl,ch,(timecnt*NumB)+baselinecnt) = mapValue(ch,bs);
290 : // read existing flags
291 0 : tfl = chunk.npass() ? flag.anyFlagged(ch,chunk.ifrNum(bs)) : flag.preFlagged(ch,chunk.ifrNum(bs));
292 : // sync with rowflag
293 0 : if( rowflags( (timecnt*NumB)+baselinecnt ) ) tfl=true;
294 : // ignore previous flags....
295 0 : if(IgnorePreflags) tfl=false;
296 :
297 : // Fill in the NumT sized flag array
298 0 : flagc(pl,ch,(timecnt*NumB)+baselinecnt) = tfl; //flag.anyFlagged(ch,ifrs(bs));
299 :
300 : // Fill in the chunk sized flag array
301 0 : chunkflags(pl,ch,(itime*NumB)+baselinecnt) = tfl; //flag.anyFlagged(ch,ifrs(bs));
302 :
303 : // Counters
304 0 : countpnts++;
305 0 : if(tfl) countflags ++;
306 : }
307 : // if(countflags>0) cout << "Time : " << itime << " Preflags for baseline : " << bs << " (" << ant1(bs) << "," << ant2(bs) << ") : " << countflags << " out of " << countpnts << " " << corrlist << endl;
308 : }
309 : }
310 :
311 0 : timecnt++;
312 0 : iterTimecnt++; // running count of visbuffers going by.
313 :
314 :
315 : /* BEGIN TIME-FLAGGING ALGORITHM HERE */
316 :
317 : /* After accumulating NumT timestamps, start processing this block */
318 : /////if(iterTimecnt > 0 && (timecnt==NumT || iterTimecnt == (vi.nRowChunk()/NumB)))
319 0 : if(iterTimecnt > 0 && (timecnt==NumT || itime==(num(TIME)-1) ))
320 : {
321 : //ut << " timecnt : " << timecnt << " itime : " << itime << " iterTimecnt : " << iterTimecnt << " NumT : " << NumT << endl;
322 : // casacore::Int ctimes = timecnt;
323 0 : casacore::Int ctimes = NumT; // User-specified time-interval
324 0 : NumT = timecnt; // Available time-interval. Usually same as NumT - but could be less.
325 : //ut << " NumT going into all functions : " << NumT << endl;
326 :
327 0 : FlagZeros();
328 :
329 0 : RunTFCrop();
330 :
331 0 : if(ShowFlagPlots() == RFA::STOP)
332 0 : return RFA::STOP;
333 :
334 0 : ExtendFlags();
335 :
336 0 : FillChunkFlags();
337 :
338 : // reset NumT to the user-specified time-interval
339 0 : NumT = ctimes;
340 :
341 : // reset the NumT time counter !!!
342 0 : timecnt=0;
343 : }
344 : // flag.advance(itime,true);
345 0 : return RFA::CONT;
346 : ////RFAFlagCubeBase::iterTime(itime);
347 0 : }
348 :
349 :
350 : /* Called for each ROW - each baseline in a VisBuf */
351 : // Fill in the visibilities and flags here.
352 : // flag.advance() has to get called in iterTime, for iterRow to see the correct values.
353 0 : RFA::IterMode RFATimeFreqCrop :: iterRow (casacore::uInt irow )
354 : {
355 0 : AlwaysAssert( irow <= NumT*NumB , casacore::AipsError);
356 : /* DUMMY CALL */
357 0 : return RFAFlagCubeBase::iterRow(irow);
358 : }
359 :
360 :
361 :
362 : /* If any data points are exactly zero, make sure corresponding flags
363 : are set. For the baseline mapper - this is an indicator of which baselines
364 : are missing (RowFlags */
365 0 : void RFATimeFreqCrop :: FlagZeros()
366 : {
367 0 : casacore::Float temp=0;
368 : //casacore::Bool flg=false;
369 0 : baselineflags=false;
370 :
371 : /* Check if the data in all channels are filled with zeros.
372 : If so, set the flags to zero */
373 : /* Also, if rowflags are set, set flags to zero. */
374 : /* Also, if all chans and times are flagged for a baseline, set the baselineflag
375 : baselineflag is used internally to skip unnecessary baselines */
376 :
377 0 : for(casacore::uInt pl=0;pl<NumP;pl++)
378 : {
379 0 : for(casacore::uInt bs=0;bs<NumB;bs++)
380 : {
381 0 : casacore::Bool bflag=true; // default is flagged. If anything is unflagged, this will change to false
382 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
383 : {
384 : // If rowflag is set, flag all chans in it
385 0 : if(rowflags(tm*NumB+bs))
386 : {
387 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
388 0 : flagc(pl,ch,tm*NumB+bs) = true;
389 : }
390 :
391 : // Count flags across channels, and also count the data.
392 0 : temp=0;
393 : //flg=true;
394 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
395 : {
396 0 : temp += visc(pl,ch,tm*NumB+bs);
397 : //flg &= flagc(pl,ch,tm*NumB+bs);
398 : }
399 :
400 : // If data is zero for all channels (not read in), set flags to zero.
401 0 : if(temp==0)
402 : {
403 0 : rowflags(tm*NumB+bs)=true;
404 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
405 0 : flagc(pl,ch,tm*NumB+bs)=true;
406 : }
407 :
408 : // Count flags across channels and time,,,,,
409 : // If any flag is false, bflag will become false
410 0 : for(casacore::uInt ch=0; ch<NumC; ch++)
411 0 : bflag &= flagc(pl,ch,tm*NumB+bs);
412 :
413 : }// for tm
414 : // If all times/chans are flagged for this baseline, set the baselineflag.
415 0 : if(bflag) baselineflags(bs)=true;
416 0 : else baselineflags(bs)=false;
417 : }// for bs
418 0 : casacore::Int ubs=0;
419 0 : for(casacore::uInt bs=0;bs<NumB;bs++)
420 0 : if(!baselineflags(bs)) ubs++;
421 0 : if(ShowPlots) cout << "Working with " << ubs << " unflagged baseline(s). " << endl;
422 : }// for pl
423 0 : }// end of FlagZeros()
424 :
425 :
426 :
427 :
428 0 : void RFATimeFreqCrop :: RunTFCrop()
429 : {
430 : casacore::uInt a1,a2;
431 :
432 0 : meanBP = 0;
433 0 : cleanBP = 0;
434 :
435 0 : for(casacore::uInt pl=0;pl<NumP;pl++)
436 : {
437 0 : for(casacore::uInt bs=0;bs<NumB;bs++)
438 : {
439 0 : if( !baselineflags(bs) )
440 : {
441 0 : Ants(bs,&a1,&a2);
442 0 : if((CorrChoice==0 && a1 == a2)||(CorrChoice!=0 && a1 != a2))
443 : {
444 :
445 0 : FlagTimeSeries(pl,bs);
446 :
447 0 : FitCleanBandPass(pl,bs);
448 :
449 0 : FlagBandPass(pl,bs);
450 :
451 0 : GrowFlags(pl,bs);
452 : }// if corrchoice
453 : }// if baseline is not flagged
454 : }// end for bs
455 : }// end for pl
456 :
457 0 : }// end runTFCrop
458 :
459 :
460 :
461 :
462 : /* Flag in time, and build the average bandpass */
463 : /* Grow flags by one timestep, check for complete flagged baselines. */
464 0 : void RFATimeFreqCrop :: FlagTimeSeries(casacore::uInt pl, casacore::uInt bs)
465 : {
466 : //casacore::Float mn=0;
467 0 : casacore::Float sd=0,temp=0,flagcnt=0,tol=0;
468 : casacore::uInt a1,a2;
469 : //casacore::Bool flg=false;
470 : /* For each Channel - fit lines to 1-D data in time - flag according
471 : * to them and build up the mean bandpass */
472 :
473 0 : casacore::Float rmean=0;
474 0 : Ants(bs,&a1,&a2);
475 : // cout << " Antennas : " << a1 << " & " << a2 << endl;
476 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
477 : {
478 0 : tempTS=0;flagTS=false;
479 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
480 : {
481 0 : tempTS[tm] = visc(pl,ch,((tm*NumB)+bs));
482 0 : flagTS[tm] = flagc(pl,ch,((tm*NumB)+bs));
483 : }//for tm
484 :
485 :
486 0 : rmean += UMean(tempTS,flagTS);
487 :
488 0 : temp=0;
489 0 : for(int loop=0;loop<5;loop++)
490 : {
491 : // UUU : HERE - give choice of PolyFit in time...
492 0 : LineFit(tempTS,flagTS,fitTS,0,tempTS.nelements()-1);
493 0 : sd = UStd(tempTS,flagTS,fitTS);
494 :
495 0 : for(casacore::uInt i=0;i<NumT;i++)
496 0 : if(flagTS[i]==false && fabs(tempTS[i]-fitTS[i]) > T_TOL*sd)
497 : {
498 0 : flagTS[i]=true ;flagcnt++;
499 : }
500 0 : if(fabs(temp-sd) < 0.1)break;
501 0 : temp=sd;
502 : }
503 :
504 : // If sum of 2 adjacent flags also crosses threshold, flag
505 : /*
506 : for(casacore::uInt i=1;i<NumT-1;i++)
507 : {
508 : if(flagTS[i])
509 : {
510 : if( ( fabs(tempTS[i-1]-fitTS[i-1]) + fabs(tempTS[i+1]-fitTS[i+1]) ) > T_TOL*sd )
511 : {flagTS[i-1]=true; flagTS[i+1]=true;}
512 : }
513 : }
514 : */
515 :
516 0 : meanBP(pl,bs,ch) = UMean(tempTS,flagTS) ;
517 :
518 : /* write flags to local flag cube */
519 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
520 0 : flagc(pl,ch,((tm*NumB)+bs))=flagTS[tm];
521 :
522 : }//for ch
523 :
524 :
525 : /* Check for completely flagged ants/bs */
526 : if(1)
527 : {
528 0 : if((CorrChoice==0 && a1 == a2)||(CorrChoice!=0 && a1 != a2))
529 : {
530 0 : if(CorrChoice==0)tol=ANT_TOL;
531 0 : else tol=BASELN_TOL;
532 0 : if(fabs(rmean/float(NumC)) < tol)
533 : {
534 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
535 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
536 0 : flagc(pl,ch,((tm*NumB)+bs))=true;
537 0 : if(CorrChoice==0)
538 0 : cout << "Antenna Flagged : " << a1 << endl;
539 : else
540 0 : cout << "Mean : " << rmean/NumC << " : Baseline Flagged : " << a1 << ":" << a2 << endl;
541 : }
542 :
543 : }
544 : }///if(0);
545 0 : }// end of FlagTimeSeries
546 :
547 : /* Fit a smooth bandpass to the mean bandpass and store it
548 : * one for each baseline */
549 :
550 : // matpos : NumP. NumB. NumC
551 :
552 0 : void RFATimeFreqCrop :: FitCleanBandPass(casacore::uInt pl, casacore::uInt bs)
553 : {
554 : //casacore::Float mn=0,sd=0,temp=0,flagcnt=0,tol=0;
555 : casacore::uInt a1,a2;
556 0 : casacore::Bool flg=false;
557 :
558 0 : Ants(bs,&a1,&a2);
559 : /* Fit a smooth bandpass */
560 0 : flg=true;
561 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
562 : {
563 0 : tempBP[ch] = meanBP(pl,bs,ch);
564 0 : if(tempBP[ch] != 0) flg=false;
565 : }
566 :
567 0 : if(flg==false)
568 : {
569 : /* Piecewise Poly Fit to the meanBP */
570 0 : if(!FreqLineFit)
571 : {
572 0 : CleanBand(tempBP,fitBP);
573 : }
574 : else
575 : {
576 : /* LineFit to flag off a line fit in frequency */
577 0 : flagBP=false;
578 0 : for(casacore::uInt ch=0;ch<tempBP.nelements();ch++)
579 0 : if(tempBP[ch]==0) flagBP[ch]=true;
580 0 : LineFit(tempBP,flagBP,fitBP,0,tempBP.nelements()-1);
581 : }
582 :
583 : #ifdef PLOT
584 : if(CorrChoice==0)
585 : cout<<" Antenna : "<<bs<<" Polzn : "<<pl<<endl;
586 : else
587 : {
588 : Ants(bs,&a1,&a2);
589 : cout << " Baseline : " << a1 << ":" << a2 << " Polzn : " << pl << endl;
590 : }
591 : Plot_ds9(tempBP.nelements(), tempBP,fitBP); // Plot the band
592 : #endif
593 : }
594 : /* else
595 : {
596 : Ants(bs,&a1,&a2);
597 : emptylist += casacore::String::toString(a1)+"-"+casacore::String::toString(a2)+" ";
598 : // cout << "meanBP is filled with zeros : baseline : " << a1 << "-" << a2 << endl;
599 : }
600 : */
601 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
602 : {
603 0 : if(flg==false) cleanBP(pl,bs,ch)= fitBP[ch];
604 0 : else cleanBP(pl,bs,ch)=0;
605 : }
606 :
607 0 : }// end FitCleanBandPass
608 :
609 : /* FLAGGING IN FREQUENCY */
610 0 : void RFATimeFreqCrop :: FlagBandPass(casacore::uInt pl, casacore::uInt bs)
611 : {
612 0 : casacore::Float mn=0,sd=0,temp=0,flagcnt=0;
613 : casacore::uInt a1,a2;
614 : //casacore::Bool flg=false;
615 :
616 :
617 0 : Ants(bs,&a1,&a2);
618 :
619 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
620 : {
621 : /* Divide (or subtract) out the clean bandpass */
622 0 : tempBP=0,flagBP=0;
623 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
624 : {
625 0 : flagBP[ch] = flagc(pl,ch,((tm*NumB)+bs));
626 0 : if(flagBP[ch]==false)
627 0 : tempBP[ch] = visc(pl,ch,((tm*NumB)+bs))/cleanBP(pl,bs,ch);
628 : }//for ch
629 :
630 : /* Flag outliers */
631 0 : temp=0;
632 0 : for(casacore::Int loop=0;loop<5;loop++)
633 : {
634 0 : mn=1;
635 0 : sd = UStd(tempBP,flagBP,mn);
636 :
637 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
638 : {
639 0 : if(flagBP[ch]==false && fabs(tempBP[ch]-mn) > F_TOL*sd)
640 : {
641 0 : flagBP[ch]=true ;flagcnt++;
642 : }
643 : }
644 0 : if(fabs(temp-sd) < 0.1)break;
645 0 : temp=sd;
646 : }
647 :
648 : /* If sum of power in two adjacent channels is more than thresh, flag both side chans */
649 0 : if(FlagLevel>0)
650 : {
651 0 : for(casacore::uInt ch=1;ch<NumC-1;ch++)
652 : {
653 0 : if(flagBP[ch])
654 : {
655 0 : if( ( fabs(tempBP[ch-1]-mn) + fabs(tempBP[ch+1]-mn) ) > F_TOL*sd )
656 0 : {flagBP[ch-1]=true; flagBP[ch+1]=true;}
657 : }
658 : }
659 : }
660 :
661 : /* Fill the flags into the visbuffer array */
662 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
663 0 : flagc(pl,ch,((tm*NumB)+bs))=flagBP[ch];
664 :
665 :
666 : }//for tm
667 :
668 0 : }// end FlagBandPass
669 :
670 : /* APPLY FLAG HEURISTICS ON THE FLAGS FOR ALL AUTOCORRELATIONS */
671 0 : void RFATimeFreqCrop :: GrowFlags(casacore::uInt pl, casacore::uInt bs)
672 : {
673 : casacore::uInt a1,a2;
674 : //casacore::Bool flg=false;
675 :
676 0 : if(FlagLevel > 0)
677 : {
678 0 : Ants(bs,&a1,&a2);
679 :
680 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
681 : {
682 : // casacore::uInt fsum=0;
683 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
684 : {
685 0 : if(FlagLevel>1) // flag level 2 and above...
686 : {
687 : // flagging one timestamp before and after
688 0 : if(tm>0)
689 0 : if(flagc(pl,ch,((tm*NumB+bs)))==true)
690 0 : flagc(pl,ch,(((tm-1)*NumB+bs)))=true;
691 :
692 0 : if((NumT-tm)<NumT-1)
693 0 : if(flagc(pl,ch,(((NumT-tm)*NumB+bs)))==true)
694 0 : flagc(pl,ch,(((NumT-tm+1)*NumB+bs)))=true;
695 : }
696 :
697 0 : if(FlagLevel>1) // flag level 2 and above
698 : {
699 : // flagging one channel before and after
700 0 : if(ch>0)
701 0 : if(flagc(pl,ch,((tm*NumB+bs)))==true)
702 0 : flagc(pl,ch-1,((tm*NumB+bs)))=true;
703 :
704 0 : if((NumC-ch)<NumC-1)
705 0 : if(flagc(pl,(NumC-ch),(tm*NumB+bs))==true)
706 0 : flagc(pl,(NumC-ch+1),(tm*NumB+bs))=true;
707 : }
708 :
709 0 : if(FlagLevel>0) // flag level 1 and above
710 : {
711 : /* If previous and next channel are flagged, flag it */
712 0 : if(ch>0 && ch < NumC-1)
713 : {
714 0 : if( flagc(pl,ch-1,(tm*NumB+bs) ) == true
715 0 : && flagc(pl,ch+1,(tm*NumB+bs) ) == true )
716 0 : flagc(pl,ch,(tm*NumB+bs) ) = true;
717 : }
718 : /* If previous and next timestamp are flagged, flag it */
719 0 : if(tm>0 && tm < NumT-1)
720 : {
721 0 : if( flagc(pl,ch,((tm-1)*NumB+bs) ) == true
722 0 : && flagc(pl,ch,((tm+1)*NumB+bs) ) == true )
723 0 : flagc(pl,ch,(tm*NumB+bs) ) = true;
724 : }
725 : }
726 :
727 0 : if(FlagLevel>1) // flag level 2 and above
728 : {
729 : /* If next two channels are flagged, flag it */
730 0 : if(ch < NumC-2)
731 0 : if( flagc(pl,ch+1,(tm*NumB+bs)) == true
732 0 : && flagc(pl,ch+2,(tm*NumB+bs) ) == true )
733 0 : flagc(pl,ch,(tm*NumB+bs) ) = true;
734 : }
735 :
736 : }//for tm
737 :
738 : // if more than 60% of the timetange flagged - flag whole timerange for that channel
739 : //if(fsum < 0.4*NumT && FlagLevel > 1) // flag level 2
740 : // for(casacore::uInt tm=0;tm<NumT;tm++)
741 : // flagc(pl,ch,((tm*NumB+bs)))=true;
742 : }//for ch
743 :
744 0 : if(FlagLevel>0) // flag level 1 and above
745 : {
746 : /* If more than 4 surrounding points are flagged, flag it */
747 0 : casacore::uInt fsum2=0;
748 0 : for(casacore::uInt ch=1;ch<NumC-1;ch++)
749 : {
750 0 : for(casacore::uInt tm=1;tm<NumT-1;tm++)
751 : {
752 0 : fsum2 = (casacore::uInt)(flagc(pl,ch-1,(((tm-1)*NumB+bs)))) + (casacore::uInt)(flagc(pl,ch-1,((tm*NumB+bs)))) +
753 0 : (casacore::uInt)(flagc(pl,ch-1,(((tm+1)*NumB+bs)))) + (casacore::uInt)(flagc(pl,ch,(((tm-1)*NumB+bs)))) +
754 0 : (casacore::uInt)(flagc(pl,ch,(((tm+1)*NumB+bs)))) + (casacore::uInt)(flagc(pl,ch+1,(((tm-1)*NumB+bs)))) +
755 0 : (casacore::uInt)(flagc(pl,ch+1,((tm*NumB+bs)))) + (casacore::uInt)(flagc(pl,ch+1,(((tm+1)*NumB+bs))));
756 0 : if(fsum2 > 4) flagc(pl,ch,((tm*NumB+bs))) = true;
757 : } // for tm
758 : }// for ch
759 : }// if FlagLevel>0
760 :
761 0 : if(FlagLevel>0) // flaglevel = 1 and above
762 : {
763 0 : casacore::uInt fsum2=0;
764 : /* Grow flags in time */
765 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
766 : {
767 0 : fsum2=0;
768 : /* count unflagged points for this channel (all times) */
769 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
770 0 : fsum2 += (flagc(pl,ch,((tm*NumB+bs)))==true)?0:1 ;
771 : /*if more than 50% of the timetange flagged - flag whole timerange for that channel */
772 0 : if(fsum2 < 0.5*NumT)
773 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
774 0 : flagc(pl,ch,((tm*NumB+bs)))=true;
775 : }// for ch
776 : }// if flaglevel>0
777 :
778 : }//if flag_level
779 :
780 : // Count flags
781 : /*
782 : casacore::Float runningcount=0, runningflag=0;
783 : runningcount=0;
784 : runningflag=0;
785 : for(int pl=0;pl<NumP;pl++)
786 : {
787 : for(casacore::uInt bs=0;bs<NumB;bs++)
788 : {
789 : Ants(bs,&a1,&a2);
790 : if(CorrChoice==0)
791 : {if(a1 != a2) continue;} // choose autocorrelations
792 : else
793 : {if(a1==a2) continue;} // choose cross correlations
794 :
795 : for(int ch=0;ch<NumC;ch++)
796 : {
797 : for(casacore::uInt tm=0;tm<NumT;tm++)
798 : {
799 : runningflag += casacore::Float(flagc(pl,ch,(tm*NumB)+bs));
800 : runningcount++ ;
801 : }// for tm
802 : }//for ch
803 : }// for bs
804 : }// for pl
805 :
806 :
807 : cout << " Flagged : " << 100 * runningflag/runningcount << " % for timesteps " << iterTimecnt-NumT << " - " << iterTimecnt << endl;
808 : */
809 0 : }// end of GrowFlags
810 :
811 : /* GRAY SCALE DISPLAYS ON ds9 */
812 0 : RFA::IterMode RFATimeFreqCrop :: ShowFlagPlots()
813 : {
814 : casacore::uInt a1,a2;
815 0 : if(ShowPlots)
816 : {
817 :
818 : // cout << "About to display : allocating cubes" << endl;
819 : /*
820 : casacore::Float **dispdat=NULL,**flagdat=NULL;
821 : dispdat = (casacore::Float **) malloc(sizeof(casacore::Float *) * NumT + sizeof(casacore::Float)*NumC*NumT);
822 : for(casacore::uInt i=0;i<NumT;i++)
823 : dispdat[i] = (casacore::Float *) (dispdat + NumT) + i * NumC;
824 :
825 : flagdat = (casacore::Float **) malloc(sizeof(casacore::Float *) * NumT + sizeof(casacore::Float)*NumC*NumT);
826 : for(casacore::uInt i=0;i<NumT;i++)
827 : flagdat[i] = (casacore::Float *) (flagdat + NumT) + i * NumC;
828 : */
829 :
830 0 : casacore::IPosition shp(2),tshp(2); shp(0)=NumC; shp(1)=NumT;
831 0 : casacore::Matrix<casacore::Float> dispdat(shp), flagdat(shp);
832 :
833 : // cout << "About to display : allocated. " << endl;
834 :
835 0 : char choice = 'a';
836 :
837 0 : casacore::Float runningsum=0, runningflag=0, oldrunningflag=0;
838 0 : for(casacore::uInt pl=0;pl<NumP;pl++)
839 : {
840 0 : if(choice == 's') continue;
841 0 : for(casacore::uInt bs=0;bs<NumB;bs++)
842 : {
843 0 : if(choice == 's') continue;
844 0 : if(baselineflags(bs)) continue;
845 :
846 0 : runningsum=0;
847 0 : runningflag=0;
848 0 : oldrunningflag=0;
849 0 : Ants(bs,&a1,&a2);
850 0 : if(CorrChoice==0)
851 0 : {if(a1 != a2) continue;} // choose autocorrelations
852 : else
853 0 : {if(a1==a2) continue;} // choose cross correlations
854 :
855 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
856 : {
857 0 : tempBP[ch] = meanBP(pl,bs,ch);
858 0 : fitBP[ch] = cleanBP(pl,bs,ch);
859 : }
860 :
861 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
862 : {
863 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
864 : {
865 : // casacore::Data with pre-flags
866 0 : dispdat(ch,tm) = visc(pl,ch,(((tm*NumB)+bs))) * (!chunkflags(pl,ch,((tm+iterTimecnt-NumT)*NumB)+bs));
867 : // casacore::Data with all flags (pre and new)
868 0 : flagdat(ch,tm) = dispdat(ch,tm)*(!flagc(pl,ch,(tm*NumB)+bs));
869 : // Sum of the visibilities (all of them, flagged and unflagged)
870 0 : runningsum += visc(pl,ch,(((tm*NumB)+bs)));
871 : /*
872 : // Count of all flags
873 : runningflag += (casacore::Float)(flagc(pl,ch,(tm*NumB)+bs));
874 : // Count of only pre-flags
875 : oldrunningflag += (casacore::Float)(chunkflags(pl,ch,((tm+iterTimecnt-NumT)*NumB)+bs));
876 : */ ////// CHECK that iterTimecnt is correct, and a valid part of chunkflags is being read !!!!!!!!!
877 : // Count of all flags
878 0 : if( (flagc(pl,ch,(tm*NumB)+bs)) ) runningflag++;
879 : // Count of only pre-flags
880 0 : if(chunkflags(pl,ch,((tm+iterTimecnt-NumT)*NumB)+bs)) oldrunningflag++;
881 : }//for tm
882 : }//for ch
883 :
884 : //Plot on ds9 !!
885 :
886 : // cout << "Antenna1 : " << a1 << " Antenna2 : " << a2 << " Polarization : " << pl << endl;
887 : // cout << "Vis sum : " << runningsum << " Flag % : " << 100 * runningflag/(NumC*NumT) << endl;
888 : // cout << " Flagged : " << 100 * runningflag/(NumC*NumT) << " %" << endl;
889 0 : cout << " Flagged : " << 100 * runningflag/(NumC*NumT) << " % (Pre-Flag : " << 100 * oldrunningflag/(NumC*NumT) << " %) on " << Expr << " for timesteps " << iterTimecnt-NumT << " - " << iterTimecnt << " on baseline " << a1 << "-" << a2;
890 :
891 0 : if(!runningsum)
892 : {
893 0 : cout << " : No non-zero data !" << endl;
894 : }
895 : else
896 : {
897 0 : cout << endl;
898 :
899 0 : Display_ds9(NumC,NumT,dispdat,1);
900 0 : Display_ds9(NumC,NumT,flagdat,2);
901 0 : Plot_ds9(tempBP.nelements(), tempBP, fitBP);
902 : //UPlot(tempBP,fitBP,0,tempBP.nelements()-1); // Plot the band
903 :
904 0 : cout << "Press <c> to continue display, <s> to stop display but continue flagging, <q> to quit." << endl;
905 : //getchar();
906 0 : cin >> choice;
907 : // cout << " Choice : " << choice << endl;
908 0 : switch(choice)
909 : {
910 0 : case 'q':
911 0 : ShowPlots = false;
912 0 : StopAndExit = true;
913 0 : cout << "Exiting flagger" << endl;
914 0 : return RFA::STOP;
915 0 : case 's':
916 0 : ShowPlots = false;
917 0 : cout << "Stopping display. Continuing flagging." << endl;
918 0 : break;
919 0 : default:
920 0 : break;
921 : }
922 : }
923 : }//for bs
924 : }//for pl
925 :
926 0 : }// end of if ShowPlots
927 0 : return RFA::CONT;
928 : }// end ShowFlagPlots
929 :
930 : /* Extend Flags
931 : (1) If flagging on self-correlations, extend to cross-corrs.
932 : (2) If requested, extend across polarization (fiddle with corrmask)
933 : (3) If requested, extend across baseline/antenna
934 : */
935 0 : void RFATimeFreqCrop :: ExtendFlags()
936 : {
937 : casacore::uInt a1,a2;
938 :
939 : /* FLAG BASELINES FROM THE SELF FLAGS */
940 0 : if(CorrChoice ==0)
941 : {
942 0 : cout << " Flagging Cross correlations from self correlation flags " << endl;
943 0 : for(casacore::uInt pl=0;pl<NumP;pl++)
944 : {
945 0 : for(casacore::uInt bs=0;bs<NumB;bs++)
946 : {
947 0 : if(baselineflags(bs)) continue;
948 0 : Ants(bs,&a1,&a2);
949 0 : if(a1 == a2) continue; // choose cross correlations
950 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
951 : {
952 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
953 : {
954 0 : flagc(pl,ch,((tm*NumB+bs))) = flagc(pl,ch,((tm*NumB)+SELF(a1))) | flagc(pl,ch,((tm*NumB)+SELF(a1)));
955 : }//for tm
956 : }//for ch
957 : }//for bs
958 : }//for pl
959 : }
960 :
961 0 : }// end Extend Flags
962 :
963 0 : void RFATimeFreqCrop :: FillChunkFlags()
964 : {
965 : //cout << " Diagnostics on flag cube " << endl;
966 :
967 0 : for(casacore::uInt pl=0;pl<NumP;pl++)
968 : {
969 0 : for(casacore::uInt bs=0;bs<NumB;bs++)
970 : {
971 : //if(RowFlags(pl,(tm*NumB)+bs)==true)
972 : // continue;
973 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
974 : {
975 0 : for(casacore::uInt tm=0;tm<NumT;tm++)
976 : {
977 0 : if(flagc(pl,ch,((tm*NumB)+bs))==true )
978 : {
979 0 : chunkflags(pl,ch,((tm+iterTimecnt-NumT)*NumB)+bs) = true;
980 : }
981 : }// for tm
982 : }//for ch
983 : }// for bs
984 : }// for pl
985 :
986 : // what is this ??
987 0 : timecnt=0;
988 :
989 0 : }// end of FillChunkFlags
990 :
991 :
992 :
993 : /*RedFlagger::run - ends loop for all agents with endData
994 : * Calls endData once at the end of each PASS */
995 0 : RFA::IterMode RFATimeFreqCrop :: endData ()
996 : {
997 : // cout << " In End Data. Ending timecnt : " << timecnt << endl;
998 :
999 0 : RFAFlagCubeBase::endData();
1000 0 : return RFA::STOP;
1001 : }
1002 :
1003 :
1004 : /* Called at the beginning of each PASS */
1005 0 : void RFATimeFreqCrop :: startFlag (bool verbose)
1006 : {
1007 : // corrmask = RFDataMapper::corrMask(chunk.visIter());
1008 : // cout << "StartFlag : corrmask : " << corrmask << endl;
1009 :
1010 0 : RFAFlagCubeBase::startFlag(verbose);
1011 :
1012 0 : }
1013 :
1014 :
1015 : /* Write Flags to casacore::MS */
1016 0 : void RFATimeFreqCrop :: iterFlag(casacore::uInt itime)
1017 : {
1018 : // cout << "iterFlag :: Set flags for time : " << itime << endl;
1019 :
1020 : // FLAG DATA
1021 :
1022 0 : if(!DryRun)
1023 : {
1024 0 : flag.advance(itime);
1025 0 : corrmask = RFDataMapper::corrMask(chunk.visIter());
1026 :
1027 0 : const casacore::Vector<casacore::Int> &ifrs( chunk.ifrNums() );
1028 0 : if(ant1.shape() != (vb.antenna1()).shape())
1029 0 : ant1.resize((vb.antenna1()).shape());
1030 0 : ant1 = vb.antenna1();
1031 :
1032 0 : if(ant2.shape() != (vb.antenna2()).shape())
1033 0 : ant2.resize((vb.antenna2()).shape());
1034 0 : ant2 = vb.antenna2();
1035 :
1036 0 : casacore::uInt nbs = ant1.nelements();
1037 0 : for(casacore::uInt pl=0;pl<NumP;pl++)
1038 : {
1039 0 : for(casacore::uInt bs=0;bs<nbs;bs++)
1040 : {
1041 0 : casacore::Bool bflag=true;
1042 0 : casacore::uInt baselinecnt = BaselineIndex(bs,ant1[bs],ant2[bs]);
1043 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
1044 : {
1045 0 : if(chunkflags(pl,ch,(itime*NumB)+baselinecnt)==(casacore::Bool)true)
1046 0 : flag.setFlag(ch,ifrs(bs));
1047 0 : bflag &= chunkflags(pl,ch,(itime*NumB)+baselinecnt);
1048 : }
1049 0 : if(bflag) flag.setRowFlag(ifrs(bs),itime);
1050 : }
1051 : }
1052 :
1053 0 : flag.setMSFlags(itime);
1054 :
1055 : }// if not dry-run
1056 : else
1057 : {
1058 0 : RFAFlagCubeBase::iterFlag(itime);
1059 : }
1060 :
1061 : /// FLAG ROWS
1062 : /*
1063 : ant1 = vb.antenna1();
1064 : ant2 = vb.antenna2();
1065 : casacore::uInt npols = (chunkflags.shape())[0];
1066 : casacore::uInt nbs = ant1.nelements();
1067 :
1068 : vi.flag(ff);
1069 : vi.flagRow(fr);
1070 :
1071 : for(casacore::uInt bs=0;bs<nbs;bs++)
1072 : {
1073 : casacore::Bool rowflag=true;
1074 : for(casacore::uInt pl=0;pl<npols;pl++)
1075 : {
1076 : for(casacore::uInt ch=StartChan;ch<=EndChan;ch++)
1077 : {
1078 : ff(pl,ch,bs) = chunkflags(pl,ch-StartChan,(itime*NumB)+baselinecnt);
1079 : rowflag &= ff(pl,ch,bs);
1080 : }
1081 : }
1082 : if(rowflag) fr[bs]=true;
1083 : }
1084 :
1085 : //vi.setFlag(ff);
1086 : //vi.setFlagRow(fr);
1087 :
1088 : chunk.visIter().setFlag(ff);
1089 : chunk.visIter().setFlagRow(fr);
1090 :
1091 : */
1092 0 : }
1093 :
1094 : /*RedFlagger::run - calls 'endChunk()' on all agents. */
1095 :
1096 0 : void RFATimeFreqCrop :: endChunk ()
1097 : {
1098 0 : casacore::LogIO os(casacore::LogOrigin("tfcrop","endChunk","WHERE"));
1099 : // cout << "endChunk : counting flags" << endl;
1100 : // Count flags
1101 0 : if(!StopAndExit)
1102 : {
1103 0 : casacore::Float runningcount=0, runningflag=0;
1104 :
1105 : //const casacore::Vector<casacore::Int> &ifrs( chunk.ifrNums() );
1106 0 : if(ant1.shape() != (vb.antenna1()).shape())
1107 0 : ant1.resize((vb.antenna1()).shape());
1108 0 : ant1 = vb.antenna1();
1109 :
1110 0 : if(ant2.shape() != (vb.antenna2()).shape())
1111 0 : ant2.resize((vb.antenna2()).shape());
1112 0 : ant2 = vb.antenna2();
1113 :
1114 0 : casacore::uInt nbs = ant1.nelements();
1115 0 : for(casacore::uInt pl=0;pl<NumP;pl++)
1116 : {
1117 0 : for(casacore::uInt bs=0;bs<nbs;bs++)
1118 : {
1119 0 : casacore::uInt baselinecnt = BaselineIndex(bs,ant1[bs],ant2[bs]);
1120 0 : for(casacore::uInt ch=0;ch<NumC;ch++)
1121 : {
1122 0 : for(casacore::uInt tm=0;tm<num(TIME);tm++)
1123 : {
1124 0 : if (chunkflags(pl,ch,((tm)*NumB)+baselinecnt)) runningflag++;
1125 0 : runningcount++;
1126 : }
1127 : }
1128 : }
1129 : }
1130 :
1131 : /* Polarizations/correlations */
1132 0 : corrmask = RFDataMapper::corrMask(chunk.visIter());
1133 0 : casacore::Vector<casacore::Int> corrlist(num(POLZN));
1134 0 : for(casacore::uInt corr=0; corr<num(POLZN); corr++)
1135 0 : corrlist[corr] = (casacore::Int) ( (corrmask >> corr) & 1 );
1136 :
1137 : // cout << "--> Flagged " << 100 * runningflag/runningcount << " % on " << Expr << ". Applying to correlations : " << corrlist << endl;
1138 0 : os << "TFCROP : Flagged " << 100 * runningflag/runningcount << " % on " << Expr << ". Applying to corrs : " << corrlist;
1139 0 : if(DryRun) os << " (Not writing flags to MS)" << casacore::LogIO::POST;
1140 0 : else os << " (Writing flags to MS)" << casacore::LogIO::POST;
1141 0 : }
1142 0 : RFAFlagCubeBase::endChunk();
1143 : /// (chunk.visIter()).setRowBlocking(0); //reset to default
1144 : // cout << " End of endChunk !!" << endl;
1145 0 : }
1146 :
1147 :
1148 : /* Destructor for RFATimeFreqCrop */
1149 :
1150 :
1151 0 : RFATimeFreqCrop :: ~RFATimeFreqCrop ()
1152 : {
1153 : //cout << "destructor for RFATimeFreqCrop" << endl;
1154 0 : }
1155 :
1156 :
1157 :
1158 :
1159 : /* Calculate the MEAN of 'vect' ignoring values flagged in 'flag' */
1160 0 : casacore::Float RFATimeFreqCrop :: UMean(casacore::Vector<casacore::Float> vect, casacore::Vector<casacore::Bool> flag)
1161 : {
1162 0 : casacore::Float mean=0;
1163 0 : int cnt=0;
1164 0 : for(int i=0;i<(int)vect.nelements();i++)
1165 0 : if(flag[i]==false)
1166 : {
1167 0 : mean += vect[i];
1168 0 : cnt++;
1169 : }
1170 0 : if(cnt==0) cnt=1;
1171 0 : return mean/cnt;
1172 : }
1173 :
1174 :
1175 : /* Calculate the STANDARD DEVN. of 'vect' w.r.to a given 'fit'
1176 : * ignoring values flagged in 'flag' */
1177 0 : casacore::Float RFATimeFreqCrop :: UStd(casacore::Vector<casacore::Float> vect, casacore::Vector<casacore::Bool> flag, casacore::Vector<casacore::Float> fit)
1178 : {
1179 0 : casacore::Float std=0;
1180 0 : int n=0,cnt=0;
1181 0 : n = vect.nelements() < fit.nelements() ? vect.nelements() : fit.nelements();
1182 0 : for(int i=0;i<n;i++)
1183 0 : if(flag[i]==false)
1184 : {
1185 0 : cnt++;
1186 0 : std += (vect[i]-fit[i])*(vect[i]-fit[i]);
1187 : }
1188 0 : if(cnt==0) cnt=1;
1189 0 : return sqrt(std/cnt);
1190 : }
1191 :
1192 :
1193 : /* Calculate the STANDARD DEVN. of 'vect' w.r.to a given mean
1194 : * ignoring values flagged in 'flag' */
1195 0 : casacore::Float RFATimeFreqCrop :: UStd(casacore::Vector<casacore::Float> vect, casacore::Vector<casacore::Bool> flag, casacore::Float mean)
1196 : {
1197 0 : casacore::Float std=0;
1198 0 : int cnt=0;
1199 0 : for(int i=0;i<(int)vect.nelements();i++)
1200 0 : if(flag[i]==false)
1201 : {
1202 0 : cnt++;
1203 0 : std += (vect[i]-mean)*(vect[i]-mean);
1204 : }
1205 0 : return sqrt(std/cnt);
1206 : }
1207 :
1208 : /* Fit Piecewise polynomials to 'data' and get the 'fit' */
1209 0 : void RFATimeFreqCrop :: CleanBand(casacore::Vector<casacore::Float> data,casacore::Vector<casacore::Float> fit)
1210 : {
1211 : // casacore::Int step=0,ind=0;
1212 0 : casacore::Int deg=0; //start=0;
1213 0 : casacore::Int left=0,right=0;
1214 : // casacore::Int le=0,ri=0;
1215 0 : casacore::Float sd,TOL=3;
1216 0 : casacore::Vector<casacore::Float> tdata;
1217 0 : casacore::Vector<casacore::Bool> tfband;
1218 :
1219 0 : tfband.resize(data.nelements());
1220 0 : tdata.resize(data.nelements());
1221 :
1222 0 : tfband = false;
1223 0 : tdata = data;
1224 :
1225 : /* replace empty data values by adjacent values */
1226 0 : for(casacore::uInt i=0;i<tdata.nelements();i++)
1227 : {
1228 0 : if(tdata[i]==0)
1229 : {
1230 0 : if(i==0)// find first non-zero value and set to that.
1231 : {
1232 0 : casacore::Int ind=0;
1233 0 : for(casacore::uInt j=1;j<tdata.nelements();j++)
1234 0 : if(tdata[j]!=0){ind=j;break;}
1235 0 : if(ind==0) tdata[i]=0;
1236 0 : else tdata[i]=tdata[ind];
1237 : }
1238 : else// find next non-zero value and interpolate.
1239 : {
1240 0 : casacore::Int indr=0;
1241 0 : for(casacore::uInt j=i+1;j<tdata.nelements();j++)
1242 0 : if(tdata[j]!=0){indr=j;break;}
1243 0 : casacore::Int indl=-1;
1244 0 : for(int j = i ; j >= 0 ; j--)
1245 0 : if(tdata[j]!=0){indl=j;break;}
1246 :
1247 0 : if(indl==-1 && indr==0) tdata[i]=0;
1248 0 : if(indl>-1 && indr==0) tdata[i]=tdata[indl];
1249 0 : if(indl==-1 && indr>0) tdata[i]=tdata[indr];
1250 0 : if(indl>-1 && indr>0) tdata[i]=(tdata[indl]+tdata[indr])/2.0;
1251 : }
1252 : }
1253 : }
1254 :
1255 :
1256 : /* If there still are empty points (entire spectrum is flagged) flag it. */
1257 0 : for(casacore::uInt i=0;i<tdata.nelements();i++)
1258 0 : if(tdata[i]==0)
1259 : {
1260 : //cout << "chan " << i << " is blank" << endl;
1261 0 : tfband[i]=true;
1262 : }
1263 :
1264 0 : fit = tdata;
1265 :
1266 0 : casacore::Int psize=1;
1267 0 : casacore::Int leftover=1,leftover_front=0,npieces=1;
1268 :
1269 0 : deg=1;
1270 0 : npieces=1;
1271 :
1272 0 : for(casacore::uInt j=0;j<=4;j++)
1273 : {
1274 : // if(j==0) {deg = 1;npieces=1;}
1275 : // if(j==1) {deg = 1;npieces=5;}
1276 : // if(j==2) {deg = 2;npieces=6;}
1277 : // if(j==3) {deg = 3;npieces=7;}
1278 : // if(j==4) {deg = 3;npieces=8;}
1279 :
1280 0 : npieces = MIN(2*j+1, MaxNPieces);
1281 0 : if(j>1) {deg=2;}
1282 0 : if(j>2) {deg=3;}
1283 :
1284 0 : psize = (int)(tdata.nelements()/npieces);
1285 : // cout << "Iter : " << j << " with Deg : " << deg << " and Piece-size : " << psize << endl;
1286 :
1287 0 : leftover = (int)(tdata.nelements() % npieces);
1288 :
1289 0 : leftover_front = (int)(leftover/2.0);
1290 :
1291 0 : left=0; right=tdata.nelements()-1;
1292 0 : for(casacore::Int p=0;p<npieces;p++)
1293 : {
1294 0 : if(npieces>1)
1295 : {
1296 0 : left = leftover_front + p*psize;
1297 0 : right = left + psize;
1298 :
1299 0 : if(p==0) {left = 0; right = leftover_front + psize;}
1300 0 : if(p==npieces-1) {right = tdata.nelements()-1;}
1301 : }
1302 0 : if(deg==1)
1303 0 : LineFit(tdata,tfband,fit,left,right);
1304 : else
1305 0 : PolyFit(tdata,tfband,fit,left,right,deg);
1306 : }
1307 :
1308 : /* Now, smooth the fit - make this nicer later */
1309 :
1310 0 : int winstart=0, winend=0;
1311 0 : float winsum=0.0;
1312 0 : int offset=2;
1313 0 : for(casacore::uInt i=offset;i<tdata.nelements()-offset;i++)
1314 : {
1315 0 : winstart = i-offset;
1316 0 : winend = i+offset;
1317 0 : if(winstart<0)winstart=0;
1318 0 : if(static_cast<casacore::uInt>(winend)>=tdata.nelements()) winend=tdata.nelements()-1;
1319 0 : if(winend <= winstart) break;
1320 0 : winsum=0.0;
1321 0 : for(casacore::uInt wi=winstart;wi<=static_cast<casacore::uInt>(winend);++wi)
1322 0 : winsum += fit[wi];
1323 0 : fit[i] = winsum/(winend-winstart+1);
1324 : }
1325 :
1326 :
1327 : /* Calculate the STD of the fit */
1328 0 : sd = UStd(tdata,tfband,fit);
1329 0 : if(j>=2) TOL=2;
1330 0 : else TOL=3;
1331 :
1332 : /* Display the Fit and the data */
1333 : #ifdef UPLOT
1334 : Plot_ds9(data.nelements(),data,fit1);
1335 : #endif
1336 :
1337 : /* Detect outliers */
1338 0 : for(casacore::uInt i=0;i<tdata.nelements();i++)
1339 : {
1340 0 : if(tdata[i]-fit[i] > TOL*sd)
1341 0 : tfband[i]=true;
1342 : }
1343 :
1344 : } // for j
1345 :
1346 0 : } // end of CleanBand
1347 :
1348 :
1349 :
1350 : /* Fit a polynomial to 'data' from lim1 to lim2, of given degree 'deg',
1351 : * taking care of flags in 'flag', and returning the fitted values in 'fit' */
1352 0 : void RFATimeFreqCrop :: PolyFit(casacore::Vector<casacore::Float> data,casacore::Vector<casacore::Bool> flag, casacore::Vector<casacore::Float> fit, casacore::uInt lim1, casacore::uInt lim2,casacore::uInt deg)
1353 : {
1354 0 : static casacore::Vector<casacore::Double> x;
1355 0 : static casacore::Vector<casacore::Double> y;
1356 0 : static casacore::Vector<casacore::Double> sig;
1357 0 : static casacore::Vector<casacore::Double> solution;
1358 :
1359 0 : casacore::uInt cnt=0;
1360 0 : for(casacore::uInt i=lim1;i<=lim2;i++)
1361 0 : if(flag[i]==false) cnt++;
1362 :
1363 0 : if(cnt <= deg)
1364 : {
1365 0 : LineFit(data,flag,fit,lim1,lim2);
1366 0 : return;
1367 : }
1368 :
1369 :
1370 0 : casacore::LinearFit<casacore::Double> fitter;
1371 0 : casacore::Polynomial<casacore::AutoDiff<casacore::Double> > combination(deg);
1372 :
1373 :
1374 0 : combination.setCoefficient(0,0.0);
1375 0 : if (deg >= 1) combination.setCoefficient(1, 0.0);
1376 0 : if (deg >= 2) combination.setCoefficient(2, 0.0);
1377 0 : if (deg >= 3) combination.setCoefficient(3, 0.0);
1378 0 : if (deg >= 4) combination.setCoefficient(4, 0.0);
1379 :
1380 0 : x.resize(lim2-lim1+1);
1381 0 : y.resize(lim2-lim1+1);
1382 0 : sig.resize(lim2-lim1+1);
1383 0 : solution.resize(deg+1);
1384 :
1385 0 : for(casacore::uInt i=lim1;i<=lim2;i++)
1386 : {
1387 0 : x[i-lim1] = i+1;
1388 0 : y[i-lim1] = data[i];
1389 0 : sig[i-lim1] = (flag[i]==true)?0:1;
1390 : }
1391 :
1392 0 : fitter.asWeight(true);
1393 :
1394 0 : fitter.setFunction(combination);
1395 0 : solution = fitter.fit(x,y,sig);
1396 :
1397 0 : for(casacore::uInt i=lim1;i<=lim2;i++)
1398 : {
1399 0 : fit[i]=0;
1400 0 : for(casacore::uInt j=0;j<deg+1;j++)
1401 0 : fit[i] += solution[j]*pow((double)(x[i-lim1]),(double)j);
1402 : }
1403 :
1404 0 : }
1405 :
1406 :
1407 :
1408 : /* Fit a LINE to 'data' from lim1 to lim2, taking care of flags in
1409 : * 'flag', and returning the fitted values in 'fit' */
1410 0 : void RFATimeFreqCrop :: LineFit(casacore::Vector<casacore::Float> data, casacore::Vector<casacore::Bool> flag, casacore::Vector<casacore::Float> fit, casacore::uInt lim1, casacore::uInt lim2)
1411 : {
1412 0 : float Sx = 0, Sy = 0, Sxx = 0, Sxy = 0, S = 0, a, b, sd, mn;
1413 :
1414 0 : mn = UMean(data, flag);
1415 0 : sd = UStd (data, flag, mn);
1416 :
1417 0 : for (casacore::uInt i = lim1; i <= lim2; i++)
1418 : {
1419 0 : if (flag[i] == false) // if unflagged
1420 : {
1421 0 : S += 1 / (sd * sd);
1422 0 : Sx += i / (sd * sd);
1423 0 : Sy += data[i] / (sd * sd);
1424 0 : Sxx += (i * i) / (sd * sd);
1425 0 : Sxy += (i * data[i]) / (sd * sd);
1426 : }
1427 : }
1428 0 : a = (Sxx * Sy - Sx * Sxy) / (S * Sxx - Sx * Sx);
1429 0 : b = (S * Sxy - Sx * Sy) / (S * Sxx - Sx * Sx);
1430 :
1431 0 : for (casacore::uInt i = lim1; i <= lim2; i++)
1432 0 : fit[i] = a + b * i;
1433 :
1434 0 : }
1435 :
1436 : /* Return antenna numbers from baseline number - upper triangle storage */
1437 0 : void RFATimeFreqCrop :: Ants(casacore::uInt bs, casacore::uInt *a1, casacore::uInt *a2)
1438 : {
1439 0 : casacore::uInt sum=0,cnt=0;
1440 0 : for(casacore::uInt i=(NumAnt);i>1;i--)
1441 : {
1442 0 : sum += i;
1443 0 : if(sum<=bs) cnt++;
1444 0 : else break;
1445 : }
1446 0 : *a1 = cnt;
1447 :
1448 0 : sum = (NumAnt)*((NumAnt)+1)/2 - ((NumAnt)-(*a1))*((NumAnt)-(*a1)+1)/2;
1449 :
1450 0 : *a2 = bs - sum + (*a1);
1451 0 : }
1452 :
1453 : /* Return baseline index from a pair of antenna numbers - upper triangle storage */
1454 0 : casacore::uInt RFATimeFreqCrop :: BaselineIndex(casacore::uInt /*row*/, casacore::uInt a1, casacore::uInt a2)
1455 : {
1456 0 : return ( (NumAnt)*((NumAnt)+1)/2 - ((NumAnt)-a1)*((NumAnt)-a1+1)/2 + (a2 - a1) );
1457 : }
1458 :
1459 :
1460 :
1461 : /* Display a 2D data set on DS9 in gray scale */
1462 0 : void RFATimeFreqCrop :: Display_ds9(casacore::Int xdim, casacore::Int ydim, casacore::Matrix<casacore::Float> &data, casacore::Int frame)
1463 : {
1464 :
1465 0 : FILE *SAOout = NULL;
1466 : char tmp[100];
1467 0 : char server[100] = "ds9";
1468 0 : int bitpix = -32;
1469 0 : char xpa[100] = "";
1470 :
1471 0 : strcpy (xpa, "xpaset");
1472 :
1473 0 : sprintf (tmp, "%s %s \"array [xdim=%d,ydim=%d,bitpix=%d]\"\n",
1474 : xpa, server, xdim, ydim, bitpix);
1475 0 : SAOout = (FILE *) popen (tmp, "w");
1476 :
1477 0 : if (frame > 0)
1478 : {
1479 0 : sprintf (tmp, "echo \"frame %d \" | %s %s \n", frame, xpa, server);
1480 0 : system (tmp);
1481 : }
1482 :
1483 0 : if (SAOout != NULL)
1484 : {
1485 : // for (i = 0; i < ydim; i++)
1486 : // fwrite (data[i], sizeof (float) * xdim, 1, SAOout);
1487 0 : casacore::Bool deleteit=false;
1488 0 : casacore::Float *dataptr = data.getStorage(deleteit);
1489 0 : fwrite(dataptr, sizeof(casacore::Float)*xdim*ydim, 1, SAOout);
1490 0 : pclose (SAOout);
1491 : //if(deleteit) data.freeStorage(dataptr,true);
1492 : }
1493 : else
1494 : {
1495 0 : perror ("Error in opening SAO - ds9 \n");
1496 : }
1497 :
1498 0 : if (frame > 0)
1499 : {
1500 0 : system ("xpaset -p ds9 zoom to fit");
1501 : }
1502 :
1503 : //cout << " Press enter to continue... " << endl;
1504 : //getchar();
1505 :
1506 :
1507 0 : }
1508 :
1509 : /* Display a line plot in DS9 !!! */
1510 0 : void RFATimeFreqCrop :: Plot_ds9(casacore::Int dim, casacore::Vector<casacore::Float> data1, casacore::Vector<casacore::Float> data2)
1511 : {
1512 :
1513 : // FILE *SAOout = NULL;
1514 : //char tmp[100];
1515 : //char server[100] = "ds9";
1516 : //int bitpix = -32;
1517 : int i;
1518 0 : char xpa[100] = "";
1519 :
1520 : //static casacore::Bool firstentry=true;
1521 0 : strcpy (xpa, "xpaset");
1522 :
1523 0 : casacore::String cmd("");
1524 :
1525 : {
1526 0 : cmd = "xpaset -p ds9 plot flagger clear \n";
1527 0 : system(cmd.data());
1528 : // SAOout = (FILE *) popen (tmp, "w");
1529 : }
1530 :
1531 :
1532 0 : cmd = "echo '";
1533 0 : for(i=0;i<dim;i++)
1534 : {
1535 0 : cmd += casacore::String::toString(i) + " " + casacore::String::toString(data2[i]) + " " ;
1536 : }
1537 0 : cmd += "\n' | xpaset ds9 plot flagger data xy\n";
1538 0 : cmd += "xpaset -p ds9 plot flagger color linear blue\n";
1539 0 : cmd += "xpaset -p ds9 plot flagger line linear width 2.0\n";
1540 0 : system(cmd.data());
1541 :
1542 0 : cmd = "echo '";
1543 0 : for(i=0;i<dim;i++)
1544 : {
1545 0 : cmd += casacore::String::toString(i) + " " + casacore::String::toString(data1[i]) + " " ;
1546 : }
1547 0 : cmd += "\n' | xpaset ds9 plot flagger data xy\n";
1548 0 : cmd += "xpaset -p ds9 plot flagger color linear red\n";
1549 0 : cmd += "xpaset -p ds9 plot flagger line linear width 2.0\n";
1550 0 : system(cmd.data());
1551 :
1552 0 : }// end of plot_ds9
1553 :
1554 :
1555 : } //# NAMESPACE CASA - END
1556 :
|
