Line data Source code
1 : //# RFADiffBase.cc: this defines RFADiffBase
2 : //# Copyright (C) 2000,2001
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/RFADiffBase.h>
28 : #include <msvis/MSVis/VisibilityIterator.h>
29 : #include <msvis/MSVis/VisBuffer.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/MaskArrMath.h>
32 : #include <casacore/casa/Arrays/ArrayLogical.h>
33 : #include <casacore/casa/Arrays/Slice.h>
34 :
35 : using namespace casacore;
36 : namespace casa { //# NAMESPACE CASA - BEGIN
37 :
38 : Bool RFADiffBase::dummy_Bool;
39 :
40 : // -----------------------------------------------------------------------
41 : // RFADiffBase constructor
42 : // Construct from a Record of parameters
43 : // -----------------------------------------------------------------------
44 0 : RFADiffBase::RFADiffBase ( RFChunkStats &ch,const RecordInterface &parm ) :
45 : RFAFlagCubeBase(ch,parm),
46 0 : clip_level(parm.asDouble(RF_THR)),
47 0 : row_clip_level(parm.asDouble(RF_ROW_THR)),
48 0 : disable_row_clip(parm.asBool(RF_ROW_DISABLE)),
49 0 : rowclipper(chunk,flag,row_clip_level,parm.asInt(RF_ROW_HW))
50 : {
51 0 : }
52 :
53 : // -----------------------------------------------------------------------
54 : // Destructor
55 : // -----------------------------------------------------------------------
56 0 : RFADiffBase::~RFADiffBase ()
57 : {
58 0 : }
59 :
60 : // -----------------------------------------------------------------------
61 : // RFADiffBase::getDefaults
62 : // Returns record of default parameters
63 : // -----------------------------------------------------------------------
64 0 : const RecordInterface & RFADiffBase::getDefaults ()
65 : {
66 0 : static Record rec;
67 : // create record description on first entry
68 0 : if( !rec.nfields() )
69 : {
70 0 : rec = RFAFlagCubeBase::getDefaults();
71 0 : rec.define(RF_THR,(Double)5);
72 0 : rec.define(RF_ROW_THR,(Double)10);
73 0 : rec.define(RF_ROW_HW,(Int)6);
74 0 : rec.define(RF_ROW_DISABLE,false);
75 0 : rec.setComment(RF_THR,"Pixel flagging threshold, in AADs");
76 0 : rec.setComment(RF_ROW_THR,"Row flagging threshold, in AADs");
77 0 : rec.setComment(RF_ROW_HW,"Row flagging, half-window of sliding median");
78 0 : rec.setComment(RF_ROW_DISABLE,"Disable row-based flagging");
79 : }
80 0 : return rec;
81 : }
82 :
83 : // -----------------------------------------------------------------------
84 : // RFADiffBase::getDesc
85 : // Returns description of parameters
86 : // -----------------------------------------------------------------------
87 0 : String RFADiffBase::getDesc ()
88 : {
89 : char s[128];
90 0 : if( disable_row_clip )
91 0 : sprintf(s,"%s=%.1f",RF_THR,clip_level);
92 : else
93 0 : sprintf(s,"%s=%.1f %s=%.1f",RF_THR,clip_level,RF_ROW_THR,row_clip_level);
94 0 : String str(s);
95 0 : return str+" "+RFAFlagCubeBase::getDesc();
96 : }
97 :
98 0 : uInt RFADiffBase::estimateMemoryUse ()
99 : {
100 0 : return diff.estimateMemoryUse(num(CHAN),num(IFR),num(TIME)) +
101 0 : RFAFlagCubeBase::estimateMemoryUse() + 2;
102 : }
103 :
104 : // -----------------------------------------------------------------------
105 : // RFADiffBase::newChunk
106 : // Sets up for new chunk of data
107 : // -----------------------------------------------------------------------
108 0 : Bool RFADiffBase::newChunk (Int &maxmem)
109 : {
110 : // compute correlations mask
111 0 : corrmask = newCorrMask();
112 0 : if( !corrmask )
113 : {
114 0 : os<<LogIO::WARN<<"missing selected correlations, ignoring this chunk\n"<<LogIO::POST;
115 0 : return active=false;
116 : }
117 : // memory management. Estimate the max memory use for the diff
118 : // lattice, plus a 5% margin
119 0 : Int mmdiff = (Int)(1.05*diff.estimateMemoryUse(num(CHAN),num(IFR),num(TIME)));
120 : // sufficient memory? reserve it
121 0 : if( maxmem>mmdiff )
122 0 : maxmem -= mmdiff;
123 : else // insufficient memory: use disk-based lattice
124 : {
125 0 : mmdiff = 0;
126 0 : maxmem -= 2; // reserve 2Mb for the iterator anyway
127 : }
128 : // init flag cube
129 0 : RFAFlagCubeBase::newChunk(maxmem);
130 : // create a temp lattice to hold nchan x num(IFR) x ntime diff-medians
131 0 : diff.init(num(CHAN),num(IFR),num(TIME),num(CORR), nAgent, mmdiff,2);
132 : // init the row-clipper object
133 0 : rowclipper.init(num(IFR),num(TIME));
134 0 : diffrow.resize(num(CHAN));
135 :
136 : // if rows are too short, there's no point in flagging them in toto
137 : // based on their noise level
138 0 : clipping_rows = !disable_row_clip;
139 0 : if( num(CHAN)<10 )
140 0 : clipping_rows = false;
141 :
142 0 : return active=true;
143 : }
144 :
145 : // -----------------------------------------------------------------------
146 : // RFADiffBase::endChunk
147 : // Resets at end of chunk
148 : // -----------------------------------------------------------------------
149 0 : void RFADiffBase::endChunk ()
150 : {
151 0 : RFAFlagCubeBase::endChunk();
152 0 : diff.cleanup();
153 0 : rowclipper.cleanup();
154 0 : diffrow.resize();
155 0 : }
156 :
157 : // -----------------------------------------------------------------------
158 : // RFADiffBase::startData
159 : // Prepares for an data pass over a VisIter chunk
160 : // -----------------------------------------------------------------------
161 0 : void RFADiffBase::startData (bool verbose)
162 : {
163 0 : RFAFlagCubeBase::startData(verbose);
164 0 : diff.reset(chunk.npass()>0,true);
165 0 : rowclipper.reset();
166 :
167 0 : pflagiter = &flag.iterator();
168 0 : }
169 :
170 : // -----------------------------------------------------------------------
171 : // RFADiffBase::startDry
172 : // Prepares for an dry pass
173 : // -----------------------------------------------------------------------
174 0 : void RFADiffBase::startDry (bool verbose)
175 : {
176 0 : RFAFlagCubeBase::startDry(verbose);
177 0 : diff.reset(chunk.npass()>0,false);
178 0 : rowclipper.reset();
179 0 : pflagiter = &flag.iterator();
180 0 : }
181 :
182 :
183 : // -----------------------------------------------------------------------
184 : // RFADiffBase::iterTime
185 : // Default version of iter time just keeps the diff and flag lattices
186 : // in sync with the time slot.
187 : // -----------------------------------------------------------------------
188 0 : RFA::IterMode RFADiffBase::iterTime (uInt it)
189 : {
190 0 : RFAFlagCubeBase::iterTime(it);
191 0 : diff.advance(it);
192 0 : return RFA::CONT;
193 : }
194 :
195 : // -----------------------------------------------------------------------
196 : // RFADiffBase::iterDry
197 : // Dry run iterator: recomputes the AAD and does flagging
198 : // -----------------------------------------------------------------------
199 0 : RFA::IterMode RFADiffBase::iterDry ( uInt it )
200 : {
201 0 : RFAFlagCubeBase::iterDry(it);
202 0 : diff.advance(it);
203 :
204 0 : for( uInt ifr=0; ifr<num(IFR); ifr++ ) // outer loop over IFRs
205 : {
206 0 : if( flag.rowFlagged(ifr,it) ) // skip if whole row is flagged
207 : {
208 0 : continue;
209 : }
210 0 : Float thr = clip_level*rowclipper.sigma0(ifr,it);
211 0 : idiffrow=0;
212 0 : Bool updated=false;
213 0 : for( uInt ich=0; ich<num(CHAN); ich++ ) // loop over channels
214 : {
215 0 : if( flag.preFlagged(ich,ifr) ) // skip pixel if pre-flagged
216 0 : continue;
217 :
218 0 : Float d = diff(ich,ifr);
219 0 : if( d > thr ) // should be clipped?
220 : {
221 0 : Bool res = flag.setFlag(ich,ifr);
222 0 : updated |= res;
223 : }
224 : else
225 : {
226 0 : diffrow(idiffrow++) = d;
227 : }
228 : } // for(ich)
229 : // update the noise level, if any changes in flags
230 0 : if( updated )
231 0 : rowclipper.setSigma(ifr,it,idiffrow ? median( diffrow(Slice(0,idiffrow)) ) : -1 );
232 : } // for(ifr)
233 0 : return CONT;
234 : }
235 :
236 : // -----------------------------------------------------------------------
237 : // RFADiffBase::endData
238 : // After a data pass, we always request one more dry pass
239 : // -----------------------------------------------------------------------
240 0 : RFA::IterMode RFADiffBase::endData ()
241 : {
242 0 : RFAFlagCubeBase::endData();
243 : uInt dum;
244 0 : rowclipper.updateSigma(dum, dum, true, false);
245 :
246 0 : return RFA::DRY;
247 : }
248 :
249 : // after a dry pass - see if AAD has managed to update itself
250 : // significantly
251 0 : RFA::IterMode RFADiffBase::endDry ()
252 : {
253 0 : RFAFlagCubeBase::endDry();
254 : // update the reference AAD
255 : uInt ifrmax,itmax;
256 0 : Float dmax = rowclipper.updateSigma(ifrmax,itmax, true, false);
257 :
258 0 : dprintf(os,"Max diff (%f) at ifr %d (%s), it %d: new sigma is %f\n",
259 0 : dmax,ifrmax,chunk.ifrString(ifrmax).chars(),itmax,rowclipper.sigma0(ifrmax,itmax));
260 :
261 : // no significant change this pass? Then we're really through with it.
262 0 : if( dmax <= RFA_AAD_CHANGE )
263 0 : return RFA::STOP;
264 : // else try another dry pass
265 : // NB: perhaps request a data pass, if too many flags?
266 0 : return RFA::DRY;
267 : }
268 :
269 0 : void RFADiffBase::startDataRow (uInt)
270 : {
271 0 : idiffrow=0;
272 0 : }
273 :
274 0 : void RFADiffBase::endDataRow (uInt ifr)
275 : {
276 : // if( !idiffrow )
277 : // dprintf(os,"No data points at ifr %d\n",ifr);
278 0 : Float sigma = idiffrow ? median( diffrow( Slice(0,idiffrow) ) ) : -1;
279 0 : uInt it = diff.position();
280 0 : rowclipper.setSigma(ifr,it,sigma);
281 0 : }
282 :
283 : // -----------------------------------------------------------------------
284 : // RFADiffBase::setDiff
285 : // Meant to be called during a data pass. Updates the difference
286 : // lattice; flags things if appropriate.
287 : // This assumes caller has already checked existing pre-flags,
288 : // and that the current data point is NOT flagged.
289 : // Returns the threshold used for flagging, or 0 if no threshold is
290 : // yet available.
291 : // -----------------------------------------------------------------------
292 0 : Float RFADiffBase::setDiff ( uInt ich,uInt ifr,Float d,Bool &flagged )
293 : {
294 0 : Float thr = 0;
295 0 : d = abs(d);
296 0 : uInt it = diff.position();
297 :
298 : // write diff to lattice
299 0 : diff(ich,ifr) = d;
300 :
301 0 : flagged=false;
302 0 : if( chunk.npass() && rowclipper.sigma0(ifr,it)>0 )
303 : {
304 0 : thr = rowclipper.sigma0(ifr,it);
305 :
306 0 : if( d > thr )
307 : {
308 0 : if( flag.setFlag(ich,ifr,*pflagiter) )
309 0 : rowclipper.markSigma(ifr);
310 0 : flagged=true;
311 : }
312 : }
313 0 : if( !flagged )
314 0 : diffrow(idiffrow++) = d;
315 :
316 0 : return thr;
317 : }
318 :
319 : // -----------------------------------------------------------------------
320 : // class RFADiffMapBase
321 : //
322 : //
323 : //
324 : // -----------------------------------------------------------------------
325 :
326 : // -----------------------------------------------------------------------
327 : // RFADiffMapBase constructor
328 : // Construct from Record of parameters
329 : // -----------------------------------------------------------------------
330 0 : RFADiffMapBase::RFADiffMapBase ( RFChunkStats &ch,const RecordInterface &parm )
331 : : RFADiffBase(ch,parm),
332 0 : RFDataMapper(parm.asArrayString(RF_EXPR),parm.asString(RF_COLUMN))
333 : {
334 0 : }
335 :
336 : // -----------------------------------------------------------------------
337 : // RFADiffMapBase destructor
338 : // -----------------------------------------------------------------------
339 0 : RFADiffMapBase::~RFADiffMapBase ()
340 : {
341 0 : }
342 :
343 : // -----------------------------------------------------------------------
344 : // RFADiffMapBase::getDefaults
345 : // Returns record of default paramaters
346 : // -----------------------------------------------------------------------
347 0 : const RecordInterface & RFADiffMapBase::getDefaults ()
348 : {
349 0 : static Record rec;
350 : // create record description on first entry
351 0 : if( !rec.nfields() )
352 : {
353 0 : rec = RFADiffBase::getDefaults();
354 0 : rec.define(RF_COLUMN,"DATA");
355 0 : rec.define(RF_EXPR,"+ ABS XX YY");
356 0 : rec.setComment(RF_COLUMN,"Use column: [DATA|MODEL|CORRected]");
357 0 : rec.setComment(RF_EXPR,"Expression for deriving value (e.g. \"ABS XX\", \"+ ABS XX YY\")");
358 : }
359 0 : return rec;
360 : }
361 :
362 : // -----------------------------------------------------------------------
363 : // RFADiffMapBase::iterTime
364 : // Sets up mapper
365 : // -----------------------------------------------------------------------
366 0 : RFA::IterMode RFADiffMapBase::iterTime (uInt it)
367 : {
368 0 : setupMapper();
369 0 : return RFADiffBase::iterTime(it);
370 : }
371 :
372 : // -----------------------------------------------------------------------
373 : // RFADiffMapBase::getDesc
374 : // Returns short description of parameters
375 : // -----------------------------------------------------------------------
376 0 : String RFADiffMapBase::getDesc ()
377 : {
378 0 : return RFDataMapper::description()+"; "+RFADiffBase::getDesc();
379 : }
380 :
381 : } //# NAMESPACE CASA - END
382 :
|
