Line data Source code
1 : //# MsAverager.cc: Implementation of MsAverager.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: casa-feedback@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 :
29 :
30 : #include <msvis/MSVis/SelectAverageSpw.h>
31 : #include <msvis/MSVis/MsAverager.h>
32 : #include <casacore/casa/Exceptions/Error.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/Slice.h>
35 : #include <casacore/casa/OS/Timer.h>
36 : #include <iomanip>
37 :
38 : #include <stdio.h>
39 : #include <stdlib.h>
40 : #include <errno.h>
41 : #include <iostream>
42 : #include <iomanip>
43 : #include <cmath>
44 :
45 : #define LOG2 0
46 :
47 : #if LOG2
48 : #define IfLog2(x) x
49 : #else
50 : #define IfLog2(x) /*x*/
51 : #endif
52 :
53 : using namespace casacore;
54 : namespace casa {
55 :
56 : const String MsAverager::clname = "MsAverager";
57 :
58 : const String MsAverager::DataColumn[] = {"DATA", "CORRECTEDDATA",
59 : "MODELDATA", "RESIDUAL"};
60 :
61 0 : MsAverager::MsAverager(MS* ms, OutputMode mode) {
62 0 : aveOK = false;
63 :
64 : //log = SLog::slog();
65 0 : reset(ms, mode);
66 0 : }
67 :
68 0 : void MsAverager::reset(MS* ms, OutputMode mode) {
69 : //cout << "MS=" << *ms << endl;
70 0 : aveOK = false;
71 0 : outputMode = mode;
72 :
73 0 : if (!ms) {
74 0 : LogIO os(LogOrigin("MsAverager", "reset"));
75 0 : os << LogIO::WARN << "Could not reset MsAverager: input MS is NULL" << LogIO::POST;
76 : //SLog::slog()->out("Could not reset MsAverager: input MS is NULL",
77 : //"MsAverager", clname, LogMessage::WARN, true);
78 0 : return;
79 0 : }
80 :
81 : //SLog::slog()->out(String("Number of selected rows is ") +
82 : // String::toString(ms->nrow()),
83 : // "MsAverager", clname, LogMessage::DEBUG1);
84 :
85 0 : pMS = ms;
86 0 : vs = 0;
87 0 : pAveBuff = 0;
88 :
89 0 : spw.resize(0);
90 :
91 : //Block<int> sort(4);
92 : //sort[0] = MS::DATA_DESC_ID;
93 : //sort[1] = MS::FIELD_ID;
94 : //sort[2] = MS::TIME;
95 : //sort[3] = MS::ARRAY_ID;
96 :
97 : //Matrix<Int> allChannels;
98 : //Double intrvl = 2;
99 : //vs = new VisSet(*pMS, sort, allChannels, intrvl);
100 :
101 0 : TableDesc td = pMS->tableDesc();
102 0 : String aveTable=String("casapy.scratch-0.table");
103 0 : if (!Table::isOpened(aveTable)) {
104 0 : SetupNewTable aveSetup(aveTable, td, Table::Scratch);
105 0 : StManAipsIO stm;
106 0 : aveSetup.bindAll(stm);
107 :
108 : //if (mspointer!=NULL) delete mspointer;
109 : //mspointer=0;
110 : //MS* mspointer= SubMS::setupMS(aveTable, 10, 4, "VLA");
111 : //aMS=*mspointer;
112 :
113 0 : aMS = MS();
114 0 : aMS = MS(aveSetup, 0, false);
115 0 : msc = new MSColumns(aMS);
116 : //cout << "---aMS=" << aMS << endl;
117 0 : }
118 :
119 0 : msName = pMS->tableName();
120 : //msdv = new MSDerivedValues;
121 : //msdv->setMeasurementSet(*pMS);
122 : //msdv->setVelocityReference(MDoppler::RADIO);
123 : //msdv->setFrequencyReference(MFrequency::LSRK);
124 0 : }
125 :
126 0 : MsAverager::~MsAverager() {
127 0 : cleanup();
128 0 : }
129 :
130 0 : void MsAverager::cleanup() {
131 0 : if (pMS) {
132 : //delete pMS;
133 0 : pMS = 0;
134 : }
135 0 : if (vs) {
136 : //delete vs;
137 0 : vs = 0;
138 : }
139 0 : if (msc) {
140 : //delete msc;
141 0 : msc = 0;
142 : }
143 0 : if (pAveBuff) {
144 : //delete pAveBuff;
145 0 : pAveBuff = 0;
146 : }
147 0 : }
148 :
149 :
150 0 : void MsAverager::setAverager(
151 : const Matrix<Int>& cList,
152 : const Matrix<Int>& bls,
153 : Double aveT, Int aveC,
154 : const String& col,
155 : const String& aveM,
156 : const Bool& aveF,
157 : const Bool& crossscan,
158 : const Bool& crossbline,
159 : const Bool& crossarray,
160 : const Bool& aveVelo,
161 : const String& restfre,
162 : const String& fram,
163 : const String& doppl) {
164 0 : String fnname = "average";
165 : //static uInt loop = 0;
166 :
167 0 : for (Int i = aMS.nrow() - 1; i >= 0; i--) {
168 0 : aMS.removeRow(i);
169 : }
170 : //cout << "aMS rows=" << aMS.nrow() << endl;
171 : //showColumnNames();
172 0 : LogIO os(LogOrigin("MsAverager", "setAverager"));
173 0 : if (upcase(col)=="MODEL")
174 0 : if (!hasColumn("MODEL_DATA")) {
175 0 : os << LogIO::WARN << "The MS does not have MODEL_DATA column" << LogIO::POST;
176 0 : return;
177 : }
178 0 : if (upcase(col)=="CORRECTED")
179 0 : if (!hasColumn("CORRECTED_DATA")) {
180 0 : os << LogIO::WARN << "The MS does not have CORRECTED_DATA column" << LogIO::POST;
181 0 : return;
182 : }
183 0 : if (upcase(col)=="RESIDUAL") {
184 0 : if (!hasColumn("CORRECTED_DATA")) {
185 0 : os << LogIO::WARN << "The MS does not have CORRECTED_DATA column" << LogIO::POST;
186 0 : return;
187 : }
188 0 : if (!hasColumn("MODEL_DATA")) {
189 0 : os << LogIO::WARN << "The MS does not have MODEL_DATA column" << LogIO::POST;
190 0 : return;
191 : }
192 : }
193 :
194 0 : column = col;
195 0 : if (column == "CORRECTED") {
196 0 : column = "CORRECTEDDATA";
197 : }
198 0 : if (column == "MODEL") {
199 0 : column = "MODELDATA";
200 : }
201 0 : if (column == "RESIDUAL") {
202 0 : column = "RESIDUAL";
203 : }
204 :
205 0 : crossSpws = false;
206 0 : if (aveC==1234567) {
207 0 : crossSpws = true;
208 : }
209 :
210 0 : Matrix<Int> allChannels;
211 0 : Double intrvl = 2;
212 0 : if (crossSpws) {
213 : //cout << "crossSpws=true sort by time" << endl;
214 0 : Block<int> sort(2);
215 0 : sort[0] = MS::TIME;
216 : //sort[1] = MS::FIELD_ID;
217 0 : sort[1] = MS::DATA_DESC_ID;
218 : //sort[3] = MS::ARRAY_ID;
219 0 : vs = new VisSet(*pMS, sort, allChannels, intrvl);
220 0 : }
221 : else {
222 : //cout << "crossSpws=false sort by spw" << endl;
223 0 : Block<int> sort(4);
224 0 : sort[0] = MS::DATA_DESC_ID;
225 0 : sort[1] = MS::FIELD_ID;
226 0 : sort[2] = MS::TIME;
227 0 : sort[3] = MS::ARRAY_ID;
228 0 : vs = new VisSet(*pMS, sort, allChannels, intrvl);
229 0 : }
230 :
231 0 : aveFlag = aveF;
232 0 : aveMode = aveM;
233 0 : aveChan = aveC;
234 0 : aveTime = aveT;
235 0 : crossScans = crossscan;
236 0 : crossBlines = crossbline;
237 0 : crossArrays = crossarray;
238 0 : aveVel = aveVelo;
239 0 : sorryVel = true;
240 :
241 0 : baselines.resize(bls.shape()[0], bls.shape()[1]);
242 0 : baselines = bls;
243 :
244 0 : restfreq = restfre;
245 0 : frame = fram;
246 0 : doppler = doppl;
247 :
248 : //cout
249 : // << "aveFlag=" << aveFlag << " aveMode=" << aveMode
250 : // << " aveChan=" << aveChan << " aveTime=" << aveTime
251 : // << " aveC=" << aveC << " crossScans=" << crossscan
252 : // << " crossBlines=" << crossbline
253 : // << " crossArrays=" << crossarray
254 : // << " crossSpws=" << crossSpws
255 : // << " aveVel=" << aveVel
256 : // << " baselines=" << baselines
257 : // << " restfreq=" << restfreq
258 : // << " frame=" << frame
259 : // << " doppler=" << doppler
260 : // << endl;
261 :
262 0 : nAvePol = -1;
263 0 : if (aveChan < 1)
264 0 : aveChan = 1;
265 0 : if (aveTime < 0)
266 0 : aveTime = 0;
267 :
268 0 : msRow = 0;
269 :
270 : //showColumnNames();
271 : //if (!hasColumn(upcase(col))) {
272 : // cout << "The MS does not have '" + col + "' column" << endl;
273 : // return;
274 : //}
275 :
276 : //if (aveMode == "SCALAR" &&
277 : // !((isDataColumn(xcol) && xval == "AMP") ||
278 : // (isDataColumn(ycol) && yval == "AMP"))) {
279 : // SLog::slog()->out("Scalar everage is for amplitude of visibility only",
280 : // fnname, clname, LogMessage::WARN);
281 : // return;
282 : //}
283 :
284 0 : if (aveTime == 0 && aveChan == 1) {
285 : //SLog::slog()->out("No averaging",
286 : //fnname, clname, LogMessage::NORMAL5);
287 0 : aveOK = false;
288 0 : return;
289 : }
290 :
291 : try {
292 0 : VisIter& vi(vs->iter());
293 0 : VisBuffer vb(vi);
294 :
295 : //aveRow = 0;
296 :
297 0 : Int iChunk = 0; //count the total number of chunks
298 0 : Int iIter = 0; //count the total number of iterations
299 0 : Int iRow = 0; //count the total number of input rows
300 :
301 0 : Int corrs = -1;
302 0 : Int chans = -1;
303 0 : Int corrChange = 0;
304 0 : Int chanChange = 0;
305 : //cout << "before checking the shape" << endl;
306 0 : vi.origin();
307 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
308 0 : for (vi.origin(); vi.more(); vi++) {
309 : //cout << "corrs=" << vi.nCorr()
310 : // << " chans=" << vb.nChannel() << endl;
311 0 : if (corrs != vi.nCorr()) {
312 0 : corrs = vi.nCorr();
313 0 : corrChange++;
314 : }
315 0 : if (chans != vb.nChannel()) {
316 0 : chans = vb.nChannel();
317 0 : chanChange++;
318 : //cout << "chanChange=" << chanChange << endl;
319 : }
320 : }
321 : }
322 : //cout << "after checking the shape" << endl;
323 0 : if (chanChange > 1 || corrChange > 1) {
324 0 : LogIO os(LogOrigin("MsAverager", "setAverager"));
325 : os << LogIO::WARN << "Average over variable shape of "
326 : << "channel/polarization is not supported"
327 0 : << LogIO::POST;
328 0 : aveOK = false;
329 0 : return;
330 0 : }
331 :
332 0 : Int nChan = vb.nChannel();
333 0 : aveChan = max(1, aveChan);
334 0 : aveChan = min(aveChan, nChan);
335 :
336 0 : nAveChan =
337 0 : SAS::selectAverageChan(pMS, cList, spw, aveChan);
338 : //SAS::showSASC(spw);
339 : //cout << "nAveChan=" << nAveChan << endl;
340 :
341 0 : SAS::chanMap(aveChanMap, spw);
342 : //cout << "aveChanMap=" << aveChanMap ;
343 :
344 0 : aveRowMap.resize(pMS->nrow(), 3);
345 : //pMS->antenna().nrow() = vs->numberAnt()
346 : //due to the baseline/antenna selection, not all the
347 : //antenna present in the pMS
348 : //Int nAnt = (pMS->antenna().nrow());//vs->numberAnt();
349 : //cout << "nAnt=" << nAnt << endl;
350 :
351 0 : Int nAnt = baselines.shape()[0];
352 0 : nAntenna = nAnt;
353 : //cout << "nAnt=" << nAnt << endl;
354 :
355 0 : aveTime = max(0.01, aveTime);
356 : //cout << "aveTime=" << aveTime << endl;
357 0 : vi.setInterval(aveTime);
358 :
359 0 : nAveTime = 0; //number of bins of time averaged buffer
360 0 : Int nTime = 0; //number of integrations in the current cumulating bin
361 0 : nAveRow = 0;
362 :
363 0 : Vector<Double> curTime;
364 : //Vector<Int> ant1;
365 : //Vector<Int> ant2;
366 :
367 0 : Double bufTime = -1;
368 0 : Bool newTime = true;
369 :
370 0 : Double sumTime = 0.;
371 0 : Double nextTime = -1;
372 0 : Int nTotalTime = 0;
373 :
374 0 : Int bufField = -1;
375 0 : Bool newField = true;
376 :
377 0 : Int bufScan = -1;
378 0 : Bool newScan = true;
379 :
380 0 : Int bufArray = -1;
381 0 : Bool newArray = true;
382 :
383 0 : Int bufDesc = -1;
384 0 : Bool newDesc = true;
385 :
386 : //timer for profile - cumulate time for patAveTable
387 : //double usr = 0.;
388 : //double rea = 0.;
389 : //double sys = 0.;
390 : //Timer tmr2;
391 : //
392 :
393 : //timer for the total averaging time
394 0 : Timer tmr3;
395 0 : vi.origin();
396 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
397 : #if LOG2
398 : cout << ">>>>>>>>>chunk=" << iChunk << endl;
399 : #endif
400 :
401 0 : for (vi.origin(); vi.more(); vi++) {
402 0 : if (nAvePol == -1)
403 0 : nAvePol = vi.nCorr();
404 :
405 0 : if (nAvePol != vi.nCorr())
406 0 : cout << "nAvePol=" << nAvePol << endl;
407 :
408 0 : vi.time(curTime);
409 : //vi.antenna1(ant1);
410 : //vi.antenna2(ant2);
411 :
412 0 : Int nRow = curTime.nelements();
413 : //Int nRow = vi.nRow();
414 : //Int nRow = vb.nRow();
415 :
416 : #if LOG2
417 : cout << ">>>>>>iIter=" << iIter << " nRow=" << nRow << endl;
418 : cout << std::setprecision(12);
419 : #endif
420 :
421 0 : Cube<Complex>& vc = vb.visCube();
422 0 : if (!upcase(column).compare("CORRECTEDDATA") ||
423 0 : !upcase(column).compare("CORRECTED_DATA")) {
424 0 : vc = vb.correctedVisCube();
425 : }
426 0 : if (!upcase(column).compare("MODELDATA") ||
427 0 : !upcase(column).compare("MODEL_DATA")) {
428 0 : vc = vb.modelVisCube();
429 : }
430 0 : if (!upcase(column).compare("RESIDUAL") ||
431 0 : !upcase(column).compare("CORRECTEDMODEL_DATA")) {
432 0 : vc = vb.correctedVisCube() - vb.modelVisCube();
433 : }
434 : Int i1, i2, i3;
435 0 : vc.shape(i1, i2, i3);
436 : //cout << "i1=" << i1 << " i2=" << i2 << " i3=" << i3 << endl;
437 0 : Cube<Complex> aveV(i1, nAveChan, 1);
438 0 : Cube<Bool> aveF(i1, nAveChan, 1);
439 :
440 0 : Int field = vb.fieldId();
441 0 : Int arry = vb.arrayId();
442 :
443 0 : Vector<rownr_t> rowNumber = vb.rowIds();
444 : //cout << "rowIds=" << rowNumber << endl;
445 0 : for (int row = 0; row < nRow; row++) {
446 : //skip flagged rows when selecting good data
447 : //this is done to make sure for good data average
448 : //each timebin start with good data
449 0 : while(row < nRow && vb.flagRow()(row) && aveFlag == 0) {
450 0 : row++;
451 : };
452 0 : if (row == nRow)
453 0 : break;
454 :
455 : //skip rows that are not in the selected spw
456 0 : Int dId = vi.dataDescriptionId();
457 0 : Int spwSelected = SAS::spwByDesc(dId, spw);
458 0 : while (row < nRow && spwSelected == -1) {
459 0 : row++;
460 : }
461 : //cout << "row=" << row
462 : // << " spwSelected=" << spwSelected << endl;
463 0 : if (row == nRow)
464 0 : break;
465 :
466 : //should not be, but careful
467 0 : if (spwSelected == -1)
468 0 : spwSelected = 0;
469 :
470 : //cout << "row=" << row << " nRow=" << vb.nRow()
471 : // << " dId=" << dId << endl;
472 :
473 : //Double thisTime = curTime(row);
474 0 : Double thisTime = vb.time()(row);
475 0 : newTime = (thisTime - bufTime > aveTime - 0.005);
476 : //cout << "thisTime=" << std::setprecision(20) << thisTime
477 : // << " bufTime=" << bufTime << endl;
478 :
479 :
480 0 : newField = field - bufField;
481 0 : newArray = arry - bufArray;
482 :
483 0 : Int scan = vb.scan()(row);
484 0 : newScan = (scan - bufScan);
485 :
486 0 : newDesc = (dId - bufDesc);
487 :
488 : //if (newDesc)
489 : // cout << "newDesc=" << dId << endl;
490 :
491 0 : if (newField || (!crossScans && newScan) ||
492 0 : newTime || (!crossArrays && newArray) ||
493 0 : (!crossSpws && newDesc)) {
494 :
495 0 : if (nTotalTime > 0)
496 0 : sumTime /= nTotalTime;
497 : //cout << "-nTotalTime=" << nTotalTime
498 : // << " sumTime=" << sumTime << endl;
499 :
500 : //tmr2.mark();
501 0 : if (outputMode == MsAverager::TableMS)
502 0 : putAveTable(bufTime, bufField, bufScan, bufArray,
503 0 : bufDesc, *pAveBuff, nTime, sumTime);
504 : else
505 0 : putAveBuffer(bufTime, bufField, bufScan, bufArray,
506 0 : *pAveBuff, nTime);
507 : //usr += tmr2.user();
508 : //rea += tmr2.real();
509 : //sys += tmr2.system();
510 :
511 0 : if (pAveBuff) {
512 0 : delete pAveBuff;
513 0 : pAveBuff = 0;
514 : }
515 :
516 0 : bufTime = thisTime;
517 0 : bufField = field;
518 0 : bufScan = scan;
519 0 : bufArray = arry;
520 0 : bufDesc = dId;
521 :
522 0 : pAveBuff = new VisBuffer;
523 0 : initAveBuffer(bufTime, *pAveBuff, nAnt, nAveChan);
524 0 : nAveTime++;
525 0 : nTime = 0;
526 0 : sumTime = 0;
527 0 : nTotalTime = 0;
528 : }
529 :
530 0 : if (nextTime != thisTime) {
531 0 : sumTime += (thisTime - bufTime);
532 0 : nTotalTime++;
533 : }
534 :
535 : //cout << "row=" << row << " thisTime="
536 : // << std::setprecision(12) << thisTime
537 : // << " bufTime=" << bufTime
538 : // << " sumTime=" << sumTime
539 : // << " nTime=" << nTime << endl;
540 :
541 0 : Int ant1 = vb.antenna1()(row);
542 0 : Int ant2 = vb.antenna2()(row);
543 :
544 0 : Int orow = baselineRow(ant1, ant2);
545 0 : if (orow < 0)
546 0 : continue;
547 :
548 : //Int orow = baselineRow(nAnt, ant1, ant2);
549 : //cout << "row=" << row << " orow=" << orow
550 : // << " nAnt=" << nAnt
551 : // << " ant1=" << ant1 << " ant2=" << ant2 << endl;
552 : //aveBuff.antenna1()(orow) = ant1;
553 : //aveBuff.antenna2()(orow) = ant2;
554 :
555 0 : pAveBuff->time()(orow) = bufTime;
556 0 : pAveBuff->feed1()(orow) = vb.feed1()(row);
557 : //pAveBuff->feed2()(orow) = vb.feed2()(row);
558 :
559 0 : Complex xxx = 0.;
560 :
561 : //cout << "flagrow=" << vb.flagRow()(row) << endl;
562 : //showVisRow(vc, row);
563 : //cout << vb.flagCube();
564 :
565 0 : Int sw = SAS::spwIndexByDesc(dId, spw);
566 0 : if (sw < 0) sw = 0;
567 0 : pAveBuff->sigma()(orow) = spw(sw).rFreq;
568 :
569 0 : if (aveFlag) {
570 : //select flagged
571 :
572 : //init every cell unflagged
573 0 : for (Int pol = 0; pol < nAvePol; pol++) {
574 0 : for (Int chn = 0; chn < nAveChan; chn++) {
575 0 : aveV(pol, chn, 0) = 0.;
576 0 : aveF(pol, chn, 0) = 0;
577 : }
578 : }
579 :
580 0 : if (vb.flagRow()(row)) {
581 : //row flagged = everyone is flagged, should select
582 : //cout << "flagged row" << endl;
583 :
584 0 : for (Int pol = 0; pol < nAvePol; pol++) {
585 :
586 0 : xxx = 0.;
587 0 : Int p = 0;
588 0 : Int chn = 0;
589 :
590 0 : Int sChan = 0;
591 0 : for (chn = 0; chn < Int(spw(sw).chans.size()); chn++) {
592 :
593 0 : Int chNum = spw(sw).chans(chn);
594 0 : p++;
595 0 : if (aveMode == "SCALAR") {
596 0 : xxx += Complex(fabs(vc(pol, chNum, row)),0.0);
597 : }
598 : else {
599 0 : xxx += vc(pol, chNum, row);
600 : }
601 0 : if (p == aveChan) {
602 0 : aveV(pol, sChan, 0) = xxx / p;
603 0 : aveF(pol, sChan, 0) = 1;
604 0 : xxx = 0;
605 0 : p = 0;
606 0 : sChan++;
607 : }
608 : }
609 0 : if (p > 0) {
610 0 : aveV(pol, sChan, 0) = xxx / p;
611 0 : aveF(pol, sChan, 0) = 1;
612 : }
613 0 : sChan++;
614 : }
615 : }
616 : else {
617 : //only select flagged ;
618 : //cout << " not row flagged" << endl;
619 0 : for (Int pol = 0; pol < nAvePol; pol++) {
620 :
621 0 : xxx = 0.;
622 0 : Int p = 0;
623 0 : Int sx = 0;
624 0 : Int chn = 0;
625 :
626 0 : Int sChan = 0;
627 0 : for (chn = 0; chn < Int(spw(sw).chans.size()); chn++) {
628 :
629 0 : Int chNum = spw(sw).chans(chn);
630 0 : p++;
631 0 : if (vb.flagCube()(pol, chNum, row)) {
632 0 : if (aveMode == "SCALAR") {
633 0 : xxx += Complex(fabs(vc(pol, chNum, row)),0.0);
634 : }
635 : else {
636 0 : xxx += vc(pol, chNum, row);
637 : }
638 0 : sx++;
639 : }
640 0 : if (p == aveChan) {
641 0 : if (sx > 0) {
642 0 : aveV(pol, sChan, 0) = xxx / sx;
643 0 : aveF(pol, sChan, 0) = 1;
644 : }
645 0 : xxx = 0;
646 0 : p = 0;
647 0 : sx = 0;
648 0 : sChan++;
649 : }
650 : }
651 0 : if (p > 0) {
652 0 : if (sx > 0) {
653 0 : aveV(pol, sChan, 0) = xxx / sx;
654 0 : aveF(pol, sChan, 0) = 1;
655 : }
656 0 : sChan++;
657 : }
658 : }
659 : }
660 : //showVisRow(aveV, row);
661 : //cout << aveF.xyPlane(0) << endl;
662 : }
663 : else {
664 : //select non-flagged=good data
665 : // cout << "row=" << rowNumber(row)
666 : // << " flagrow=" << vb.flagRow()(row)
667 : // << " " << vb.flagCube();
668 :
669 : //init every cell flagged
670 0 : for (Int pol = 0; pol < nAvePol; pol++) {
671 0 : for (Int chn = 0; chn < nAveChan; chn++) {
672 0 : aveV(pol, chn, 0) = 0.;
673 0 : aveF(pol, chn, 0) = 1;
674 : }
675 : }
676 :
677 0 : if (!vb.flagRow()(row)) {
678 : //only select non-flagged
679 0 : for (Int pol = 0; pol < nAvePol; pol++) {
680 :
681 0 : xxx = 0.;
682 0 : Int p = 0;
683 0 : Int sx = 0;
684 :
685 0 : Int sChan = 0;
686 : Int chn;
687 0 : for (chn = 0; chn < Int(spw(sw).chans.size()); chn++){
688 :
689 0 : Int chNum = spw(sw).chans(chn);
690 0 : p++;
691 0 : if (!(vb.flagCube()(pol, chNum, row))) {
692 0 : if (aveMode == "SCALAR") {
693 0 : xxx += Complex(fabs(vc(pol, chNum, row)),0.0);
694 : }
695 : else {
696 0 : xxx += vc(pol, chNum, row);
697 : }
698 0 : sx++;
699 : }
700 0 : if (p == aveChan) {
701 0 : if (sx > 0) {
702 0 : aveV(pol, sChan, 0) = xxx / sx;
703 0 : aveF(pol, sChan, 0) = 0;
704 : }
705 0 : xxx = 0;
706 0 : p = 0;
707 0 : sx = 0;
708 0 : sChan++;
709 : }
710 : }
711 0 : if (p > 0) {
712 0 : if (sx > 0) {
713 0 : aveV(pol, sChan, 0) = xxx / sx;
714 0 : aveF(pol, sChan, 0) = 0;
715 : }
716 0 : sChan++;
717 : }
718 : }
719 : }
720 : }
721 :
722 : //cout << "orow=" << orow << endl;
723 : //weight is on polarization, should be a vector, but
724 : //defined in visBuffer as a scalar, good enugh for averaging
725 0 : Float wt = vb.weight()(row);
726 :
727 0 : pAveBuff->timeInterval()(orow) +=
728 0 : vb.timeInterval()(row) * wt;
729 0 : RigidVector<Double, 3> wtuvw = vb.uvw()(row) * Double(wt);
730 0 : pAveBuff->uvw()(orow) += wtuvw;
731 0 : pAveBuff->weight()(orow) += wt;
732 :
733 0 : for (Int pol = 0; pol < nAvePol; pol++) {
734 0 : for (Int chn = 0; chn < nAveChan; chn++) {
735 : //cout << "pol=" << pol << " chn=" << chn
736 : // << " wt=" << wt
737 : // << " cellFlag=" << aveF(pol, chn, 0)
738 : // << " aveFlag=" << aveFlag
739 : // << " Buff="
740 : // << pAveBuff->visCube()(pol, chn, orow);
741 0 : if (aveF(pol, chn, 0) == aveFlag) {
742 0 : pAveBuff->flagCube()(pol, chn, orow) = aveFlag;
743 : // DComplex wtvis =
744 : // aveV(pol, chn, 0) * Double(wt);
745 : //cout << " aveV=" << wt * aveV(pol, chn, 0) ;
746 0 : pAveBuff->visCube()(pol, chn, orow) +=
747 0 : wt * aveV(pol, chn, 0);
748 0 : pAveBuff->weightCube()(pol, chn, orow) += wt;
749 : }
750 : //cout << " ==>"
751 : // << pAveBuff->visCube()(pol, chn, orow);
752 :
753 : }
754 : }
755 : //showVisRow(pAveBuff->visCube(), orow);
756 : //showVisRow(aveV, 0);
757 : //cout << pAveBuff->flagCube().xyPlane(orow);
758 :
759 0 : if (crossBlines)
760 0 : aveRowMap(iRow, 0) = msRow;
761 : else
762 0 : aveRowMap(iRow, 0) = msRow + orow;
763 0 : aveRowMap(iRow, 1) = rowNumber(row);
764 0 : aveRowMap(iRow, 2) = spwSelected;
765 :
766 : //cout << "row=" << row << " "
767 : // << aveRowMap(iRow, 0) << " "
768 : // << aveRowMap(iRow, 1) << " "
769 : // << aveRowMap(iRow, 2) << endl;
770 :
771 0 : iRow++;
772 : //cout << "iRow=" << iRow << endl;
773 : }
774 0 : nTime++;
775 : //showVisRow(vc, 0);
776 0 : iIter++;
777 : //cout << "nTime=" << nTime << " iIter=" << iIter << endl;
778 0 : }
779 0 : iChunk++;
780 : }
781 0 : tmr3.show(" Averaging:");
782 :
783 0 : if (nTotalTime > 0)
784 0 : sumTime /= nTotalTime;
785 : //cout << "nTotalTime=" << nTotalTime
786 : // << " sumTime=" << sumTime << endl;
787 :
788 : //tmr2.mark();
789 0 : if (outputMode == MsAverager::TableMS)
790 0 : putAveTable(bufTime, bufField, bufScan, bufArray, bufDesc,
791 0 : *pAveBuff, nTime, sumTime);
792 : else
793 0 : putAveBuffer(bufTime, bufField, bufScan, bufArray,
794 0 : *pAveBuff, nTime);
795 : //usr += tmr2.user();
796 : //rea += tmr2.real();
797 : //sys += tmr2.system();
798 :
799 0 : if (pAveBuff) {
800 0 : delete pAveBuff;
801 0 : pAveBuff = 0;
802 : }
803 : //showAveMap(aveRowMap, aveChanMap);
804 :
805 : {
806 0 : ostringstream os;
807 : os << "TOTAL "
808 0 : << "Chunks=" << iChunk << " Iters=" << iIter
809 0 : << " Rows=" << iRow << " nAveTime=" << nAveTime
810 0 : << " nAveChan=" << nAveChan << " nAveRow=" << nAveRow
811 0 : << " nAvePol=" << nAvePol
812 0 : << endl;
813 : //SLog::slog()->out(os, fnname, clname, LogMessage::NORMAL5);
814 0 : }
815 :
816 : //cout << "iRow=" << iRow << endl;
817 : //cout << "pMS->nrow()=" << pMS->nrow() << endl;
818 0 : aveRowMap.resize(iRow, 3, true);
819 : //showAveMap(aveRowMap, aveChanMap);
820 :
821 : {
822 0 : ostringstream os;
823 0 : os << aveRowMap;
824 : //SLog::slog()->out(os, fnname, clname, LogMessage::DEBUG1);
825 0 : }
826 :
827 : //cout << "putAveTable:" << setw(11) << rea << " real "
828 : // << setw(11) << usr << " user "
829 : // << setw(11) << sys << " system" << endl;
830 :
831 0 : }
832 0 : catch (const AipsError &x) {
833 : //SLog::slog()->out(String("Error: ") + x.getMesg(),
834 : //fnname, clname, LogMessage::WARN);
835 0 : LogIO os(LogOrigin("MsAverager", "setAverager"));
836 : os << LogIO::WARN << "Error: "
837 : << x.getMesg()
838 0 : << LogIO::POST;
839 0 : aveOK = false;
840 0 : return ;
841 0 : }
842 : //catch (...) {cout << "problem------" << endl;}
843 0 : aveOK = true;
844 0 : return;
845 0 : }
846 :
847 0 : void MsAverager::putAveTable(Double bufTime, Int bufField, Int bufScan,
848 : Int bufArray, Int bufDesc, VisBuffer& aveBuff,
849 : Int nTime, Double timeShift) {
850 0 : if (nTime < 1)
851 0 : return;
852 :
853 0 : if (bufTime == -1 && bufField == -1 && bufScan == -1 && bufArray == -1)
854 0 : return;
855 :
856 :
857 : #if LOG2
858 : cout << " putAveTable: bufTime=" << bufTime
859 : << " bufField=" << bufField
860 : << " bufScan=" << bufScan
861 : << " bufArray=" << bufArray
862 : << " bufDesc=" << bufDesc
863 : << " nRow=" << aveBuff.nRow()
864 : << " nTime=" << nTime
865 : << " nChannel=" << aveBuff.nChannel()
866 : << endl;
867 : #endif
868 0 : if (nAveRow < 1)
869 0 : nAveRow = aveBuff.nRow();
870 : //cout << "nAveRow=" << nAveRow << endl;
871 : //ScalarColumn<Double> timeCol(aMS, "TIME");
872 :
873 : //sigma stores refFreq, verify it containing only one nozero value
874 0 : Float sigVal = 0;
875 0 : Int sigValChange = 0;
876 0 : for (Int row = 0; row < nAveRow; row++) {
877 0 : Float sig = aveBuff.sigma()(row);
878 0 : if (sig && sig != sigVal) {
879 0 : sigVal = sig;
880 0 : sigValChange++;
881 : }
882 : //cout << "sigma=" << aveBuff.sigma()(row)
883 : // << " row=" << row << endl;
884 : }
885 : //cout << "sigValChange=" << sigValChange << endl;
886 : //
887 0 : Vector<Float> refVal(1);
888 0 : refVal= sigVal;
889 :
890 :
891 0 : Vector<Complex> frqvel(nAveChan);
892 0 : Vector<Complex> sxsave(nAveChan);
893 0 : sxsave = 0;
894 0 : frqvel = 0;
895 :
896 0 : Int idx = SAS::spwIndexByDesc(bufDesc, spw);
897 0 : if (idx < 0) return;
898 :
899 : //calculate averaged velocities, store into velo
900 0 : Vector<Double> velo(nAveChan);
901 :
902 0 : if (aveVel)
903 0 : SAS::averageVelocity(sorryVel, pMS, spw, velo, idx, bufField,
904 0 : restfreq, frame, doppler);
905 : //cout << "velo=" << velo << " freq=" << spw(idx).aveFreqs << endl;
906 :
907 0 : for (uInt j = 0; j < spw(idx).aveFreqs.size(); j++) {
908 0 : frqvel(j) = Complex(velo(j), spw(idx).aveFreqs(j));
909 0 : sxsave(j) = Complex(spw(idx).aveChans(j), spw(idx).sxsChans(j));
910 : }
911 :
912 :
913 0 : if (crossBlines) {
914 : //cout << " crossBlines=true" << endl;
915 0 : aMS.addRow();
916 :
917 0 : Int tRow = msRow;
918 0 : msc->time().put(tRow, bufTime);
919 0 : msc->antenna1().put(tRow, aveBuff.antenna1()(0));
920 0 : msc->antenna2().put(tRow, aveBuff.antenna2()(0));
921 :
922 0 : msc->fieldId().put(tRow, bufField);
923 0 : msc->scanNumber().put(tRow, bufScan);
924 0 : msc->arrayId().put(tRow, bufArray);
925 0 : msc->feed1().put(tRow, aveBuff.feed1()(0));
926 0 : msc->dataDescId().put(tRow, bufDesc);
927 : //msc->feed2().put(tRow, aveBuff.feed2()(0));
928 0 : msc->sigma().put(tRow, refVal);
929 :
930 :
931 0 : Vector<Float> w(nAvePol < 1 ? 1 : nAvePol);
932 0 : for (Int row = 0; row < nAveRow; row++) {
933 0 : for (Int pol = 0; pol < nAvePol; pol++) {
934 0 : for (Int chn = 0; chn < nAveChan; chn++) {
935 0 : Float wc = aveBuff.weightCube()(pol, chn, row);
936 0 : if (wc > 0.) {
937 0 : aveBuff.visCube()(pol, chn, row) /= wc;
938 : }
939 : else {
940 0 : aveBuff.visCube()(pol, chn, row) = 0.;
941 : }
942 : }
943 : }
944 :
945 0 : Float wt = aveBuff.weight()(row);
946 0 : if (wt == 0.0f) {
947 0 : w = wt;
948 0 : wt = 1;
949 : }
950 : else {
951 0 : w = wt / nTime;
952 : }
953 0 : aveBuff.weight()(row) = wt;
954 0 : aveBuff.uvw()(row) = aveBuff.uvw()(row) * Double(1. / wt);
955 0 : aveBuff.timeInterval()(row) /= wt;
956 : }
957 :
958 0 : Cube<Complex> blV(nAvePol, nAveChan, 1);
959 0 : Cube<Complex> blU(nAvePol, nAveChan, 1);
960 0 : Cube<Bool> blF(nAvePol, nAveChan, 1);
961 0 : for (Int pol = 0; pol < nAvePol; pol++) {
962 0 : for (Int chn = 0; chn < nAveChan; chn++) {
963 0 : blV(pol, chn, 0) = 0.;
964 0 : blF(pol, chn, 0) = !aveFlag;
965 : }
966 : }
967 :
968 0 : for (Int pol = 0; pol < nAvePol; pol++) {
969 0 : for (Int chn = 0; chn < nAveChan; chn++) {
970 0 : Int blW = 0;
971 0 : for (Int row = 0; row < nAveRow; row++) {
972 0 : if (aveBuff.flagCube()(pol, chn, row) == aveFlag) {
973 0 : blV(pol, chn, 0) +=
974 0 : aveBuff.visCube()(pol, chn, row);
975 0 : blW++;
976 : }
977 : }
978 0 : if (blW > 0) {
979 0 : blV(pol, chn, 0) /= blW;
980 0 : blF(pol, chn, 0) = aveFlag;
981 : }
982 : }
983 : }
984 :
985 0 : Int nFlags = 0;
986 0 : for (Int pol = 0; pol < nAvePol; pol++) {
987 0 : for (Int chn = 0; chn < nAveChan; chn++) {
988 0 : if (!blF(pol, chn, 0)) {
989 0 : nFlags++;
990 : }
991 : }
992 : }
993 0 : if (nFlags)
994 0 : msc->flagRow().put(tRow, 0);
995 : else
996 0 : msc->flagRow().put(tRow, 1);
997 :
998 0 : msc->data().put(tRow, blV.xyPlane(0));
999 0 : msc->flag().put(tRow, blF.xyPlane(0));
1000 0 : RigidVector<Double, 3> bluvw = aveBuff.uvw()(0) * Double(1.);
1001 0 : msc->uvw().put(tRow, bluvw.vector());
1002 :
1003 : //showVisRow(blV, 0);
1004 : //showMsRow(msc, tRow);
1005 :
1006 0 : msc->weight().put(tRow, w);
1007 :
1008 0 : for (Int pol = 0; pol < nAvePol; pol++) {
1009 0 : for (Int chn = 0; chn < nAveChan; chn++) {
1010 0 : blV(pol, chn, 0) = frqvel(chn);
1011 0 : blU(pol, chn, 0) = sxsave(chn);
1012 : }
1013 : }
1014 : //msc->modelData().put(tRow, blV.xyPlane(0));
1015 : //msc->correctedData().put(tRow, blU.xyPlane(0));
1016 :
1017 0 : if (timeShift > 0.) {
1018 : //cout << " bufTime=" << std::setprecision(12)
1019 : // << bufTime << " timeShift=" << timeShift << endl;
1020 0 : msc->time().put(tRow, bufTime + timeShift);
1021 : }
1022 : else {
1023 0 : msc->time().put(tRow,
1024 0 : bufTime + 0.5 * aveBuff.timeInterval()(0) * (nTime - 1));
1025 : }
1026 0 : }
1027 : else {
1028 :
1029 : //cout << " crossBlines=false" << endl;
1030 :
1031 0 : for (Int row = 0; row < nAveRow; row++) {
1032 0 : for (Int pol = 0; pol < nAvePol; pol++) {
1033 0 : for (Int chn = 0; chn < nAveChan; chn++) {
1034 0 : aveBuff.modelVisCube()(pol, chn, row) = frqvel(chn);
1035 0 : aveBuff.correctedVisCube()(pol, chn, row) = sxsave(chn);
1036 : }
1037 : }
1038 : }
1039 :
1040 0 : for (Int row = 0; row < nAveRow; row++) {
1041 :
1042 0 : Int nFlags = 0;
1043 : //Int tFlags = nAvePol * nAveChan;
1044 0 : for (Int pol = 0; pol < nAvePol; pol++) {
1045 0 : for (Int chn = 0; chn < nAveChan; chn++) {
1046 0 : Float wc = aveBuff.weightCube()(pol, chn, row);
1047 : //flags should be correct, no need to fiddle
1048 0 : if (wc > 0.) {
1049 0 : aveBuff.visCube()(pol, chn, row) /= wc;
1050 : }
1051 : else {
1052 0 : aveBuff.visCube()(pol, chn, row) = 0.;
1053 : }
1054 0 : if (!aveBuff.flagCube()(pol, chn, row))
1055 0 : nFlags++;
1056 : //if (wc > 0)
1057 : // cout << " wc=" << wc << " row=" << row << " >>>"
1058 : // << aveBuff.visCube()(pol, chn, row);
1059 : }
1060 : }
1061 :
1062 0 : Int tRow = msRow + row;
1063 0 : aMS.addRow();
1064 :
1065 0 : msc->time().put(tRow, bufTime);
1066 0 : msc->antenna1().put(tRow, aveBuff.antenna1()(row));
1067 0 : msc->antenna2().put(tRow, aveBuff.antenna2()(row));
1068 : //cout << "sigma=" << aveBuff.sigma()(row) << " row=" << row << endl;
1069 0 : msc->sigma().put(tRow, refVal);
1070 :
1071 :
1072 0 : if (nFlags)
1073 0 : msc->flagRow().put(tRow, 0);
1074 : else
1075 0 : msc->flagRow().put(tRow, 1);
1076 :
1077 :
1078 0 : Vector<Float> w(nAvePol < 1 ? 1 : nAvePol);
1079 0 : Float wt = aveBuff.weight()(row);
1080 0 : if (wt == 0.0f) {
1081 0 : w = wt;
1082 0 : wt = 1;
1083 : }
1084 : else {
1085 0 : w = wt / nTime;
1086 : }
1087 :
1088 0 : msc->flag().put(tRow, aveBuff.flagCube().xyPlane(row));
1089 0 : msc->weight().put(tRow, w);
1090 0 : msc->uvw().put(tRow, (aveBuff.uvw()(row) * Double(1. / wt)).vector());
1091 :
1092 : //msc->data().put(tRow, aveBuff.visCube().xyPlane(row) / wt);
1093 0 : msc->data().put(tRow, aveBuff.visCube().xyPlane(row));
1094 0 : msc->fieldId().put(tRow, bufField);
1095 0 : msc->scanNumber().put(tRow, bufScan);
1096 0 : msc->arrayId().put(tRow, bufArray);
1097 0 : msc->feed1().put(tRow, aveBuff.feed1()(row));
1098 0 : msc->dataDescId().put(tRow, bufDesc);
1099 : //msc->feed2().put(tRow, aveBuff.feed2()(row));
1100 : //msc->modelData().put(tRow, aveBuff.visCube().xyPlane(row));
1101 :
1102 0 : if (timeShift > 0.) {
1103 : //cout << " bufTime=" << std::setprecision(12)
1104 : // << bufTime << " timeShift=" << timeShift << endl;
1105 0 : msc->time().put(tRow, bufTime + timeShift);
1106 : }
1107 : else {
1108 0 : Float itvl = aveBuff.timeInterval()(row) / wt;
1109 0 : msc->interval().put(tRow, itvl);
1110 : //cout << "interval=" << itvl << endl;
1111 0 : msc->time().put(tRow, bufTime + 0.5 * itvl * (nTime - 1));
1112 : }
1113 :
1114 : //showVisRow(aveBuff.visCube(), row);
1115 : //showMsRow(msc, tRow);
1116 :
1117 : // the flag and vis for the resulting row
1118 : // cout << "row=" << row
1119 : // << " flagCube=" << aveBuff.flagCube().xyPlane(row)
1120 : // << " visCube=" << aveBuff.visCube().xyPlane(row)
1121 : // << endl;
1122 : //msc->modelData().put(tRow, aveBuff.modelVisCube().xyPlane(row));
1123 : //msc->correctedData().put(tRow, aveBuff.correctedVisCube().xyPlane(row));
1124 0 : }
1125 :
1126 : }
1127 :
1128 : //MSSpWindowColumns spwColumn(aMS.spectralWindow());
1129 : //ScalarColumn<Int> numCol = spwColumn.numChan();
1130 : //spwColumn.numChan().put(idx, nAveChan);
1131 :
1132 0 : msRow = aMS.nrow();
1133 : //cout << "msRow=" << msRow << endl;
1134 :
1135 0 : return;
1136 0 : }
1137 :
1138 0 : void MsAverager::putAveBuffer(Double IfLog2 (bufTime), Int IfLog2 (bufField), Int IfLog2 (bufScan),
1139 : Int IfLog2 (bufArray), VisBuffer& aveBuff, Int nTime) {
1140 : #if LOG2
1141 : cout << " putAveBuffer: bufTime=" << bufTime
1142 : << " bufField=" << bufField
1143 : << " bufScan=" << bufScan
1144 : << " bufArray=" << bufArray
1145 : << " nRow=" << aveBuff.nRow()
1146 : << " nTime=" << nTime
1147 : << " nChannel=" << aveBuff.nChannel()
1148 : << endl;
1149 : #endif
1150 :
1151 0 : if (nTime < 1)
1152 0 : return;
1153 :
1154 0 : if (nAveRow < 1)
1155 0 : nAveRow = aveBuff.nRow();
1156 :
1157 0 : for (Int row = 0; row < nAveRow; row++) {
1158 :
1159 : //aveBuff.time()(row) = bufTime;
1160 :
1161 0 : Float wt = aveBuff.weight()(row);
1162 0 : if (wt == 0.0f) {
1163 0 : wt = 1.;
1164 : }
1165 :
1166 0 : aveBuff.timeInterval()(row) /= wt;
1167 0 : aveBuff.uvw()(row) *= 1.0f / wt;
1168 0 : for (Int pol = 0; pol < nAvePol; pol++) {
1169 0 : for (Int chn = 0; chn < nAveChan; chn++) {
1170 0 : aveBuff.visCube()(pol, chn, row) *= 1.0f / wt;
1171 : }
1172 : }
1173 : }
1174 :
1175 :
1176 0 : aveList.push_back(pAveBuff);
1177 :
1178 0 : return;
1179 : }
1180 :
1181 0 : Int MsAverager::baselineRow(const Int& nAnt, const Int& a1, const Int& a2)
1182 : {
1183 : Int index;
1184 0 : index = nAnt * a1 - (a1 * (a1 - 1)) / 2 + a2 - a1;
1185 0 : if (a1 > a2)
1186 0 : index = nAnt * a2 - (a2 * (a2 - 1)) / 2 + a1 - a2;
1187 0 : return index;
1188 : }
1189 :
1190 0 : Int MsAverager::baselineRow(const Int& a1, const Int& a2)
1191 : {
1192 : //cout << "bls=" << bls << " bls.shape=" << bls.shape()
1193 : // << " a1=" << a1 << " a2=" << a2 << endl;
1194 0 : if (a1 < 0 && a2 < 0)
1195 0 : return baselines.shape()[0];
1196 0 : Int index = -1;
1197 0 : for (Int i = 0; i < baselines.shape()[0]; i++) {
1198 0 : if (a1 == baselines(i, 0) && a2 == baselines(i, 1))
1199 0 : index = i;
1200 : }
1201 0 : return index;
1202 : }
1203 :
1204 0 : void MsAverager::showMsRow(MSMainColumns* msc, Int row) {
1205 0 : cout << "row=" << row
1206 0 : << "\ntime=" << std::setprecision(12) << msc->time()(row)
1207 0 : << " rowflag=" << msc->flagRow()(row)
1208 : << std::setprecision(5)
1209 0 : << "\nant1=" << msc->antenna1()(row)
1210 0 : << " ant2=" << msc->antenna2()(row)
1211 0 : << "\ndata=" << msc->data()(row)
1212 0 : << "\nflag=" << msc->flag()(row)
1213 0 : << "\nsigma=" << msc->sigma()(row)
1214 0 : << "\ndesc=" << msc->dataDescId()(row)
1215 0 : << endl;
1216 0 : }
1217 :
1218 0 : void MsAverager::showVisRow(Cube<Complex>& vc, Int row) {
1219 0 : IPosition ipos = vc.shape();
1220 0 : Int nPol = ipos(0);
1221 0 : Int nChan = ipos(1);
1222 : cout << std::setprecision(8)
1223 0 : << " vis (" << nPol << " pol x " << nChan << " chan):\n";
1224 0 : for (Int chn = 0; chn < nChan; chn++) {
1225 0 : for (Int pol = 0; pol < nPol; pol++) {
1226 0 : cout << " " << vc(pol, chn, row);
1227 : }
1228 0 : cout << endl;
1229 : }
1230 0 : }
1231 :
1232 0 : void MsAverager::showColumnNames() {
1233 0 : cout << pMS->tableDesc().columnNames() << endl;
1234 : //UVW, FLAG, FLAG_CATEGORY, WEIGHT, SIGMA, ANTENNA1, ANTENNA2,
1235 : //ARRAY_ID, DATA_DESC_ID, EXPOSURE, FEED1, FEED2, FIELD_ID,
1236 : //FLAG_ROW, INTERVAL, OBSERVATION_ID, PROCESSOR_ID, SCAN_NUMBER,
1237 : //STATE_ID, TIME, TIME_CENTROID, DATA, WEIGHT_SPECTRUM,
1238 : //MODEL_DATA, CORRECTED_DATA
1239 0 : }
1240 :
1241 0 : void MsAverager::showAveMap(Matrix<Int> &/*rMap*/, Matrix<Int> &/*cMap*/) {
1242 : //cout << "aveRowMap=" << std::setprecision(8) << rMap;
1243 : //cout << "aveChanMap=" << std::setprecision(12) << cMap;
1244 0 : }
1245 :
1246 0 : Bool MsAverager::hasColumn(casacore::String const& col) {
1247 0 : Vector<String> cols = pMS->tableDesc().columnNames();
1248 0 : for (uInt i = 0; i < cols.nelements(); i++) {
1249 0 : if (cols(i) == col)
1250 0 : return true;
1251 : }
1252 0 : LogIO os(LogOrigin("MsAverager", "hasColumn"));
1253 0 : os << LogIO::WARN << String("No column '") + col + "' in the MS"
1254 0 : << LogIO::POST;
1255 :
1256 0 : return false;
1257 0 : }
1258 :
1259 0 : Bool MsAverager::isDataColumn(casacore::String const& col) {
1260 0 : return col == "DATA"
1261 0 : || col == "CORRECTEDDATA"
1262 0 : || col == "MODELDATA"
1263 : ;
1264 : }
1265 :
1266 0 : void MsAverager::getMap(Matrix<Int>& rowMap, Matrix<Int>& chanMap) {
1267 0 : rowMap.resize(0, 0);
1268 0 : rowMap = aveRowMap;
1269 :
1270 0 : chanMap.resize(0, 0);
1271 0 : chanMap = aveChanMap;
1272 0 : }
1273 :
1274 0 : void MsAverager::getMS(MS& ms) {
1275 : //cout << "getMS aMS=" << aMS << endl;
1276 0 : if (outputMode != MsAverager::TableMS) {
1277 : //SLog::slog()->out(String("MS is not available in 'ListBuffer' mode"),
1278 : //"getMS", clname, LogMessage::WARN);
1279 0 : LogIO os(LogOrigin("MsAverager", "getMS"));
1280 0 : os << LogIO::WARN << String("MS is not available in 'ListBuffer' mode")
1281 0 : << LogIO::POST;
1282 0 : return;
1283 0 : }
1284 : //ms = aMS;
1285 0 : ms = MS(aMS);
1286 : }
1287 :
1288 0 : void MsAverager::getXY(Vector<Double>& x, casacore::Vector<Double>& y,
1289 : Vector<Int>& f, Int pol) {
1290 :
1291 0 : if (outputMode != MsAverager::ListBuffer) {
1292 : //SLog::slog()->out(String("MS is not available in 'TableMS' mode"),
1293 : //"getXY", clname, LogMessage::WARN);
1294 0 : LogIO os(LogOrigin("MsAverager", "getXY"));
1295 0 : os << LogIO::WARN << String("MS is not available in 'TableMS' mode")
1296 0 : << LogIO::POST;
1297 0 : return;
1298 0 : }
1299 :
1300 0 : VisBuffer* pAveBuff = aveList.front( );
1301 0 : IPosition ip = pAveBuff->visCube().shape();
1302 :
1303 0 : Int nBuff = aveList.size( ); // nAveTime
1304 0 : Int nRow = ip(2);
1305 : #if LOG2
1306 : Int nChan = ip(1); // nAveChan
1307 : Int nPol = ip(0); // nAvePol
1308 :
1309 : cout << "getXY nBuff=" << nBuff << " nPol=" << nPol
1310 : << " nChan=" << nChan << " nRow=" << nRow
1311 : << " len=" << nChan * nRow << endl;
1312 : #endif
1313 0 : Int len = nBuff * nRow * nAveChan;
1314 :
1315 0 : x.resize(len);
1316 0 : y.resize(len);
1317 0 : f.resize(len);
1318 :
1319 0 : int i = 0;
1320 0 : for ( auto pb : aveList ) {
1321 0 : if (pb != 0) {
1322 0 : for (int row = 0; row < nRow; row++) {
1323 0 : for (Int chn = 0; chn < nAveChan; chn++) {
1324 0 : Complex xxx = pb->visCube()(pol, chn, row);
1325 0 : y(i) = xxx.real();
1326 0 : x(i) = pb->time()(row);
1327 0 : f(i) = pb->flagCube()(pol, chn, row);
1328 0 : i++;
1329 : }
1330 : //showVisRow(aveV, 0);
1331 : }
1332 : }
1333 : }
1334 :
1335 0 : }
1336 :
1337 0 : void MsAverager::initAveBuffer(Double bufTime,
1338 : VisBuffer& aveBuff, Int /*nAnt*/, Int nChan)
1339 : {
1340 :
1341 :
1342 0 : Int nRowAdd = baselineRow ();
1343 : //Int nRowAdd = baselineRow (nAnt, nAnt - 1, nAnt - 1) + 1;
1344 0 : aveBuff.nRow();
1345 0 : aveBuff.nRow() = nRowAdd;
1346 0 : aveBuff.nChannel() = nChan;
1347 : //cout << "nRowAdd=" << nRowAdd << endl;
1348 : #if LOG2
1349 : cout << "initAveBuffer: bufTime=" << std::setprecision(12) << bufTime
1350 : << " nRow=" << aveBuff.nRow()
1351 : << " nChan=" << aveBuff.nChannel()
1352 : << endl;
1353 : #endif
1354 :
1355 0 : Int nRow = aveBuff.nRow();
1356 0 : aveBuff.antenna1().resize(nRow);
1357 0 : aveBuff.antenna2().resize(nRow);
1358 0 : aveBuff.time().resize(nRow);
1359 0 : aveBuff.timeInterval().resize(nRow);
1360 0 : aveBuff.uvw().resize(nRow);
1361 : //aveBuff.visibility().resize(nChan, nRow);
1362 : //aveBuff.flag().resize(nChan, nRow);
1363 0 : aveBuff.visCube().resize(nAvePol, nChan, nRow);
1364 0 : aveBuff.flagCube().resize(nAvePol, nChan, nRow);
1365 0 : aveBuff.weightCube().resize(nAvePol, nChan, nRow);
1366 0 : aveBuff.modelVisCube().resize(nAvePol, nChan, nRow);
1367 0 : aveBuff.correctedVisCube().resize(nAvePol, nChan, nRow);
1368 0 : aveBuff.weight().resize(nRow);
1369 0 : aveBuff.flagRow().resize(nRow);
1370 0 : aveBuff.feed1().resize(nRow);
1371 : //aveBuff.feed2().resize(nRow);
1372 0 : aveBuff.sigma().resize(nRow);
1373 :
1374 0 : Int row = 0;
1375 : //for (Int ant1 = 0; ant1 < nAnt; ant1++) {
1376 : // for (Int ant2 = ant1; ant2 < nAnt; ant2++) {
1377 : // aveBuff.antenna1()(row) = ant1;
1378 : // aveBuff.antenna2()(row) = ant2;
1379 : // row++;
1380 : // }
1381 : //}
1382 :
1383 0 : for (row = 0; row < nRow; row++) {
1384 0 : aveBuff.antenna1()(row) = baselines(row, 0);
1385 0 : aveBuff.antenna2()(row) = baselines(row, 1);
1386 0 : aveBuff.time()(row) = bufTime;
1387 0 : aveBuff.timeInterval()(row) = 0.0;
1388 0 : aveBuff.uvw()(row) = 0.0;
1389 0 : aveBuff.feed1()(row) = 0;
1390 0 : aveBuff.sigma()(row) = 0;
1391 : //aveBuff.feed2()(row) = 0;
1392 : //for (Int chn = 0; chn < nChan; chn++) {
1393 : // aveBuff.visibility()(chn, row) = CStokesVector();
1394 : // aveBuff.flag()(chn, row) = true;
1395 : //};
1396 0 : for (Int pol = 0; pol < nAvePol; pol++) {
1397 0 : for (Int chn = 0; chn < nChan; chn++) {
1398 0 : aveBuff.flagCube()(pol, chn, row) = !aveFlag;
1399 0 : aveBuff.visCube()(pol, chn, row) = 0;
1400 0 : aveBuff.weightCube()(pol, chn, row) = 0.;
1401 0 : aveBuff.modelVisCube()(pol, chn, row) = 0.;
1402 0 : aveBuff.correctedVisCube()(pol, chn, row) = 0.;
1403 : }
1404 : };
1405 0 : aveBuff.weight()(row) = 0.0f;
1406 0 : aveBuff.flagRow()(row) = !aveFlag;
1407 : };
1408 : //cout << "initAve: nRow=" << aveBuff.nRow()
1409 : // << " nChannel=" << aveBuff->nChannel() << endl;
1410 :
1411 :
1412 :
1413 0 : }
1414 :
1415 :
1416 : }
1417 :
|
