Line data Source code
1 : //# RFChunkStats.cc: this defines RFChunkStats
2 : //# Copyright (C) 2000,2001,2002,2003
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 <casacore/casa/Arrays/ArrayMath.h>
28 : #include <casacore/casa/Exceptions/Error.h>
29 : #include <flagging/Flagging/Flagger.h>
30 : #include <flagging/Flagging/RFChunkStats.h>
31 : #include <msvis/MSVis/VisibilityIterator.h>
32 : #include <msvis/MSVis/VisBuffer.h>
33 :
34 : using namespace casacore;
35 : namespace casa { //# NAMESPACE CASA - BEGIN
36 :
37 0 : RFChunkStats::RFChunkStats( VisibilityIterator &vi,VisBuffer &vb,Flagger &rf )
38 : : visiter(vi),
39 : visbuf(vb),
40 0 : flagger(rf)
41 : {
42 0 : chunk_no=0;
43 0 : counts[ANT] = flagger.numAnt();
44 0 : counts[IFR] = flagger.numIfr();
45 0 : counts[FEED] = flagger.numFeed();
46 0 : counts[FEEDCORR] = flagger.numFeedCorr();
47 0 : }
48 :
49 0 : const MeasurementSet & RFChunkStats::measSet () const
50 : {
51 0 : return flagger.measSet();
52 : }
53 :
54 0 : const String RFChunkStats::msName () const
55 : {
56 0 : return flagger.measSet().tableName();
57 : }
58 :
59 0 : const Vector<String> & RFChunkStats::antNames () const
60 : {
61 0 : return flagger.antNames();
62 : }
63 :
64 0 : const Vector<Double> & RFChunkStats::frequency ()
65 : {
66 0 : return visiter.frequency(freq);
67 : }
68 :
69 0 : void RFChunkStats::newChunk(bool init_quack)
70 : {
71 0 : itime=-1;
72 0 : chunk_no++;
73 : // setup shapes
74 0 : visshape = visiter.visibilityShape(); //4
75 0 : counts[POLZN] = visshape(0);
76 0 : counts[CHAN] = visshape(1);
77 :
78 0 : counts[TIME] = visiter.nSubInterval();
79 0 : counts[ROW] = visiter.nRowChunk();
80 :
81 : // get correlation types
82 0 : visiter.corrType(corrtypes);
83 : // reset min/max time slots
84 0 : start_time = end_time = current_time = 0;
85 :
86 : // setups statistics
87 :
88 0 : nrf_time.resize(num(TIME));
89 0 : nrf_time.set(0);
90 0 : nrf_ifr.resize(num(IFR));
91 0 : nrf_ifr.set(0);
92 0 : rows_per_ifr.resize(num(IFR));
93 0 : rows_per_ifr.set(0);
94 :
95 0 : nf_ifr_time.resize(num(IFR),num(TIME));
96 0 : nf_ifr_time.set(0);
97 0 : nf_chan_ifr.resize(num(CHAN),num(IFR));
98 0 : nf_chan_ifr.set(0);
99 : // nf_corr_ifr.resize(num(CORR),num(IFR));
100 : // nf_corr_ifr.set(0);
101 : // nf_chan_corr.resize(num(CHAN),num(CORR));
102 : // nf_chan_corr.set(0);
103 : // nf_chan_time.resize(num(CHAN),num(TIME));
104 : // nf_chan_time.set(0);
105 : // nf_corr_time.resize(num(CORR),num(TIME));
106 : // nf_corr_time.set(0);
107 :
108 : if (0) {
109 : cout << "Ifr : " << num(IFR) << " Chan : " << num(CHAN) << " Corr : " << num(CORR) << " Time : " << num(TIME) << endl;
110 : }
111 :
112 : // build up description of correlations
113 0 : corr_string = "";
114 0 : for( uInt i=0; i<corrtypes.nelements(); i++ )
115 0 : corr_string += " " + Stokes::name( Stokes::type(corrtypes(i)) );
116 : char s[256];
117 :
118 0 : sprintf(s,"Chunk %d : [field: %d, spw: %d] %s, %d channels, %d time slots, %d baselines, %d rows\n", chunk_no,visBuf().fieldId(),visBuf().spectralWindow(),corr_string.chars(),num(CHAN),num(TIME),num(IFR),num(ROW));
119 : // Flagger::logSink()<<s<<LogIO::POST;
120 :
121 0 : if (init_quack) {
122 : // figure out this scan's start and end times
123 : // for use in quack mode
124 0 : for (visiter.origin();
125 0 : visiter.more();
126 0 : visiter++) {
127 :
128 0 : int scan = visbuf.scan()(0);
129 :
130 0 : if (scan < 0) {
131 :
132 0 : std::string s;
133 0 : std::stringstream ss;
134 0 : ss << scan;
135 0 : s = ss.str();
136 0 : throw AipsError("Scan number must be non-negative, is " + s);
137 : }
138 0 : const Vector<Double> ×(visbuf.time());
139 :
140 0 : double t0 = times(0);
141 :
142 : /* Figure out if anything is flagged in this buffer */
143 0 : Bool any_unflagged = false;
144 : casacore::Bool dataIsAcopy;
145 0 : casacore::Bool * flags = visbuf.flagCube().getStorage(dataIsAcopy);
146 :
147 0 : unsigned n = visbuf.flagCube().nelements();
148 0 : for (unsigned i = 0; i < n; i++) {
149 : //cout << "flags[" << i << " = " << flags[i] << endl;
150 0 : if (!flags[i]) {
151 0 : any_unflagged = true;
152 0 : break;
153 : }
154 : }
155 0 : visbuf.flagCube().putStorage(flags, dataIsAcopy);
156 :
157 : //cout << "flagCube() = " << visbuf.flagCube() << endl;
158 : //cout << "flag() = " << visbuf.flag() << endl;
159 : //cout << "flagRow() = " << visbuf.flagRow() << endl;
160 :
161 0 : if (scan_start.size() < (unsigned) scan+1) {
162 : // initialize to -1
163 0 : scan_start .resize(scan+1, -1.0);
164 0 : scan_start_flag.resize(scan+1, -1.0);
165 0 : scan_end .resize(scan+1, -1.0);
166 0 : scan_end_flag.resize(scan+1, -1.0);
167 : }
168 0 : if (scan_start[scan] < 0 ||
169 0 : t0 < scan_start[scan]) {
170 0 : scan_start[scan] = t0;
171 : }
172 0 : if (scan_end[scan] < 0 ||
173 0 : t0 > scan_end[scan-1]) {
174 0 : scan_end[scan] = t0;
175 : }
176 :
177 : if (0) cout << "any_unflags = " << any_unflagged << endl;
178 0 : if (any_unflagged) {
179 0 : if (scan_start_flag[scan] < 0 ||
180 0 : t0 < scan_start_flag[scan]) {
181 0 : scan_start_flag[scan] = t0;
182 : }
183 0 : if (scan_end_flag[scan] < 0 ||
184 0 : t0 > scan_end_flag[scan-1]) {
185 0 : scan_end_flag[scan] = t0;
186 : }
187 : }
188 : }
189 : if (0) for (unsigned int i = 0; i < scan_start.size(); i++) {
190 :
191 : cerr << "scan " << i << " start = " <<
192 : MVTime( scan_start[i]/C::day).string(MVTime::DMY,7)
193 : << " end = " <<
194 : MVTime( scan_end[i]/C::day).string(MVTime::DMY,7) << endl;
195 :
196 : cerr << "scan " << i << " start flag = " <<
197 : MVTime( scan_start_flag[i]/C::day).string(MVTime::DMY,7)
198 : << " end = " <<
199 : MVTime( scan_end_flag[i]/C::day).string(MVTime::DMY,7) << endl;
200 :
201 : }
202 : }
203 0 : }
204 :
205 0 : void RFChunkStats::newPass (uInt npass)
206 : {
207 0 : itime = -1;
208 0 : pass_no = npass;
209 0 : }
210 :
211 0 : void RFChunkStats::newTime ()
212 : {
213 : // setup IFR numbers for every row in time slot
214 0 : ifr_nums.resize( visbuf.antenna1().nelements() );
215 0 : ifr_nums = flagger.ifrNumbers( visbuf.antenna1(),
216 0 : visbuf.antenna2() );
217 :
218 : // setup FEED CORRELATION numbers for every row in time slot
219 0 : feed_nums.resize( visbuf.feed1().nelements() );
220 0 : feed_nums = flagger.ifrNumbers( visbuf.feed1(),visbuf.feed2() );
221 :
222 : // reset stats
223 0 : for( uInt i=0; i<ifr_nums.nelements(); i++ ) {
224 :
225 0 : Int baseline_num = ifr_nums(i);
226 :
227 0 : if (baseline_num >= (Int) rows_per_ifr.nelements()) {
228 0 : std::stringstream ss;
229 0 : ss << "Internal error: Baseline number is " << baseline_num
230 0 : << ", but there are only " << rows_per_ifr.nelements() << " baselines" << endl;
231 0 : throw AipsError(ss.str());
232 : }
233 :
234 0 : rows_per_ifr(ifr_nums(i))++;
235 : }
236 :
237 : // set start/end times
238 0 : current_time = (visbuf.time()(0))/(24*3600);
239 0 : if( current_time<start_time || start_time==0 )
240 0 : start_time = current_time;
241 0 : if( current_time>end_time )
242 0 : end_time = current_time;
243 :
244 0 : itime++;
245 : // dprintf(os,"newTime: %d\n",itime);
246 0 : }
247 :
248 0 : uInt RFChunkStats::antToIfr ( uInt ant1,uInt ant2 )
249 : {
250 0 : return flagger.ifrNumber((Int)ant1,(Int)ant2);
251 : }
252 :
253 0 : void RFChunkStats::ifrToAnt ( uInt &ant1,uInt &ant2,uInt ifr )
254 : {
255 0 : flagger.ifrToAnt(ant1,ant2,ifr);
256 0 : }
257 :
258 0 : String RFChunkStats::ifrString ( uInt ifr )
259 : {
260 : uInt a1,a2;
261 0 : ifrToAnt(a1,a2,ifr);
262 : char s[32];
263 0 : sprintf(s,"%d-%d",a1,a2);
264 0 : return String(s);
265 : }
266 :
267 0 : void RFChunkStats::printStats ()
268 : {
269 : // collapse 2D stats into 1D arrays
270 0 : Vector<uInt> nf_ifr(num(IFR)),nf_chan(num(CHAN)),nf_time(num(TIME));
271 0 : for( uInt i=0; i<num(TIME); i++)
272 0 : nf_time(i) = sum(nf_ifr_time.column(i));
273 0 : for( uInt i=0; i<num(IFR); i++)
274 0 : nf_ifr(i) = sum(nf_ifr_time.row(i));
275 0 : for( uInt i=0; i<num(CHAN); i++)
276 0 : nf_chan(i) = sum(nf_chan_ifr.row(i));
277 :
278 0 : uInt n = sum(nf_ifr);
279 0 : fprintf( stderr,"%d pixel flags have been raised\n",n );
280 0 : fprintf( stderr,"(%f%% of pixels flagged)\n",n*100.0/(num(ROW)*num(CHAN)*num(CORR)));
281 0 : n = sum(nrf_ifr);
282 0 : fprintf( stderr,"%d row flags have been raised\n",n );
283 0 : fprintf( stderr,"(%f%% of rows flagged)\n",n*100.0/num(ROW));
284 : // accumulate per-antenna flags
285 0 : Vector<uInt> flperant(num(ANT),0u),rowperant(num(ANT),0u);
286 0 : for( uInt i=0; i<num(IFR); i++ )
287 : {
288 : uInt a1,a2;
289 0 : ifrToAnt(a1,a2,i);
290 0 : if( nf_ifr(i) )
291 : {
292 0 : flperant(a1) += nf_ifr(i);
293 0 : flperant(a2) += nf_ifr(i);
294 : }
295 0 : if( nrf_ifr(i) )
296 : {
297 0 : rowperant(a1) += nrf_ifr(i);
298 0 : rowperant(a2) += nrf_ifr(i);
299 : }
300 : }
301 :
302 0 : cerr<<"Flags per antenna: "<<flperant<<endl;
303 0 : cerr<<"Flags per IFR %: "<<(nf_ifr*100u)/(num(CHAN)*num(TIME))<<endl;
304 0 : cerr<<"Flags per chan %: "<<(nf_chan*100u)/(num(IFR)*num(TIME))<<endl;
305 0 : cerr<<"Flags per time %: "<<(nf_time*100u)/(num(CHAN)*num(IFR))<<endl;
306 0 : cerr<<"Row flags per antenna: "<<rowperant<<endl;
307 0 : cerr<<"Row flags per IFR %: "<<(nrf_ifr*100u)/(num(TIME))<<endl;
308 0 : cerr<<"Row flags per time %: "<<(nrf_time*100u)/(num(IFR))<<endl;
309 0 : }
310 :
311 : // -----------------------------------------------------------------------
312 : // findCorrType
313 : // -----------------------------------------------------------------------
314 0 : Int findCorrType( Stokes::StokesTypes type,const Vector<Int> &corr )
315 : {
316 0 : if( type == Stokes::Undefined )
317 0 : return -1;
318 : // find this type in current polarizations
319 0 : for( uInt j=0; j<corr.nelements(); j++ )
320 0 : if( type == corr(j) )
321 0 : return j;
322 0 : return -1;
323 : }
324 :
325 :
326 : } //# NAMESPACE CASA - END
327 :
|
