Line data Source code
1 : //# RFASpectralRej.cc: this defines RFASpectralRej
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: 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/RFASpectralRej.h>
29 : #include <casacore/scimath/Functionals/Polynomial.h>
30 : #include <msvis/MSVis/VisibilityIterator.h>
31 : #include <msvis/MSVis/VisBuffer.h>
32 : #include <casacore/casa/Arrays/ArrayMath.h>
33 : #include <casacore/casa/Arrays/ArrayLogical.h>
34 : #include <casacore/casa/Arrays/Slice.h>
35 : #include <casacore/casa/System/PGPlotterInterface.h>
36 :
37 : using namespace casacore;
38 : namespace casa { //# NAMESPACE CASA - BEGIN
39 :
40 0 : void RFASpectralRej::addSegment ( Int spwid,Double fq0,Double fq1,Int ch0,Int ch1 )
41 : {
42 0 : Segment seg = { spwid,fq0,fq1,ch0,ch1 };
43 0 : uInt n = segments.nelements();
44 0 : segments.resize(n+1);
45 0 : segments[n] = seg;
46 0 : }
47 :
48 : // -----------------------------------------------------------------------
49 : // parseRegion
50 : // Helper function to parse one region specification. Returns
51 : // the number of ranges parsed.
52 : // -----------------------------------------------------------------------
53 0 : void RFASpectralRej::parseRegion ( const RecordInterface &parm )
54 : {
55 0 : Bool parsed = false;
56 0 : Int spwid = -1;
57 0 : if( isFieldSet(parm,RF_SPWID) && fieldType(parm,RF_SPWID,TpInt) )
58 0 : spwid = parm.asInt(RF_SPWID) - indexingBase();
59 : // figure out how channel ranges are specified - frequencies or channel #'s
60 : // First version is frequencies
61 0 : if( fieldSize(parm,RF_FREQS)>1 )
62 : {
63 0 : Array<Double> freqarr;
64 : try {
65 0 : freqarr = parm.toArrayDouble(RF_FREQS);
66 : // make sure array is of the right shape (can be reformed into 2xN)
67 0 : if( freqarr.ndim()>2 || (freqarr.nelements()%2) !=0 )
68 0 : throw( AipsError("") );
69 0 : } catch( AipsError x ) {
70 0 : os<<"Illegal \""<<RF_FREQS<<"\" array\n"<<LogIO::EXCEPTION;
71 : }
72 0 : Matrix<Double> fq( freqarr.reform(IPosition(2,2,freqarr.nelements()/2)) );
73 : // enqueue the region specs
74 0 : for( uInt i=0; i<fq.ncolumn(); i++ )
75 : {
76 0 : if( fq(0,i) >= fq(1,i) )
77 : {
78 : char s[128];
79 0 : sprintf(s,"Illegal spectral region %0.2f-%0.2f\n",fq(0,i),fq(1,i));
80 0 : os<<s<<LogIO::EXCEPTION;
81 : }
82 0 : addSegment(spwid,fq(0,i),fq(1,i),-1,-1);
83 : }
84 0 : parsed=true;
85 : }
86 : // second option is channel numbers
87 0 : if( fieldSize(parm,RF_CHANS)>1 )
88 : {
89 0 : Array<Int> arr;
90 : try {
91 0 : arr = parm.toArrayInt(RF_CHANS);
92 : // make sure array is of the right shape (can be reformed into 2xN)
93 0 : if( arr.ndim()>2 || (arr.nelements()%2) !=0 )
94 0 : throw( AipsError("") );
95 0 : } catch( AipsError x ) {
96 0 : os<<"Illegal \""<<RF_CHANS<<"\" array\n"<<LogIO::EXCEPTION;
97 : }
98 0 : Matrix<Int> ch( arr.reform(IPosition(2,2,arr.nelements()/2)) );
99 : // enqueue the region specs
100 0 : for( uInt i=0; i<ch.ncolumn(); i++ )
101 : {
102 0 : if( ch(0,i) >= ch(1,i) )
103 : {
104 : char s[128];
105 0 : sprintf(s,"Illegal spectral region #%d-%d\n",ch(0,i),ch(1,i));
106 0 : os<<s<<LogIO::EXCEPTION;
107 : }
108 0 : ch -= (Int)indexingBase();
109 0 : addSegment(spwid,0,0,ch(0,i),ch(1,i));
110 : }
111 0 : parsed=true;
112 : }
113 0 : if( !parsed )
114 0 : os<<"\""<<RF_FREQS<<"\" or \""<<RF_CHANS<<"\" must be specified\n"<<LogIO::EXCEPTION;
115 0 : }
116 :
117 0 : RFASpectralRej::RFASpectralRej ( RFChunkStats &ch,const RecordInterface &parm ) :
118 : RFAFlagCubeBase(ch,parm),
119 : RFDataMapper(parm.asArrayString(RF_EXPR),parm.asString(RF_COLUMN)),
120 0 : ndeg( parm.asInt(RF_NDEG) ),
121 0 : halfwin( parm.asInt(RF_ROW_HW) ),
122 0 : threshold( parm.asDouble(RF_ROW_THR) ),
123 0 : rowclipper(chunk,flag,threshold,halfwin)
124 : {
125 : // figure out how channel ranges are specified
126 : // if a full region record is specified, parse each element
127 : // otherwise just parse the parameter record itself
128 0 : if( isValidRecord(parm,RF_REGION) ) // full region record
129 : {
130 0 : const RecordInterface ®( parm.asRecord(RF_REGION) );
131 0 : for( uInt i=0; i<reg.nfields(); i++ )
132 : {
133 0 : if( reg.type(i) != TpRecord )
134 0 : os<<"\""<<RF_REGION<<"\" must be a record of records\n"<<LogIO::EXCEPTION;
135 0 : parseRegion(reg.asRecord(i));
136 : }
137 : }
138 : else // else assume only one region specified
139 0 : parseRegion(parm);
140 :
141 0 : if( !segments.nelements() )
142 0 : os<<"No spectral region has been specified\n"<<LogIO::EXCEPTION;
143 :
144 : // set up fitter
145 0 : Polynomial<AutoDiff<Float> > poly(ndeg);
146 0 : fitter.setFunction(poly);
147 :
148 : // set up debugging info
149 0 : if( fieldType(parm,RF_DEBUG,TpArrayInt) )
150 : {
151 0 : Vector<Int> dbg;
152 0 : parm.get(RF_DEBUG,dbg);
153 0 : if(dbg.nelements() != 2 && dbg.nelements() != 3)
154 : {
155 0 : os<<"\""<<RF_DEBUG<<"\" parameter must be [NIFR,NTIME] or [ANT1,ANT2,NTIME]"<<LogIO::EXCEPTION;
156 : }
157 : }
158 0 : }
159 :
160 : // -----------------------------------------------------------------------
161 : // newChunk
162 : // At each new chunk, figure out which channels fit into the
163 : // specified fitting regions.
164 : // -----------------------------------------------------------------------
165 0 : Bool RFASpectralRej::newChunk (Int &maxmem)
166 : {
167 : // compute correlations mask, return false if fails
168 0 : corrmask = RFDataMapper::corrMask(chunk.visIter());
169 0 : if( !corrmask )
170 : {
171 0 : os<<LogIO::WARN<<"missing selected correlations, ignoring this chunk\n"<<LogIO::POST;
172 0 : return active=false;
173 : }
174 : // figure out active channels (i.e. within specified segments)
175 0 : Int spwid = chunk.visBuf().spectralWindow();
176 0 : fitchan.resize(num(CHAN));
177 0 : fitchan.set(false);
178 0 : const Vector<Double> & fq( chunk.frequency()*1e-6 );
179 0 : for( uInt i=0; i<segments.nelements(); i++)
180 : {
181 0 : const Segment &seg ( segments[i] );
182 : // compare spectral windows, if specified
183 0 : if( seg.spwid >= 0 && seg.spwid != spwid )
184 0 : continue;
185 0 : if( seg.ch0 >= 0 ) // use channel numbers
186 : {
187 0 : if( (uInt)seg.ch0 < num(CHAN) )
188 : {
189 0 : Int ch1 = num(CHAN)-1;
190 0 : if( seg.ch1 < ch1 )
191 0 : ch1 = seg.ch1;
192 0 : fitchan(Slice(seg.ch0,ch1-seg.ch0+1)) = true;
193 : }
194 : }
195 : else // use frequencies
196 : {
197 0 : fitchan = fitchan || ( fq >= seg.fq0 && fq <= seg.fq1 );
198 : }
199 : }
200 : // count number of fitted channels
201 0 : num_fitchan = 0;
202 0 : for( uInt i=0; i<num(CHAN); i++ )
203 : {
204 0 : if( fitchan(i) )
205 : {
206 0 : xnorm = i;
207 0 : num_fitchan++;
208 : }
209 : }
210 : // return if none
211 0 : os<<num_fitchan<<" channels will be fitted in this chunk\n"<<LogIO::POST;
212 0 : if( num_fitchan<ndeg+2 )
213 : {
214 0 : os<<LogIO::WARN<<"not enough channels, ignoring chunk\n"<<LogIO::POST;
215 0 : return active=false;
216 : }
217 : // finish with init
218 0 : RFAFlagCubeBase::newChunk(maxmem-=1);
219 0 : rowclipper.init(num(IFR),num(TIME));
220 0 : return active=true;
221 : }
222 :
223 0 : void RFASpectralRej::endChunk ()
224 : {
225 0 : RFAFlagCubeBase::endChunk();
226 0 : flag.cleanup();
227 0 : rowclipper.cleanup();
228 0 : }
229 :
230 0 : void RFASpectralRej::startData (bool verbose)
231 : {
232 0 : RFAFlagCubeBase::startData(verbose);
233 0 : rowclipper.reset();
234 0 : }
235 :
236 0 : RFA::IterMode RFASpectralRej::iterTime (uInt it)
237 : {
238 0 : RFAFlagCubeBase::iterTime(it);
239 0 : RFDataMapper::setVisBuffer(chunk.visBuf());
240 0 : return RFA::CONT;
241 : }
242 :
243 0 : RFA::IterMode RFASpectralRej::iterRow ( uInt irow )
244 : {
245 : // during first pass, compute diff-median. Also keep track of the AAD.
246 0 : uInt iifr = chunk.ifrNum(irow);
247 0 : uInt it = chunk.iTime();
248 0 : Bool rowfl = chunk.npass() ? flag.rowFlagged(iifr,it)
249 0 : : flag.rowPreFlagged(iifr,it);
250 0 : if( !rowfl )
251 : {
252 0 : Vector<Float> x(num_fitchan),y(num_fitchan);
253 0 : Vector<uInt> chan(num_fitchan);
254 0 : uInt np=0;
255 : // loop over channels, collect valid (non-flagged) pixels into the x and y
256 : // vectors
257 0 : for( uInt ich=0; ich<num(CHAN); ich++ )
258 : {
259 0 : if( fitchan(ich) && !(chunk.npass() ? flag.anyFlagged(ich,iifr) : flag.preFlagged(ich,iifr)) )
260 : {
261 0 : x(np) = ich/xnorm;
262 0 : y(np) = mapValue(ich,irow);
263 0 : chan(np++) = ich;
264 : }
265 : }
266 : // check that we have enough points to constrain the fit
267 0 : if( np > ndeg+2 )
268 : {
269 : // resize x/y vectors, and do the fit
270 0 : Slice S(0,np);
271 0 : Vector<Float> sigma(np,1.0f),x1(x(S)),y1(y(S));
272 0 : Vector<Float> c = fitter.fit(x1,y1,sigma);
273 0 : Float chisq = fitter.chiSquare();
274 0 : rowclipper.setSigma(iifr,it,chisq);
275 : } // endif( np>ndeg+2)
276 : } // endif( !rowfl )
277 :
278 0 : return RFA::CONT;
279 0 : }
280 :
281 : // -----------------------------------------------------------------------
282 : // endData
283 : // Ends data pass
284 : // -----------------------------------------------------------------------
285 0 : RFA::IterMode RFASpectralRej::endData ()
286 : {
287 0 : RFAFlagCubeBase::endData();
288 : uInt dum;
289 0 : rowclipper.updateSigma(dum,dum);
290 0 : return RFA::STOP;
291 : }
292 :
293 : // -----------------------------------------------------------------------
294 : // RFASpectralRej::getDesc
295 : // Returns description of parameters
296 : // -----------------------------------------------------------------------
297 0 : String RFASpectralRej::getDesc ()
298 : {
299 0 : String desc( RFDataMapper::description()+";" );
300 : char s[256];
301 : // build up description of spectral segments
302 0 : for( uInt i=0; i<segments.nelements(); i++)
303 : {
304 0 : const Segment &seg ( segments[i] );
305 : // is spwid specified?
306 : char s1[32];
307 0 : if( seg.spwid >= 0 )
308 0 : sprintf(s1,"%d:",seg.spwid);
309 : else
310 0 : s1[0]=0;
311 0 : if( seg.ch0 >= 0 ) // use channel numbers
312 0 : sprintf(s, " %s#%d-%d",s1,seg.ch0,seg.ch1);
313 : else
314 0 : sprintf(s, " %s%.2f-%.2fMHz",s1,seg.fq0,seg.fq1);
315 0 : desc+=s;
316 : }
317 0 : desc += "; ";
318 0 : sprintf(s,
319 : "%s=%d %s=%.1f %s=%d",
320 : RF_NDEG,
321 : ndeg,
322 : RF_ROW_THR,
323 : threshold,
324 : RF_ROW_HW,
325 : halfwin);
326 0 : desc += s;
327 0 : desc += RFAFlagCubeBase::getDesc();
328 0 : return desc;
329 : }
330 :
331 : // -----------------------------------------------------------------------
332 : // RFASpectralRej::getDefaults
333 : // Returns record of default parameters
334 : // -----------------------------------------------------------------------
335 0 : const RecordInterface & RFASpectralRej::getDefaults ()
336 : {
337 0 : static Record rec;
338 : // create record description on first entry
339 0 : if( !rec.nfields() )
340 : {
341 0 : rec = RFAFlagCubeBase::getDefaults();
342 0 : rec.define(RF_NAME,"SpectralRej");
343 0 : rec.define(RF_COLUMN,"DATA");
344 0 : rec.define(RF_EXPR,"+ ABS XX YY");
345 0 : rec.define(RF_NDEG,(Int)2);
346 0 : rec.define(RF_ROW_THR,(Double)5);
347 0 : rec.define(RF_ROW_HW,(Int)6);
348 0 : rec.define(RF_REGION,false);
349 0 : rec.define(RF_SPWID,false);
350 0 : rec.define(RF_FREQS,false);
351 0 : rec.define(RF_CHANS,false);
352 0 : rec.define(RF_DEBUG,false);
353 :
354 0 : rec.setComment(RF_COLUMN,"Use column: [DATA|MODEL|CORRected]");
355 0 : rec.setComment(RF_EXPR,"Expression for deriving value (e.g. \"ABS XX\", \"+ ABS XX YY\")");
356 0 : rec.setComment(RF_NDEG,"Number of degrees for polynomial fit");
357 0 : rec.setComment(RF_ROW_THR,"Row flagging threshold, in AADs");
358 0 : rec.setComment(RF_ROW_HW,"Row flagging, half-window of sliding median (0 to use overall median)");
359 0 : rec.setComment(RF_SPWID,"Spectral window ID (F or -1 for all)");
360 0 : rec.setComment(RF_FREQS,"Range(s) of frequencies (in MHz)");
361 0 : rec.setComment(RF_CHANS,"Range(s) of channel numbers");
362 0 : rec.setComment(RF_REGION,"For several spectral regions, record of records");
363 0 : rec.setComment(RF_DEBUG,"Set to [NIFR,NTIME] or [ANT1,ANT2,NTIME] to produce debugging plots");
364 : }
365 0 : return rec;
366 : }
367 :
368 :
369 :
370 :
371 : } //# NAMESPACE CASA - END
372 :
|
