Line data Source code
1 : //# RFANewMedianClip.cc: this defines RFANewMedianClip
2 : //# Copyright (C) 2000,2001,2002,2003,2004
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 :
28 : #include <flagging/Flagging/RFANewMedianClip.h>
29 : #include <casacore/casa/Arrays/ArrayMath.h>
30 : #include <msvis/MSVis/VisBuffer.h>
31 : #include <casacore/casa/System/PGPlotterInterface.h>
32 :
33 : #include <stdio.h>
34 : using namespace casacore;
35 : namespace casa { //# NAMESPACE CASA - BEGIN
36 :
37 : // -----------------------------------------------------------------------
38 : // RFANewMedianClip
39 : // Accumulator class for computing the median per channels over time.
40 : // Internally, we store a matrix of nifr x nchan medians.
41 : // -----------------------------------------------------------------------
42 0 : RFANewMedianClip::RFANewMedianClip( RFChunkStats &ch,const RecordInterface &parm ) :
43 : RFAFlagCubeBase(ch,parm),
44 : RFDataMapper(parm.asArrayString(RF_EXPR),parm.asString(RF_COLUMN)),
45 0 : threshold( parm.asDouble(RF_THR) )
46 : {
47 0 : msl = NULL;
48 0 : }
49 :
50 0 : RFANewMedianClip::~RFANewMedianClip ()
51 : {
52 0 : if( msl ) delete [] msl;
53 0 : }
54 :
55 0 : uInt RFANewMedianClip::estimateMemoryUse ()
56 : {
57 0 : return RFAFlagCubeBase::estimateMemoryUse() +
58 0 : evalue.estimateMemoryUse(num(CHAN),num(IFR),num(TIME)) + 2;
59 : }
60 :
61 :
62 0 : const RecordInterface & RFANewMedianClip::getDefaults ()
63 : {
64 0 : static Record rec;
65 : // create record description on first entry
66 0 : if( !rec.nfields() )
67 : {
68 0 : rec = RFAFlagCubeBase::getDefaults();
69 0 : rec.define(RF_NAME,"NewTimeMedian");
70 0 : rec.define(RF_COLUMN,"DATA");
71 0 : rec.define(RF_EXPR,"+ ABS XX YY");
72 0 : rec.define(RF_THR,Double(5));
73 : // rec.define(RF_DEBUG,false);
74 0 : rec.setComment(RF_COLUMN,"Use column: [DATA|MODEL|CORRected]");
75 0 : rec.setComment(RF_EXPR,"Expression for deriving value (e.g. \"ABS XX\", \"+ ABS XX YY\")");
76 0 : rec.setComment(RF_THR,"Probability cut-off");
77 : // rec.setComment(RF_DEBUG,"Set to [CHAN,IFR] to produce debugging plots");
78 : }
79 0 : return rec;
80 : }
81 :
82 : // resets for new calculation
83 0 : Bool RFANewMedianClip::newChunk (Int &maxmem)
84 : {
85 : // if disk-based flag cube, reserve 2MB for local iterator
86 : // compute correlations mask, return false if fails
87 0 : corrmask = RFDataMapper::corrMask(chunk.visIter());
88 0 : if( !corrmask )
89 : {
90 0 : os<<LogIO::WARN<<"missing selected correlations, ignoring this chunk\n"<<LogIO::POST;
91 0 : return active=false;
92 : }
93 :
94 0 : maxmem -= 2;
95 0 : uInt halfwin = num(TIME)/2;
96 : // reserve memory for our bunch of median sliders
97 0 : maxmem -= (num(CHAN)*num(IFR)*MedianSlider::objsize(halfwin))/(1024*1024)+1;
98 :
99 0 : Int mmdiff = (Int)(1.05*evalue.estimateMemoryUse(num(CHAN),num(IFR),num(TIME))); // sufficient memory? reserve it
100 0 : if( maxmem>mmdiff )
101 0 : maxmem -= mmdiff;
102 : else // insufficient memory: use disk-based lattice
103 : {
104 0 : mmdiff = 0;
105 0 : maxmem -= 2; // reserve 2Mb for the iterator anyway
106 : }
107 :
108 : // call parent's newChunk
109 0 : if( !RFAFlagCubeBase::newChunk(maxmem) )
110 0 : return active=false;
111 :
112 : // create temp lattice for evalues
113 0 : evalue.init(num(CHAN),num(IFR),num(TIME), num(CORR), nAgent, 0, mmdiff ,2);
114 : //init stdev matrix
115 0 : stdev.resize(num(CHAN), num(IFR));
116 0 : stdev.set(0);
117 0 : stdeved = false;
118 : // create local flag iterator
119 0 : flag_iter = flag.newCustomIter();
120 0 : pflagiter = &flag_iter;
121 0 : RFAFlagCubeBase::newChunk(maxmem-=1);
122 :
123 0 : return active=true;
124 : }
125 :
126 0 : void RFANewMedianClip::endChunk ()
127 : {
128 0 : RFAFlagCubeBase::endChunk();
129 : // create local flag iterator
130 0 : flag_iter = FlagCubeIterator();
131 0 : if( msl ) delete [] msl;
132 0 : msl = NULL;
133 0 : }
134 :
135 : // startData
136 : // create new median sliders at start of data pass
137 0 : void RFANewMedianClip::startData (bool verbose)
138 : {
139 : //new added
140 0 : evalue.reset(false,true);
141 :
142 0 : RFAFlagCubeBase::startData(verbose);
143 0 : flag_iter.reset();
144 :
145 0 : pflagiter = &flag.iterator();
146 0 : if( msl ) delete [] msl;
147 0 : globalsigma = 0;
148 : // this is a workaround for a compiler bug that we occasionally see
149 0 : uInt tmpnum2 = num(CHAN)*num(IFR);
150 0 : uInt halfwin = num(TIME)/2;
151 : // create nchan x nifr median sliders
152 0 : msl = new MedianSlider[tmpnum2];
153 0 : for(uInt i=0; i<num(CHAN)*num(IFR); i++)
154 0 : msl[i] = MedianSlider(halfwin);
155 0 : globalmed = MedianSlider(tmpnum2);
156 0 : }
157 :
158 : // iterTime
159 0 : RFA::IterMode RFANewMedianClip::iterTime ( uInt it )
160 : {
161 0 : evalue.advance(it);
162 0 : RFAFlagCubeBase::iterTime(it);
163 : // gets pointer to visibilities cube
164 0 : RFDataMapper::setVisBuffer(chunk.visBuf());
165 : // Advance sync flag iterator
166 0 : flag.advance(it,true);
167 :
168 0 : return RFA::CONT;
169 : }
170 :
171 : // -----------------------------------------------------------------------
172 : // RFANewMedianClip::iterRow
173 : // Processes one row of data for per-channel medians over time
174 : // -----------------------------------------------------------------------
175 0 : RFA::IterMode RFANewMedianClip::iterRow ( uInt irow )
176 : {
177 0 : uInt iifr = chunk.ifrNum(irow);
178 0 : uInt it = chunk.iTime();
179 0 : Float val = 0;
180 :
181 0 : Bool rowfl = chunk.npass() ? flag.rowFlagged(iifr,it)
182 0 : : flag.rowPreFlagged(iifr,it);
183 0 : if( rowfl ) {
184 : // the whole row is flagged, so just advance all median sliders
185 0 : for( uInt i=0; i<num(CHAN); i++ )
186 0 : slider(i,iifr).next();
187 : }
188 : else {
189 : // loop over channels for this spw, ifr
190 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
191 : {
192 : // during first pass, look at pre-flags only. During subsequent passes,
193 : // look at all flags
194 0 : Bool fl = chunk.npass() ? flag.anyFlagged(ich,iifr) : flag.preFlagged(ich,iifr);
195 0 : val = mapValue(ich,irow);
196 0 : slider(ich,iifr).add( val,fl );
197 0 : if( !fl ) {
198 0 : evalue(ich,iifr) = val;
199 : }
200 : }
201 : }
202 0 : return RFA::CONT;
203 : }
204 :
205 :
206 : // -----------------------------------------------------------------------
207 : // RFANewMedianClip::endData
208 : // Called at end of iteration
209 : // -----------------------------------------------------------------------
210 0 : RFA::IterMode RFANewMedianClip::endData ()
211 : {
212 0 : RFAFlagCubeBase::endData();
213 0 : return RFA::DRY;
214 : }
215 :
216 :
217 0 : void RFANewMedianClip::startDry (bool verbose)
218 : {
219 0 : if(!stdeved)
220 0 : RFAFlagCubeBase::startDry(verbose);
221 : // reset lattices to read-only
222 0 : evalue.reset(true,false);
223 0 : pflagiter = &flag.iterator();
224 0 : }
225 :
226 :
227 :
228 : // -----------------------------------------------------------------------
229 : // RFANewMedianClip::iterDry
230 : // Dry run iterator: recomputes the AAD and does flagging
231 : // -----------------------------------------------------------------------
232 0 : RFA::IterMode RFANewMedianClip::iterDry ( uInt it )
233 : {
234 0 : RFAFlagCubeBase::iterDry(it);
235 0 : evalue.advance(it);
236 0 : Float m = globalmed.median();
237 : // already got standard deviation
238 0 : if(stdeved) {
239 0 : Float upperdiff = 0;
240 0 : Float bottomdiff = 0;
241 0 : Float thr = 0;
242 0 : Bool asymmetry = false;
243 0 : for( uInt ifr=0; ifr<num(IFR); ifr++ ) // outer loop over IFRs
244 : {
245 0 : for( uInt ich=0; ich<num(CHAN); ich++ ) // loop over channels
246 : {
247 : //thr = threshold * stdev(ich, ifr);
248 : // globalsigma is the average of standard deviations
249 0 : thr = threshold * globalsigma;
250 0 : if( !flag.anyFlagged(ich, ifr) ) // skip if whole row is flagged
251 : {
252 0 : Float d = evalue(ich,ifr);
253 : // Float m = slider(ich,ifr).median();
254 0 : if( abs(d - m) > thr ) // should be clipped?
255 : {
256 0 : flag.setFlag(ich,ifr, *pflagiter);
257 : }
258 : else {
259 0 : if(d-m > 0 && d - m > upperdiff)
260 0 : upperdiff = d-m;
261 0 : else if( d - m < 0 && d - m < bottomdiff)
262 0 : bottomdiff = d - m;
263 : // cout << evalue(ich, ifr) << endl;
264 : }
265 : }
266 : } // for(ich)
267 : } // for(ifr)
268 0 : if(abs(bottomdiff) > 1.2 * upperdiff || upperdiff > 1.2 * abs(bottomdiff))
269 0 : asymmetry = true;
270 0 : if(asymmetry) {
271 : // cout << " flag the asymmetry data" << endl;
272 0 : for( uInt ifr=0; ifr<num(IFR); ifr++ ) // outer loop over IFRs
273 : {
274 0 : for( uInt ich=0; ich<num(CHAN); ich++ ) // loop over channels
275 : {
276 : // Float thr = threshold * stdev(ich, ifr);
277 : // try global sigma
278 0 : if(upperdiff < abs(bottomdiff))
279 0 : thr = upperdiff;
280 : else {
281 0 : thr = abs(bottomdiff);
282 : }
283 0 : if( !flag.anyFlagged(ich, ifr) ) // skip if whole row is flagged
284 : {
285 0 : Float d = evalue(ich,ifr);
286 0 : if( abs(d - m) > thr ) // should be clipped?
287 : {
288 0 : flag.setFlag(ich,ifr, *pflagiter);
289 : }
290 : }
291 : } // for(ich)
292 : } // for(ifr)
293 : }
294 : } else { //compute the standard deviation
295 0 : for( uInt ifr=0; ifr<num(IFR); ifr++ ) // outer loop over IFRs
296 : {
297 0 : for( uInt ich=0; ich<num(CHAN); ich++ ) // loop over channels
298 : {
299 0 : Bool fl = flag.anyFlagged(ich, ifr);
300 0 : if(!fl) {
301 0 : Float diff = evalue(ich,ifr) - slider(ich,ifr).median();
302 0 : stdev(ich, ifr) += diff * diff;
303 : }
304 : }
305 : }
306 : } // end else
307 0 : return RFA::CONT;
308 : }
309 :
310 0 : RFA::IterMode RFANewMedianClip::endDry ()
311 : {
312 0 : if(stdeved) {
313 : // cout << " early return " << endl;
314 0 : return RFA::STOP;
315 : }
316 0 : Bool dummy = false;
317 0 : stdeved = true;
318 0 : for( uInt ifr=0; ifr<num(IFR); ifr++ ) // outer loop over IFRs
319 : {
320 0 : for( uInt ich=0; ich<num(CHAN); ich++ ) // loop over channels
321 : {
322 0 : if(slider(ich, ifr).nval()){
323 0 : stdev(ich, ifr) /= slider(ich, ifr).nval();
324 : // cout << "variance " << stdev(ich, ifr) << endl;
325 0 : stdev(ich, ifr) = sqrt(stdev(ich, ifr));
326 0 : globalmed.add(slider(ich, ifr).median(), dummy);
327 : } else {
328 0 : stdev(ich, ifr) = 0;
329 : }
330 : }
331 : }
332 0 : globalsigma = sum(stdev)/(num(CHAN) * num(IFR));
333 0 : return RFA::DRY;
334 : }
335 :
336 : // -----------------------------------------------------------------------
337 : // getDesc
338 : // -----------------------------------------------------------------------
339 0 : String RFANewMedianClip::getDesc ()
340 : {
341 0 : String desc( RFDataMapper::description()+"; " );
342 : char s[256];
343 0 : sprintf(s," %s=%f",RF_THR, threshold);
344 0 : desc += s;
345 0 : return RFAFlagCubeBase::getDesc()+s;
346 : }
347 :
348 :
349 : } //# NAMESPACE CASA - END
350 :
|
