Line data Source code
1 : //# RFAMedianClip.cc: this defines RFAMedianClip
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/RFAMedianClip.h>
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 :
31 : using namespace casacore;
32 : namespace casa { //# NAMESPACE CASA - BEGIN
33 :
34 : // -----------------------------------------------------------------------
35 : // RFATimeMedian
36 : // Accumulator class for computing sliding median per channels over time.
37 : // Internally, we store a cubic TempLattice of ntime x nifr x nchan
38 : // medians.
39 : // -----------------------------------------------------------------------
40 0 : RFATimeMedian::RFATimeMedian( RFChunkStats &ch,const RecordInterface &parm ) :
41 0 : RFADiffMapBase(ch,parm)
42 : {
43 0 : halfwin = (uInt)parm.asInt(RF_HW);
44 0 : msl = NULL;
45 0 : if( fieldType(parm,RF_DEBUG,TpArrayInt) )
46 : {
47 0 : Vector<Int> dbg;
48 0 : parm.get(RF_DEBUG,dbg);
49 0 : if (dbg.nelements() != 2 && dbg.nelements() != 3) {
50 0 : os<<"\""<<RF_DEBUG<<"\" parameter must be [NCH,NIFR] or [NCH,ANT1,ANT2]"<<LogIO::EXCEPTION;
51 : }
52 0 : }
53 0 : }
54 :
55 0 : RFATimeMedian::~RFATimeMedian ()
56 : {
57 0 : if( msl ) delete [] msl;
58 0 : }
59 :
60 0 : const RecordInterface & RFATimeMedian::getDefaults ()
61 : {
62 0 : static Record rec;
63 : // create record description on first entry
64 0 : if( !rec.nfields() )
65 : {
66 0 : rec = RFADiffMapBase::getDefaults();
67 0 : rec.define(RF_NAME,"TimeMedian");
68 0 : rec.define(RF_HW,10);
69 0 : rec.define(RF_DEBUG,false);
70 0 : rec.setComment(RF_HW,"Sliding window half-width");
71 0 : rec.setComment(RF_DEBUG,"Set to [CHAN,IFR] to produce debugging plots");
72 : }
73 0 : return rec;
74 : }
75 :
76 : // resets for new calculation
77 0 : Bool RFATimeMedian::newChunk (Int &maxmem)
78 : {
79 0 : if( num(TIME) < halfwin*4 )
80 : {
81 0 : os<<LogIO::WARN<<name()<<
82 0 : ": too few (" << num(TIME) << ") time slots (" << halfwin*4 << " needed), ignoring this chunk\n"<<LogIO::POST;
83 0 : return active=false;
84 : }
85 0 : maxmem -= 2;
86 : // reserve memory for our bunch of median sliders
87 0 : maxmem -= (num(CHAN)*num(IFR)*MedianSlider::objsize(halfwin))/(1024*1024)+1;
88 : // call parent's newChunk
89 0 : if( !RFADiffMapBase::newChunk(maxmem) )
90 0 : return active=false;
91 : // create local flag iterator
92 0 : flag_iter = flag.newCustomIter();
93 0 : pflagiter = &flag_iter;
94 0 : return active=true;
95 : }
96 :
97 0 : void RFATimeMedian::endChunk ()
98 : {
99 0 : RFADiffMapBase::endChunk();
100 : // create local flag iterator
101 0 : flag_iter = FlagCubeIterator();
102 0 : if( msl ) delete [] msl;
103 0 : msl = NULL;
104 0 : }
105 :
106 : // startData
107 : // create new median sliders at start of data pass
108 0 : void RFATimeMedian::startData (bool verbose)
109 : {
110 0 : RFADiffMapBase::startData(verbose);
111 0 : flag_iter.reset();
112 0 : if( msl ) delete [] msl;
113 : // this is a workaround for a compiler bug that we occasionally see
114 0 : uInt tmpnum2 = num(CHAN)*num(IFR);
115 : // create nchan x nifr median sliders
116 0 : msl = new MedianSlider[tmpnum2];
117 0 : for(uInt i=0; i<num(CHAN)*num(IFR); i++)
118 0 : msl[i] = MedianSlider(halfwin);
119 0 : }
120 :
121 : // iterTime
122 : // During data passes, keep the flag lattice lagging a halw-window behind
123 0 : RFA::IterMode RFATimeMedian::iterTime ( uInt it )
124 : {
125 : // gets pointer to visibilities cube
126 0 : setupMapper();
127 : // Advance sync flag iterator
128 0 : flag.advance(it,true);
129 : // During a data pass, keep the diff-lattice iterator lagging a half-window behind
130 : // and also maintain a custom flag iterator
131 0 : if( it >= halfwin )
132 : {
133 0 : diff.advance(it-halfwin);
134 0 : flag_iter.advance(it-halfwin);
135 : }
136 0 : return RFA::CONT;
137 : }
138 :
139 : // -----------------------------------------------------------------------
140 : // RFATimeMedian::iterRow
141 : // Processes one row of data for per-channel medians over time
142 : // -----------------------------------------------------------------------
143 0 : RFA::IterMode RFATimeMedian::iterRow ( uInt irow )
144 : {
145 : // start filling deviations when we get to time slot HW
146 0 : uInt iifr = chunk.ifrNum(irow);
147 0 : uInt it = chunk.iTime();
148 0 : Bool fill = ( chunk.iTime() >= (Int)halfwin );
149 0 : Bool rowfl = chunk.npass() ? flag.rowFlagged(iifr,it)
150 0 : : flag.rowPreFlagged(iifr,it);
151 0 : if( rowfl )
152 : {
153 : // the whole row is flagged, so just advance all median sliders
154 0 : for( uInt i=0; i<num(CHAN); i++ )
155 0 : slider(i,iifr).next();
156 : }
157 : else
158 : {
159 0 : startDataRow(iifr);
160 : // loop over channels for this spw, ifr
161 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
162 : {
163 : // get derived flags and values
164 0 : Float val = 0;
165 : // during first pass, look at pre-flags only. During subsequent passes,
166 : // look at all flags
167 0 : Bool fl = chunk.npass() ? flag.anyFlagged(ich,iifr) : flag.preFlagged(ich,iifr);
168 0 : if( fl )
169 : {
170 : }
171 : else
172 : {
173 0 : val = mapValue(ich,irow);
174 : }
175 0 : slider(ich,iifr).add( val,fl );
176 : // are we filling in the diff-median lattice already?
177 0 : if( fill )
178 : {
179 0 : Float d = slider(ich,iifr).diff(fl);
180 0 : if( !fl ) // ignore if flagged
181 0 : setDiff( ich,iifr,d );
182 : }
183 : }
184 0 : endDataRow(iifr);
185 : }
186 0 : return RFA::CONT;
187 : }
188 :
189 : // -----------------------------------------------------------------------
190 : // RFATimeMedian::endData
191 : // Called at end of iteration - fill remaining slots in the diff-lattice
192 : // -----------------------------------------------------------------------
193 0 : RFA::IterMode RFATimeMedian::endData ()
194 : {
195 0 : for( uInt it=num(TIME)-halfwin; it<num(TIME); it++ )
196 : {
197 0 : diff.advance(it);
198 0 : for( uInt i = 0; i<num(IFR); i++ )
199 : {
200 0 : startDataRow(i);
201 0 : for( uInt j = 0; j<num(CHAN); j++ )
202 : {
203 0 : slider(j,i).next();
204 : Bool fl;
205 0 : Float diff = slider(j,i).diff(fl);
206 0 : if( !fl )
207 0 : setDiff(j,i,diff);
208 :
209 : }
210 0 : endDataRow(i);
211 : }
212 : }
213 : // destroy sliders
214 0 : if( msl )
215 : {
216 0 : delete [] msl;
217 0 : msl = NULL;
218 : }
219 0 : return RFADiffMapBase::endData();
220 : }
221 :
222 :
223 : // -----------------------------------------------------------------------
224 : // getDesc
225 : // -----------------------------------------------------------------------
226 0 : String RFATimeMedian::getDesc ()
227 : {
228 : char s[128];
229 0 : sprintf(s," %s=%d",RF_HW,halfwin);
230 0 : return RFADiffMapBase::getDesc()+s;
231 : }
232 :
233 :
234 : // -----------------------------------------------------------------------
235 : // RFAFreqMedian
236 : // Accumulator class for computing sliding median per time slot over
237 : // channels.
238 : // -----------------------------------------------------------------------
239 0 : RFAFreqMedian::RFAFreqMedian( RFChunkStats &ch,const RecordInterface &parm ) :
240 0 : RFADiffMapBase(ch,parm)
241 : {
242 0 : halfwin = (uInt)parm.asInt(RF_HW);
243 0 : if( parm.isDefined(RF_DEBUG) && parm.dataType(RF_DEBUG) == TpArrayInt )
244 : {
245 0 : Vector<Int> dbg;
246 0 : parm.get(RF_DEBUG,dbg);
247 0 : if (dbg.nelements() != 2 && dbg.nelements() != 3)
248 : {
249 0 : os<<"\""<<RF_DEBUG<<"\" parameter must be [NIFR,NTIME] or [ANT1,ANT2,NTIME]"<<LogIO::EXCEPTION;
250 : }
251 0 : }
252 0 : }
253 :
254 :
255 : // -----------------------------------------------------------------------
256 : // RFAFreqMedian::newChunk
257 : //
258 : // -----------------------------------------------------------------------
259 0 : Bool RFAFreqMedian::newChunk (Int &maxmem)
260 : {
261 0 : if( num(CHAN) < halfwin*4 )
262 : {
263 0 : os<<LogIO::WARN<<name()<<": too few channels, ignoring this chunk\n"<<LogIO::POST;
264 0 : return active=false;
265 : }
266 0 : return active=RFADiffMapBase::newChunk(maxmem);
267 : }
268 :
269 : // -----------------------------------------------------------------------
270 : // RFAFreqMedian::iterRow
271 : // Processes one row of data for median over frequency
272 : // -----------------------------------------------------------------------
273 0 : RFA::IterMode RFAFreqMedian::iterRow ( uInt irow )
274 : {
275 : // during first pass, compute diff-median. Also keep track of the AAD.
276 0 : uInt iifr = chunk.ifrNum(irow);
277 0 : uInt it = chunk.iTime();
278 :
279 0 : Bool rowfl = chunk.npass() ? flag.rowFlagged(iifr,it)
280 0 : : flag.rowPreFlagged(iifr,it);
281 0 : if( rowfl )
282 : {
283 : }
284 : else // row not flagged
285 : {
286 0 : startDataRow(iifr);
287 0 : MedianSlider msl(halfwin);
288 : // loop through all channels in this window
289 0 : for( uInt i = 0; i<num(CHAN); i++ )
290 : {
291 0 : Float val = 0;
292 : // during first pass, look at pre-flags only. During subsequent passes,
293 : // look at all flags
294 0 : Bool fl = chunk.npass() ? flag.anyFlagged(i,iifr) : flag.preFlagged(i,iifr);
295 0 : if( fl )
296 : {
297 : }
298 : else
299 : {
300 0 : val = mapValue(i,irow);
301 : }
302 0 : msl.add( val,fl );
303 0 : if( i>=halfwin )
304 : {
305 0 : Float d = msl.diff(fl);
306 0 : if( !fl )
307 0 : setDiff(i-halfwin,iifr,d);
308 : }
309 : }
310 : // finish sliding the medians for remaining channels
311 0 : for( uInt i=num(CHAN)-halfwin; i<num(CHAN); i++ )
312 : {
313 0 : msl.next();
314 : Bool fl;
315 0 : Float d = msl.diff(fl);
316 0 : if( !fl )
317 0 : setDiff(i,iifr,d);
318 : }
319 0 : endDataRow(iifr);
320 0 : }
321 0 : return RFA::CONT;
322 : }
323 :
324 : // -----------------------------------------------------------------------
325 : // getDesc
326 : // -----------------------------------------------------------------------
327 0 : String RFAFreqMedian::getDesc ()
328 : {
329 : char s[128];
330 0 : sprintf(s," %s=%d",RF_HW,halfwin);
331 0 : return RFADiffMapBase::getDesc()+s;
332 : }
333 :
334 0 : const RecordInterface & RFAFreqMedian::getDefaults ()
335 : {
336 0 : static Record rec;
337 : // create record description on first entry
338 0 : if( !rec.nfields() )
339 : {
340 0 : rec = RFATimeMedian::getDefaults();
341 0 : rec.define(RF_NAME,"FreqMedian");
342 0 : rec.setComment(RF_DEBUG,"Set to [IFR,ITIME] to produce debugging plots");
343 : }
344 0 : return rec;
345 : }
346 :
347 : } //# NAMESPACE CASA - END
348 :
|
