Line data Source code
1 : //# RFAUVBinner.cc: this defines RFAUVBinner
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/RFAUVBinner.h>
28 : #include <casacore/casa/BasicMath/Math.h>
29 : #include <casacore/casa/BasicSL/Constants.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 : #include <msvis/MSVis/VisBuffer.h>
35 : #include <casacore/casa/System/PGPlotterInterface.h>
36 :
37 : using namespace casacore;
38 : namespace casa { //# NAMESPACE CASA - BEGIN
39 :
40 0 : RFAUVBinner::RFAUVBinner ( RFChunkStats &ch,const RecordInterface &parm ) :
41 : RFAFlagCubeBase(ch,parm),
42 : RFDataMapper(parm.asArrayString(RF_EXPR),parm.asString(RF_COLUMN)),
43 0 : threshold( parm.asDouble(RF_THR) ),
44 0 : min_population( parm.asInt(RF_MINPOP) )
45 : // rowclipper(chunk,flag,threshold,halfwin)
46 : {
47 : // get bin size arguments
48 0 : if( isFieldSet(parm,RF_NBINS) )
49 : {
50 0 : if( fieldType(parm,RF_NBINS,TpArrayInt) )
51 : {
52 0 : Vector<Int> binsize( parm.asArrayInt(RF_NBINS) );
53 0 : nbin_uv = binsize(0);
54 0 : nbin_y = binsize(1);
55 : }
56 0 : else if( fieldType(parm,RF_NBINS,TpInt) )
57 : {
58 0 : nbin_uv = nbin_y = parm.asInt(RF_NBINS);
59 : }
60 : }
61 : // check threshold for validity
62 0 : if( threshold >= 1 )
63 0 : os<<String("RFAUVBinner: ")+RF_THR+" must be below 1"<<endl<<LogIO::EXCEPTION;
64 0 : if( threshold==0 && !min_population )
65 0 : os<<String("RFAUVBinner: ")+RF_THR+" and/or "+RF_MINPOP+" must be specified"<<endl<<LogIO::EXCEPTION;
66 : // check if a report is requested for a specific channel
67 0 : }
68 :
69 0 : uInt RFAUVBinner::estimateMemoryUse ()
70 : {
71 0 : return RFAFlagCubeBase::estimateMemoryUse() +
72 0 : yvalue.estimateMemoryUse(num(CHAN),num(IFR),num(TIME)) +
73 0 : num(IFR)*num(TIME)*sizeof(Float)/(1024*1024) +
74 0 : nbin_uv*nbin_y*num(CHAN)*sizeof(Int)/(1024*1024);
75 : }
76 :
77 0 : Bool RFAUVBinner::newChunk (Int &maxmem)
78 : {
79 : // compute correlations mask, return false if fails
80 0 : corrmask = RFDataMapper::corrMask(chunk.visIter());
81 0 : if( !corrmask )
82 : {
83 0 : os<<LogIO::WARN<<"missing selected correlations, ignoring this chunk\n"<<LogIO::POST;
84 0 : return active=false;
85 : }
86 : // memory management.
87 : // bin counts are always in memory
88 0 : maxmem -= nbin_uv*nbin_y*num(CHAN)*sizeof(Int)/(1024*1024) +
89 0 : num(IFR)*num(TIME)*sizeof(Float)/(1024*1024);
90 : // Estimate the max memory use for the lattices, plus a 5% margin
91 0 : Int mmdiff = (Int)(1.05*yvalue.estimateMemoryUse(num(CHAN),num(IFR),num(TIME)));
92 : // sufficient memory? reserve it
93 0 : if( maxmem>mmdiff )
94 0 : maxmem -= mmdiff;
95 : else // insufficient memory: use disk-based lattice
96 : {
97 0 : mmdiff = 0;
98 0 : maxmem -= 2; // reserve 2Mb for the iterator anyway
99 : }
100 : // init flag cube
101 0 : RFAFlagCubeBase::newChunk(maxmem);
102 : // create temp lattice for yvalues
103 0 : yvalue.init(num(CHAN),num(IFR),num(TIME),num(CORR),nAgent,mmdiff,2);
104 : // init uvdist matrix
105 0 : uvdist.resize(num(IFR),num(TIME));
106 0 : uvdist.set(-1);
107 : // init min/max estimates
108 0 : ymin.resize(num(CHAN));
109 0 : ymax.resize(num(CHAN));
110 0 : ymin.set(C::flt_max);
111 0 : ymax.set(C::flt_min);
112 0 : uvmin.resize(num(CHAN));
113 0 : uvmax.resize(num(CHAN));
114 0 : uvmin.set(C::flt_max);
115 0 : uvmax.set(0);
116 0 : binned = false;
117 : // finish with init
118 0 : RFAFlagCubeBase::newChunk(maxmem-=1);
119 :
120 0 : return active=true;
121 : }
122 :
123 0 : void RFAUVBinner::endChunk ()
124 : {
125 0 : RFAFlagCubeBase::endChunk();
126 0 : yvalue.cleanup();
127 0 : uvdist.resize();
128 0 : bincounts.resize();
129 0 : ymin.resize();
130 0 : ymax.resize();
131 0 : ybinsize.resize();
132 0 : uvmin.resize();
133 0 : uvmax.resize();
134 0 : uvbinsize.resize();
135 0 : totcounts.resize();
136 : // rowclipper.cleanup();
137 0 : }
138 :
139 0 : void RFAUVBinner::startData (bool verbose)
140 : {
141 : // reset lattices to write-only
142 0 : yvalue.reset(false,true);
143 0 : RFAFlagCubeBase::startData(verbose);
144 : // rowclipper.reset();
145 0 : }
146 :
147 0 : RFA::IterMode RFAUVBinner::iterTime (uInt it)
148 : {
149 0 : yvalue.advance(it);
150 0 : RFAFlagCubeBase::iterTime(it);
151 0 : RFDataMapper::setVisBuffer(chunk.visBuf());
152 : // get UVW data from VisBuffer
153 0 : puvw = & chunk.visBuf().uvw();
154 0 : return RFA::CONT;
155 : }
156 :
157 0 : RFA::IterMode RFAUVBinner::iterRow ( uInt irow )
158 : {
159 : uInt ant1,ant2,ifr;
160 0 : chunk.ifrToAnt(ant1,ant2,ifr=chunk.ifrNum(irow));
161 : // skip auto-correlations
162 0 : if( ant1==ant2 )
163 0 : return RFA::CONT;
164 : // compute UV distance for this row
165 0 : Float uv = sqrt(square((*puvw)(irow)(0))+square((*puvw)(irow)(1)));
166 0 : uvdist(ifr,yvalue.position()) = uv;
167 : // compute yvalues for every unflagged pixel
168 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
169 : {
170 0 : if( flag.preFlagged(ich,ifr) )
171 0 : continue;
172 : // update UV range for this channel
173 0 : if( uv < uvmin(ich) )
174 0 : uvmin = uv;
175 0 : if( uv > uvmax(ich) )
176 0 : uvmax = uv;
177 : // compute y value and update y ranges
178 0 : Float yval = mapValue(ich,irow);
179 0 : yvalue(ich,ifr) = yval;
180 0 : if( yval < ymin(ich) )
181 0 : ymin(ich) = yval;
182 0 : if( yval > ymax(ich) )
183 0 : ymax(ich) = yval;
184 : }
185 0 : return RFA::CONT;
186 : }
187 :
188 0 : RFA::IterMode RFAUVBinner::endData ()
189 : {
190 : // compute bin sizes
191 0 : uvbinsize.resize();
192 0 : uvbinsize = (uvmax-uvmin)/nbin_uv;
193 0 : ybinsize.resize();
194 0 : ybinsize = (ymax-ymin)/nbin_y;
195 :
196 0 : RFAFlagCubeBase::endData();
197 : // uInt dum;
198 : // rowclipper.updateSigma(dum,dum);
199 0 : return RFA::DRY;
200 : }
201 :
202 :
203 0 : void RFAUVBinner::startDry (bool verbose)
204 : {
205 0 : RFAFlagCubeBase::startDry(verbose);
206 : // reset lattices to read-only
207 0 : yvalue.reset(true,false);
208 : // create bincounts cube, if necessary
209 0 : if( !binned )
210 : {
211 0 : bincounts.resize();
212 0 : bincounts = Cube<Int>(nbin_uv,nbin_y,num(CHAN),0);
213 0 : totcounts.resize();
214 0 : totcounts = Vector<Int>(num(CHAN),0);
215 : }
216 0 : }
217 :
218 0 : IPosition RFAUVBinner::computeBin( Float uv,Float y,uInt ich )
219 : {
220 0 : uInt i = (uInt)((uv-uvmin(ich))/uvbinsize(ich));
221 0 : uInt j = (uInt)((y -ymin(ich))/ybinsize(ich));
222 : // loss of precision near max values can sometimes put us into bin
223 : // N+1, so correct for this:
224 0 : if( i >= nbin_uv )
225 0 : i = nbin_uv-1;
226 0 : if( j >= nbin_y )
227 0 : j = nbin_y-1;
228 0 : return IPosition(3,i,j,ich);
229 : }
230 :
231 0 : RFA::IterMode RFAUVBinner::iterDry (uInt it)
232 : {
233 0 : RFAFlagCubeBase::iterDry(it);
234 0 : yvalue.advance(it);
235 : // already binned? Do flagging
236 0 : if( binned )
237 : {
238 0 : for( uInt ifr=0; ifr<num(IFR); ifr++ )
239 : {
240 0 : Float uv = uvdist(ifr,it);
241 0 : if( uv>0 )
242 : {
243 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
244 : {
245 0 : if( !flag.preFlagged(ich,ifr) )
246 : {
247 0 : Int bc = bincounts(computeBin(uv,yvalue(ich,ifr),ich));
248 0 : if( bc<0 )
249 0 : flag.setFlag(ich,ifr);
250 : // add point to plot if in low-pop bin
251 : }
252 : }
253 : }
254 : }
255 : }
256 : // else compute bins
257 : else
258 : {
259 0 : for( uInt ifr=0; ifr<num(IFR); ifr++ )
260 : {
261 0 : Float uv = uvdist(ifr,it);
262 0 : if( uv>0 )
263 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
264 0 : if( !flag.preFlagged(ich,ifr) )
265 : {
266 0 : Float y = yvalue(ich,ifr);
267 0 : IPosition binpos( computeBin(uv,y,ich) );
268 0 : bincounts(binpos)++;
269 : // bincounts( computeBin(uv,yvalue(ich,ifr),ich) )++;
270 0 : totcounts(ich)++;
271 : }
272 : }
273 : }
274 0 : return RFA::CONT;
275 : }
276 :
277 0 : RFA::IterMode RFAUVBinner::endDry ()
278 : {
279 : // already binned? then it must have been flagged, so stop
280 0 : if( binned )
281 : {
282 0 : return RFA::STOP;
283 : }
284 : // else compute bad bins
285 0 : binned = true;
286 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
287 : {
288 : // bins for this channel
289 0 : Matrix<Int> bins( bincounts.xyPlane(ich) );
290 0 : Int maxcount = max(bins);
291 0 : Vector<Int> cumul(maxcount+1,0);
292 : // compute total population for each non-zero count
293 : // (what we compute is actually the total number of points
294 : // resident in a bin of size N, that is, N*numbins{population=N})
295 0 : for( uInt i=0; i<bins.ncolumn(); i++ )
296 0 : for( uInt j=0; j<bins.nrow(); j++ )
297 0 : if( bins(j,i) )
298 0 : cumul( bins(j,i) ) += bins(j,i);
299 : // convert to cumul(N): number of points residing in a bin of size<=N
300 : // (cumul(0) is 0 by definition)
301 0 : for( Int i=1; i<=maxcount; i++ )
302 0 : cumul(i) += cumul(i-1);
303 0 : Int thr_count=0;
304 0 : if( threshold>0 )
305 : {
306 : // compute threshold based on cumulative counts
307 0 : Float pop_cutoff = totcounts(ich)*threshold;
308 : // find the first bin count value where the cumulative bin population
309 : // is higher than the threshold
310 0 : while( thr_count<=maxcount && cumul(thr_count)<=pop_cutoff )
311 0 : thr_count++;
312 : }
313 : // if the explicit bin cut-off is higher, use it instead
314 0 : if( (Int)thr_count < (Int)min_population )
315 0 : thr_count = min_population;
316 : // thr_count is now the first population value exceeding the threshold
317 : // Bins with populations up to thr_count should be flagged
318 0 : LogicalMatrix wh( bins<thr_count );
319 0 : bins(wh) = - bins(wh);
320 : }
321 : // request another dry pass to do the flags
322 0 : return RFA::DRY;
323 : }
324 :
325 : // -----------------------------------------------------------------------
326 : // RFAUVBinner::getDesc
327 : // Returns description of parameters
328 : // -----------------------------------------------------------------------
329 0 : String RFAUVBinner::getDesc ()
330 : {
331 0 : String desc( RFDataMapper::description()+"; " );
332 : char s[256];
333 0 : if( threshold>0 )
334 : {
335 0 : sprintf(s,"%s=%g ",RF_THR,threshold);
336 0 : desc += s;
337 : }
338 0 : if( min_population )
339 : {
340 0 : sprintf(s,"%s=%d ",RF_MINPOP,min_population );
341 0 : desc += s;
342 : }
343 0 : sprintf(s,"%s=%d,%d ",RF_NBINS,nbin_uv,nbin_y);
344 0 : desc += s + RFAFlagCubeBase::getDesc();
345 0 : return desc;
346 : }
347 :
348 : // -----------------------------------------------------------------------
349 : // RFAUVBinner::getDefaults
350 : // Returns record of default parameters
351 : // -----------------------------------------------------------------------
352 0 : const RecordInterface & RFAUVBinner::getDefaults ()
353 : {
354 0 : static Record rec;
355 : // create record description on first entry
356 0 : if( !rec.nfields() )
357 : {
358 0 : rec = RFAFlagCubeBase::getDefaults();
359 0 : rec.define(RF_NAME,"UVBinner");
360 0 : rec.define(RF_COLUMN,"DATA");
361 0 : rec.define(RF_EXPR,"+ ABS XX YY");
362 0 : rec.define(RF_THR,.001);
363 0 : rec.define(RF_MINPOP,0);
364 0 : rec.define(RF_NBINS,50);
365 :
366 0 : rec.setComment(RF_COLUMN,"Use column: [DATA|MODEL|CORRected]");
367 0 : rec.setComment(RF_EXPR,"Expression for deriving value (e.g. \"ABS XX\", \"+ ABS XX YY\")");
368 0 : rec.setComment(RF_THR,"Probability cut-off");
369 0 : rec.setComment(RF_MINPOP,"Bin population cut-off");
370 0 : rec.setComment(RF_NBINS,"Number of bins (single value, or [NUV,NY])");
371 : }
372 0 : return rec;
373 : }
374 :
375 : } //# NAMESPACE CASA - END
376 :
|
