Line data Source code
1 : //# VisBuffAccumulator.cc: Implementation of VisBuffAccumulator.h
2 : //# Copyright (C) 1996,1997,1998,1999,2000,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: VisBuffAccumulator.cc,v 19.7 2006/01/17 08:22:27 gvandiep Exp $
27 : //----------------------------------------------------------------------------
28 :
29 : #include <msvis/MSVis/VisBuffAccumulator.h>
30 : #include <casacore/casa/Exceptions/Error.h>
31 : #include <casacore/casa/Arrays/ArrayLogical.h>
32 : #include <casacore/casa/Logging/LogIO.h>
33 : #include <casacore/casa/Quanta/MVTime.h>
34 : #include <iomanip>
35 :
36 : #define PRTLEV_VBA 0
37 :
38 : using namespace casacore;
39 : namespace casa { //# NAMESPACE CASA - BEGIN
40 :
41 : //----------------------------------------------------------------------------
42 :
43 250 : VisBuffAccumulator::VisBuffAccumulator(const Int& nAnt, const Double& interval,
44 250 : const Bool& prenorm, const Bool fillModel)
45 : : avBuf_p(),
46 250 : nAnt_p(nAnt),
47 250 : interval_p(interval),
48 250 : prenorm_p(prenorm),
49 : prtlev_(PRTLEV_VBA),
50 : nBuf_p(0),
51 : fillModel_p(fillModel),
52 250 : tvi_debug(False)
53 : {
54 : // Construct from the number of antennas and the averaging interval
55 : // Input:
56 : // nAnt const Int& No. of antennas
57 : // interval const Double& Time interval (in seconds).
58 : // prenorm const Bool& Pre-normalization flag
59 : // fillModel const Bool Whether or not to accumulate MODEL_DATA
60 : // Output to private data:
61 : // nAnt_p Int No. of antennas
62 : // interval_p Double Time interval (in seconds).
63 : // prenorm_p Bool Pre-normalization flag
64 : //
65 :
66 250 : if (prtlev()>2) cout << "VBA::VBA()" << endl;
67 :
68 : // Interval MUST be strictly greater than zero
69 250 : if (interval_p < DBL_EPSILON)
70 4 : interval_p=0.1; // TBD: is this reasonable?
71 :
72 : // Reset the averager
73 250 : reset();
74 :
75 250 : }
76 :
77 : //----------------------------------------------------------------------------
78 :
79 250 : VisBuffAccumulator::~VisBuffAccumulator()
80 : {
81 : // Null default destructor
82 : //
83 250 : if (prtlev()>2) cout << "VBA::~VBA()" << endl;
84 250 : }
85 :
86 : //----------------------------------------------------------------------------
87 :
88 250 : void VisBuffAccumulator::reset()
89 : {
90 : // Reset the averager
91 : // Output to private data:
92 : // tStart_p Double Start time of current accumulation
93 : // firstInterval_p Bool Is this the first interval ?
94 : // nChan_p Int No. of channels in the averaging buffer
95 : // avrow_p Int Starting row of current accumulation
96 : // avBuf_p VisBuffer Averaging buffer
97 : //
98 250 : tStart_p = 0.0;
99 250 : firstInterval_p = true;
100 250 : nCorr_p = 0;
101 250 : nChan_p = 0;
102 250 : avrow_p = 0;
103 250 : aveTime_p = 0.0;
104 250 : aveTimeWt_p = 0.0;
105 250 : globalTime_p = 0.0;
106 250 : globalTimeWt_p = 0.0;
107 250 : nBuf_p = 0;
108 :
109 250 : if (prtlev()>2) cout << " VBA::reset()" << endl;
110 :
111 250 : }
112 :
113 : //----------------------------------------------------------------------------
114 :
115 4159 : void VisBuffAccumulator::accumulate (const VisBuffer& vb)
116 : {
117 : // Accumulate a VisBuffer
118 : // Input:
119 : // vb const VisBuffer& VisBuffer to accumulate
120 : // Output to private data:
121 : // tStart_p Double Start time of current accumulation
122 : // nChan_p Int No. of channels in the avg. buffer
123 : // avrow_p Int Start row of current accumulation
124 : // avBuf_p CalVisBuffer Averaging buffer
125 : //
126 :
127 4159 : if (prtlev()>2) cout << " VBA::accumulate() " << endl;
128 :
129 : // Check if avBuf_p initialization required; if so,
130 : // assign to vb to establish a two-way connection to
131 : // the underlying VisIter.
132 4159 : if (firstInterval_p) {
133 : // Initialize the averaging buffer
134 : // THIS IS OLD INEFFICIENT WAY (PRE-CALVISBUFFER)
135 : // avBuf_p = vb;
136 :
137 : // Assign main meta info only
138 : // avBuf_p.assign(vb,true);
139 250 : avBuf_p.assign(vb,false);
140 : //avBuf_p.updateCoordInfo(); // This is (simplified) CalVisBuffer version!
141 250 : avBuf_p.copyCoordInfo(vb, true); // This only goes to VisIter if vb is
142 : // missing something needed.
143 :
144 : // Zero row count
145 250 : avBuf_p.nRow();
146 250 : avBuf_p.nRow() = 0;
147 :
148 : // Immutables:
149 250 : nChan_p = vb.nChannel();
150 250 : nCorr_p = vb.corrType().nelements();
151 :
152 : // Initialize the first accumulation interval
153 : // without copy (nothing accumulated yet)
154 250 : initialize(false);
155 :
156 : }
157 :
158 : // Iterate through the current VisBuffer
159 4159 : Int row = 0;
160 181874 : while (row < vb.nRow()) {
161 : // Find the next unflagged time
162 177715 : while (row < vb.nRow() && vb.flagRow()(row)) {
163 0 : row++;
164 : }
165 177715 : if (row < vb.nRow()) {
166 177715 : Double thisTime = vb.time()(row);
167 :
168 : // Check for the first accumulation interval
169 177715 : if (firstInterval_p) {
170 250 : tStart_p = thisTime;
171 250 : firstInterval_p = false;
172 : }
173 :
174 : // Check for end of the current accumulation interval
175 :
176 177715 : if ((vb.time()(row) - tStart_p) > (interval_p-DBL_EPSILON)) {
177 : // Normalize
178 3824 : normalize();
179 : // Advance indices to the next accumulation interval
180 3824 : tStart_p = thisTime;
181 3824 : avrow_p += hashFunction(nAnt_p-1, nAnt_p-1) + 1;
182 : // Initialize the next accumulation interval
183 : // (copy prior preavg'd intervals)
184 3824 : initialize(true);
185 : }
186 :
187 : // Add the VisBuffer row to the current accumulation
188 : //
189 : // Only accumulate VisBuffers with the same number of channels
190 177715 : Int nCorr=vb.corrType().nelements();
191 177715 : if (vb.nChannel() != nChan_p || nCorr != nCorr_p) {
192 0 : throw(AipsError("VisBuffAccumulator: data shape does not conform"));
193 : }
194 :
195 177715 : Int ant1 = vb.antenna1()(row);
196 177715 : Int ant2 = vb.antenna2()(row);
197 :
198 : // Calculate row from antenna numbers with the hash function.
199 177715 : Int outrow = avrow_p + hashFunction (ant1, ant2);
200 :
201 : // Record the row in vb that corresponds to outrow in avBuf_p.
202 177715 : outToInRow_p[outrow] = row;
203 :
204 177715 : Vector<Float> wtM(vb.weightMat().column(row));
205 : // Row weight for timestamp ave:
206 : // (used to use vb.weight(), which read WEIGHT col
207 : // directly, which is wrong because it might be calibrated!)
208 : // The following is what vb.weight() does, but here with
209 : // the weigthMat info, which has been reset from SIGMA
210 177715 : Float wt = (wtM(0)+wtM(nCorr-1))/2;
211 :
212 : // (Prenormalization removed!)
213 :
214 : // Accumulate the visCube itself
215 : // TBD: Handle weights of channel-dep accumulation
216 : // better (i.e., if some channels are sometimes
217 : // flagged, they should accumulate with different
218 : // total weights, even if we will maintain only
219 : // a channel-indep weightMat. So, either use
220 : // weightSpectrum, or at least accumulate per
221 : // channel correctly!)
222 :
223 177715 : Int goodChan(0);
224 584380 : for (Int chn=0; chn<vb.nChannel(); chn++) {
225 406665 : if (!vb.flag()(chn,row) or tvi_debug) {// jagonzal: Always copy inputVis to outputVis
226 398796 : goodChan++;
227 398796 : avBuf_p.flag()(chn,outrow) = false;
228 1319472 : for (Int cor=0;cor<nCorr;cor++) {
229 920676 : avBuf_p.visCube()(cor,chn,outrow) +=
230 1841352 : (wtM(cor)*vb.visCube()(cor,chn,row));
231 920676 : if(fillModel_p)
232 920676 : avBuf_p.modelVisCube()(cor,chn,outrow) +=
233 1841352 : (wtM(cor)*vb.modelVisCube()(cor,chn,row));
234 : }
235 : }
236 : }
237 :
238 : // Only if there is any good channels this row
239 177715 : if (goodChan > 0) {
240 177715 : avBuf_p.flagRow()(outrow) = false;
241 177715 : avBuf_p.weight()(outrow) += wt;
242 873570 : for (Int cor=0;cor<nCorr;cor++)
243 695855 : avBuf_p.weightMat()(cor,outrow) += wtM(cor);
244 :
245 : // UVW (vector average, is this right?)
246 : // gcc-3.2 needs the multiplication on a separate line; gcc-3.4 can do
247 : // it as: avBuf_p.uvw()(outrow) += vb.uvw()(row) * Double(wt);
248 177715 : RigidVector<Double,3> wtuvw = vb.uvw()(row) * Double(wt);
249 177715 : avBuf_p.uvw()(outrow) += wtuvw;
250 :
251 : // Accumulate global timestamp average for this interval:
252 : // (subtract offset from time here to avoid roundoff problems)
253 177715 : aveTime_p += (vb.time()(row)-tStart_p) * Double(wt);
254 177715 : aveTimeWt_p += Double(wt);
255 : }
256 :
257 : // Increment the row number
258 177715 : row++;
259 : } // if (row < vb.nRow())
260 : } // while (row < vb.nRow())
261 4159 : ++nBuf_p;
262 :
263 : // jagonzal: Fill rowId and flagCube (flag only provides
264 : // an OR of the flags corresponding to all corrs per chan)
265 4159 : if (nBuf_p==1 and tvi_debug)
266 : {
267 0 : for (uInt outrow_idx=0;outrow_idx<outToInRow_p.size();outrow_idx++)
268 : {
269 0 : if (outToInRow_p(outrow_idx) > -1)
270 : {
271 0 : avBuf_p.rowIds()(outrow_idx) = vb.rowIds()(outToInRow_p(outrow_idx));
272 0 : avBuf_p.flagCube().xyPlane(outrow_idx) = vb.flagCube().xyPlane(outToInRow_p(outrow_idx));
273 : }
274 : else
275 : {
276 0 : avBuf_p.rowIds()(outrow_idx) = -1;
277 0 : avBuf_p.flagCube().xyPlane(outrow_idx) = True;
278 : }
279 : }
280 : }
281 4159 : }
282 :
283 : //----------------------------------------------------------------------------
284 :
285 250 : void VisBuffAccumulator::finalizeAverage ()
286 : {
287 : // Finalize the average, and return the result
288 : // Output:
289 : // avBuf VisBuffer& Averaged buffer
290 : //
291 :
292 250 : if (prtlev()>2) cout << " VBA::finalizeAverage()" << endl;
293 :
294 : // Normalize the current (final) accumulation interval
295 250 : normalize();
296 :
297 250 : }
298 :
299 : //----------------------------------------------------------------------------
300 :
301 4074 : void VisBuffAccumulator::initialize(const Bool& copydata)
302 : {
303 : // Initialize the next accumulation interval
304 : // Output to private data:
305 : // avBuf_p VisBuffer Averaging buffer
306 : //
307 :
308 4074 : if (prtlev()>2) cout << " VBA::initialize()" << endl;
309 :
310 :
311 :
312 :
313 : ///KG notes
314 : ////Really dangerous and unmaintainable code changing private variables of vb without
315 : ////really changing the internal connnections ...so if say some internal or vb or like here
316 : //// vb.modelVisCube decides to use its private variables ...its nrow is totally wrong ...
317 : /// after the next few lines !
318 : ////function decide to use nRow_p or nChannel_p it is bound to be wrong
319 : //// that is why i am calling the damn visCube and modelVisCube above because
320 : /// coders belive that vb.nRow should give the visbuffer number of rows when accessing /// a given array
321 :
322 4074 : Int nRowAdd = hashFunction (nAnt_p-1, nAnt_p-1) + 1;
323 4074 : avBuf_p.nRow() += nRowAdd;
324 4074 : avBuf_p.nChannel() = nChan_p;
325 :
326 : // Resize and initialize the VisBuffer columns used here
327 4074 : Int nRow = avBuf_p.nRow();
328 :
329 4074 : avBuf_p.antenna1().resize(nRow, copydata);
330 4074 : avBuf_p.antenna2().resize(nRow, copydata);
331 :
332 4074 : avBuf_p.time().resize(nRow, copydata);
333 4074 : avBuf_p.uvw().resize(nRow, copydata);
334 :
335 4074 : avBuf_p.visCube().resize(nCorr_p,nChan_p, nRow,copydata);
336 4074 : if(fillModel_p)
337 4074 : copydata ?
338 3824 : avBuf_p.modelVisCube().resize(nCorr_p,nChan_p, nRow,copydata):
339 4074 : avBuf_p.setModelVisCube(Cube<Complex>(nCorr_p,nChan_p, nRow));
340 :
341 4074 : avBuf_p.weight().resize(nRow, copydata);
342 4074 : avBuf_p.weightMat().resize(nCorr_p,nRow, copydata);
343 :
344 4074 : avBuf_p.flagRow().resize(nRow, copydata);
345 :
346 4074 : avBuf_p.flag().resize(nChan_p, nRow,copydata);
347 :
348 : // jagonzal: Fill rowId and flagCube (flag only provides
349 : // an OR of the flags corresponding to all corrs per chan)
350 4074 : if (tvi_debug)
351 : {
352 0 : avBuf_p.rowIds().resize(nRow, copydata);
353 0 : avBuf_p.flagCube().resize(nCorr_p,nChan_p, nRow,copydata);
354 : }
355 :
356 : // Setup the map from avBuf_p's row numbers to input row numbers.
357 4074 : outToInRow_p.resize(nRow, copydata);
358 4074 : if(!copydata)
359 250 : outToInRow_p = -1; // Unfilled rows point to -1.
360 :
361 : // Fill in the antenna numbers for all rows
362 4074 : Int row = avrow_p;
363 43914 : for (Int ant1=0; ant1 < nAnt_p; ant1++) {
364 256380 : for (Int ant2 = ant1; ant2 < nAnt_p; ant2++) {
365 216540 : avBuf_p.antenna1()(row) = ant1;
366 216540 : avBuf_p.antenna2()(row) = ant2;
367 216540 : row++;
368 : }
369 : }
370 :
371 : // Initialize everything else
372 220614 : for (row = avrow_p; row < nRow; row++) {
373 216540 : avBuf_p.time()(row) = 0.0;
374 216540 : avBuf_p.uvw()(row) = 0.0;
375 823890 : for (Int chn = 0; chn < nChan_p; chn++) {
376 607350 : avBuf_p.flag()(chn,row) = true;
377 1860090 : for (Int cor=0; cor < nCorr_p; cor++) {
378 1252740 : avBuf_p.visCube()(cor,chn,row) = Complex(0.0);
379 1252740 : if(fillModel_p)
380 1252740 : avBuf_p.modelVisCube()(cor,chn,row) = Complex(0.0);
381 : }
382 : }
383 216540 : avBuf_p.weight()(row) = 0.0f;
384 216540 : avBuf_p.weightMat().column(row) = 0.0f;
385 216540 : avBuf_p.flagRow()(row) = true;
386 : }
387 :
388 : // Init global timestamp
389 4074 : aveTime_p = 0.0;
390 4074 : aveTimeWt_p = 0.0;
391 :
392 4074 : }
393 :
394 : //----------------------------------------------------------------------------
395 :
396 4074 : void VisBuffAccumulator::normalize()
397 : {
398 : // Normalize the current accumulation interval
399 : // Output to private data:
400 : // avBuf_p VisBuffer& Averaged buffer
401 : //
402 :
403 4074 : if (prtlev()>2) cout << " VBA::normalize()" << endl;
404 :
405 : // Only if there will be a valid timestamp
406 4074 : if (aveTimeWt_p > 0.0 ) {
407 :
408 : // Contribute to global timestamp average
409 4074 : globalTime_p-=tStart_p;
410 4074 : globalTime_p*=globalTimeWt_p;
411 4074 : globalTime_p+=aveTime_p;
412 4074 : globalTimeWt_p+=aveTimeWt_p;
413 4074 : globalTime_p/=globalTimeWt_p;
414 4074 : globalTime_p+=tStart_p;
415 :
416 : // Mean timestamp for this interval
417 4074 : aveTime_p/=aveTimeWt_p;
418 4074 : aveTime_p+=tStart_p;
419 :
420 : // Divide by the weights
421 220614 : for (Int row=avrow_p; row<avBuf_p.nRow(); row++) {
422 216540 : Float wt=avBuf_p.weight()(row);
423 433080 : Vector<Float> wtM(avBuf_p.weightMat().column(row));
424 216540 : if (wt==0.0f) avBuf_p.flagRow()(row)=true;
425 216540 : if (!avBuf_p.flagRow()(row)) {
426 175930 : avBuf_p.time()(row)=aveTime_p;
427 175930 : avBuf_p.uvw()(row)*=(1.0f/wt);
428 580810 : for (Int chn=0; chn<avBuf_p.nChannel(); chn++) {
429 404880 : if (!avBuf_p.flag()(chn,row)) {
430 1315902 : for (Int cor=0;cor<nCorr_p;cor++) {
431 918891 : if (wtM(cor)>0.0f) {
432 574371 : avBuf_p.visCube()(cor,chn,row)*=1.0f/wtM(cor);
433 574371 : if(fillModel_p)
434 574371 : avBuf_p.modelVisCube()(cor,chn,row)*=1.0f/wtM(cor);
435 : }
436 : else {
437 344520 : avBuf_p.visCube()(cor,chn,row)=0.0;
438 344520 : if(fillModel_p)
439 344520 : avBuf_p.modelVisCube()(cor,chn,row)=0.0;
440 : }
441 : }
442 : }
443 : }
444 : }
445 : }
446 : }
447 : else {
448 : // Ensure this interval is entirely flagged
449 0 : for (Int row=avrow_p; row<avBuf_p.nRow(); row++) {
450 0 : avBuf_p.flagRow()(row)=true;
451 0 : avBuf_p.weight()(row)=0.0f;
452 0 : avBuf_p.weightMat().column(row)=0.0f;
453 : }
454 : }
455 4074 : }
456 :
457 : //----------------------------------------------------------------------------
458 :
459 185613 : Int VisBuffAccumulator::hashFunction (const Int& ant1, const Int& ant2)
460 : {
461 : // Compute row index in an accumulation interval for an
462 : // interferometer index (ant1, ant2).
463 : // Input:
464 : // ant1 const Int& Antenna 1
465 : // ant2 const Int& Antenna 2
466 : // Output:
467 : // hashFunction Int Row offset in current accumulation
468 : //
469 : Int index;
470 185613 : index = nAnt_p * ant1 - (ant1 * (ant1 - 1)) / 2 + ant2 - ant1;
471 185613 : return index;
472 : }
473 :
474 : //----------------------------------------------------------------------------
475 :
476 0 : void VisBuffAccumulator::throw_err(const String& origin, const String &msg)
477 : {
478 0 : LogOrigin("VisBuffAccumulator", origin);
479 0 : throw(AipsError(msg));
480 : }
481 :
482 0 : void VisBuffAccumulator::reportData() {
483 :
484 0 : CalVisBuffer& cvb(aveCalVisBuff());
485 0 : Slice corrs(0,2,3), sl;
486 :
487 0 : Vector<Int> a1(cvb.antenna1()), a2(cvb.antenna2());
488 0 : Cube<Complex> vCC(cvb.visCube()(corrs,sl,sl));
489 0 : Cube<Complex> vCM(cvb.modelVisCube()(corrs,sl,sl));
490 0 : Matrix<Float> wtS(cvb.weightMat()(corrs,sl));
491 0 : Matrix<Bool> flg(cvb.flag()(sl,sl));
492 :
493 0 : cout << "Time=" << MVTime(timeStamp()/C::day).string(MVTime::YMD,7) << endl;
494 0 : cout.precision(8);
495 0 : for (Int irow=0;irow<cvb.nRow();++irow) {
496 0 : for (Int ich=0;ich<cvb.nChannel();++ich) {
497 0 : if (flg(ich,irow)) continue;
498 0 : cout << std::setw(2) << a1(irow) << "-" << std::setw(2) << a2(irow) << " ";
499 0 : if (cvb.nChannel()>1) cout << "ich=" << ich << " ";
500 0 : cout << "fl=" << flg(Slice(ich),Slice(irow)).nonDegenerate() << " ";
501 0 : cout << "wt=" << wtS(Slice(),Slice(irow)).nonDegenerate() << " ";
502 0 : cout << "cM=" << vCM(Slice(),Slice(),Slice(irow)).nonDegenerate() << " ";
503 0 : cout << "cC=" << vCC(Slice(),Slice(),Slice(irow)).nonDegenerate() << " ";
504 0 : cout << endl;
505 : }
506 : }
507 :
508 :
509 :
510 0 : }
511 :
512 :
513 :
514 : } //# NAMESPACE CASA - END
515 :
|
