Line data Source code
1 : //# VisSet.cc: Implementation of VisSet
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2005
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 <msvis/MSVis/VisSet.h>
29 : #include <msvis/MSVis/VisSetUtil.h>
30 : #include <msvis/MSVis/VisBuffer.h>
31 : #include <msvis/MSVis/VisModelDataI.h>
32 : #include <casacore/ms/MeasurementSets/MSColumns.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/casa/Arrays/ArrayLogical.h>
35 : #include <casacore/casa/Arrays/ArrayUtil.h>
36 : #include <casacore/casa/Arrays/Matrix.h>
37 : #include <casacore/casa/Exceptions/Error.h>
38 : #include <casacore/casa/Containers/Record.h>
39 : #include <casacore/tables/Tables/ArrColDesc.h>
40 : #include <casacore/tables/Tables/ScaColDesc.h>
41 : #include <casacore/tables/Tables/TableDesc.h>
42 : #include <casacore/tables/Tables/TableRecord.h>
43 : #include <casacore/tables/DataMan/TiledDataStMan.h>
44 : #include <casacore/tables/DataMan/TiledShapeStMan.h>
45 : #include <casacore/tables/DataMan/StandardStMan.h>
46 : #include <casacore/tables/DataMan/TiledDataStManAccessor.h>
47 : #include <casacore/tables/Tables/TableIter.h>
48 : #include <casacore/tables/DataMan/CompressComplex.h>
49 : #include <casacore/tables/DataMan/CompressFloat.h>
50 : #include <casacore/casa/Arrays/Slice.h>
51 : #include <casacore/casa/Arrays/Slicer.h>
52 : #include <casacore/casa/Utilities/GenSort.h>
53 : #include <iostream>
54 :
55 : #include <casacore/casa/Logging/LogMessage.h>
56 : #include <casacore/casa/Logging/LogSink.h>
57 :
58 :
59 : using namespace casacore;
60 : namespace casa { //# NAMESPACE CASA - BEGIN
61 :
62 86 : VisSet::VisSet(MeasurementSet& ms,const Block<Int>& columns,
63 : const Matrix<Int>& chanSelection, Double timeInterval,
64 86 : Bool compress, Bool doModelData)
65 86 : :ms_p(ms)
66 : {
67 172 : LogSink logSink;
68 258 : LogMessage message(LogOrigin("VisSet","VisSet"));
69 :
70 86 : blockOfMS_p= new Block<MeasurementSet> ();
71 86 : blockOfMS_p->resize(1);
72 86 : (*blockOfMS_p)[0]=ms_p;
73 86 : multims_p=false;
74 : // sort out the channel selection
75 86 : Int nSpw=ms_p.spectralWindow().nrow();
76 172 : MSSpWindowColumns msSpW(ms_p.spectralWindow());
77 86 : if(nSpw==0)
78 0 : throw(AipsError(String("There is no valid spectral windows in "+ms_p.tableName())));
79 86 : selection_p.resize(2,nSpw);
80 : // fill in default selection
81 86 : selection_p.row(0)=0; //start
82 86 : selection_p.row(1)=msSpW.numChan().getColumn();
83 86 : for (uInt i=0; i<chanSelection.ncolumn(); i++) {
84 0 : Int spw=chanSelection(2,i);
85 0 : if (spw>=0 && spw<nSpw && chanSelection(0,i)>=0 &&
86 0 : chanSelection(0,i)+chanSelection(1,i)<=selection_p(1,spw)) {
87 : // looks like a valid selection, implement it
88 0 : selection_p(0,spw)=chanSelection(0,i);
89 0 : selection_p(1,spw)=chanSelection(1,i);
90 : }
91 : }
92 :
93 86 : Bool init=true;
94 86 : if (ms.tableDesc().isColumn("CORRECTED_DATA")) {
95 86 : init=false;
96 : }
97 : // in case model data exists and the user do not want it anymore
98 86 : if(ms.tableDesc().isColumn("MODEL_DATA") && !doModelData){
99 0 : init=true;
100 : }
101 : // Add scratch columns
102 86 : if (init) {
103 :
104 0 : VisSetUtil::removeCalSet(ms);
105 0 : addCalSet2(ms, compress, doModelData);
106 :
107 :
108 : // Force re-sort (in VisIter ctor below) by deleting current sort info
109 0 : if (ms.keywordSet().isDefined("SORT_COLUMNS"))
110 0 : ms.rwKeywordSet().removeField("SORT_COLUMNS");
111 0 : if (ms.keywordSet().isDefined("SORTED_TABLE"))
112 0 : ms.rwKeywordSet().removeField("SORTED_TABLE");
113 : }
114 :
115 :
116 86 : iter_p=new VisIter(ms_p,columns,timeInterval);
117 172 : for (uInt spw=0; spw<selection_p.ncolumn(); spw++) {
118 86 : iter_p->selectChannel(1,selection_p(0,spw),selection_p(1,spw),0,spw);
119 : }
120 :
121 : // Now initialize the scratch columns
122 86 : if (init)
123 0 : initCalSet(0);
124 :
125 86 : }
126 :
127 315 : VisSet::VisSet(MeasurementSet& ms,const Block<Int>& columns,
128 : const Matrix<Int>& chanSelection,
129 : Bool addScratch,
130 315 : Double timeInterval, Bool compress, Bool doModelData)
131 315 : :ms_p(ms)
132 : {
133 630 : LogSink logSink;
134 945 : LogMessage message(LogOrigin("VisSet","VisSet"));
135 :
136 315 : blockOfMS_p= new Block<MeasurementSet> ();
137 315 : blockOfMS_p->resize(1);
138 315 : (*blockOfMS_p)[0]=ms_p;
139 315 : multims_p=false;
140 : // sort out the channel selection
141 315 : Int nSpw=ms_p.spectralWindow().nrow();
142 630 : MSSpWindowColumns msSpW(ms_p.spectralWindow());
143 315 : if(nSpw==0)
144 0 : throw(AipsError(String("There is no valid spectral windows in "+ms_p.tableName())));
145 315 : selection_p.resize(2,nSpw);
146 : // fill in default selection
147 315 : selection_p.row(0)=0; //start
148 315 : selection_p.row(1)=msSpW.numChan().getColumn();
149 315 : for (uInt i=0; i<chanSelection.ncolumn(); i++) {
150 0 : Int spw=chanSelection(2,i);
151 0 : if (spw>=0 && spw<nSpw && chanSelection(0,i)>=0 &&
152 0 : chanSelection(0,i)+chanSelection(1,i)<=selection_p(1,spw)) {
153 : // looks like a valid selection, implement it
154 0 : selection_p(0,spw)=chanSelection(0,i);
155 0 : selection_p(1,spw)=chanSelection(1,i);
156 : }
157 : }
158 :
159 315 : Bool init=true;
160 :
161 315 : if (ms.tableDesc().isColumn("CORRECTED_DATA")) {
162 243 : init=false;
163 : }
164 : // in case model data exists and the user do not want it anymore
165 315 : if(ms.tableDesc().isColumn("MODEL_DATA") && !doModelData){
166 77 : init=true;
167 : }
168 : // cout << boolalpha << "addScratch = " << addScratch << endl;
169 :
170 : // Add scratch columns
171 315 : if (addScratch && init) {
172 :
173 0 : VisSetUtil::removeCalSet(ms);
174 0 : addCalSet2(ms, compress, doModelData);
175 :
176 :
177 : // Force re-sort (in VisIter ctor below) by deleting current sort info
178 0 : if (ms.keywordSet().isDefined("SORT_COLUMNS"))
179 0 : ms.rwKeywordSet().removeField("SORT_COLUMNS");
180 0 : if (ms.keywordSet().isDefined("SORTED_TABLE"))
181 0 : ms.rwKeywordSet().removeField("SORTED_TABLE");
182 : }
183 :
184 : // Generate the VisIter
185 315 : iter_p=new VisIter(ms_p,columns,timeInterval);
186 3806 : for (uInt spw=0; spw<selection_p.ncolumn(); spw++) {
187 3491 : iter_p->selectChannel(1,selection_p(0,spw),selection_p(1,spw),0,spw);
188 : }
189 :
190 : // Now initialize the scratch columns
191 315 : if (addScratch && init)
192 0 : initCalSet(0);
193 :
194 315 : }
195 :
196 0 : VisSet::VisSet(Block<MeasurementSet>& mss,const Block<Int>& columns,
197 : const Block< Matrix<Int> >& chanSelection, Bool addScratch,
198 : Double timeInterval,
199 0 : Bool compress, Bool addModelData)
200 : {
201 :
202 0 : multims_p=true;
203 0 : blockOfMS_p = &mss;
204 0 : Int numMS=mss.nelements();
205 0 : ms_p=mss[numMS-1];
206 :
207 0 : Block<Vector<Int> > blockNGroup(numMS);
208 0 : Block<Vector<Int> > blockStart(numMS);
209 0 : Block<Vector<Int> > blockWidth(numMS);
210 0 : Block<Vector<Int> > blockIncr(numMS);
211 0 : Block<Vector<Int> > blockSpw(numMS);
212 :
213 0 : for (Int k=0; k < numMS; ++k){
214 : // sort out the channel selection
215 0 : Int nSpw=mss[k].spectralWindow().nrow();
216 0 : if(nSpw==0)
217 0 : throw(AipsError(String("There is no valid spectral windows in "+mss[k].tableName())));
218 0 : blockNGroup[k]=Vector<Int> (nSpw,1);
219 0 : blockIncr[k]=Vector<Int> (nSpw,1);
220 : // At this stage select all spw
221 0 : blockSpw[k].resize(nSpw);
222 0 : indgen(blockSpw[k]);
223 0 : if(chanSelection[k].nelements()!=0){
224 0 : blockStart[k]=chanSelection[k].row(0);
225 0 : blockWidth[k]=chanSelection[k].row(1);
226 : }
227 0 : MSSpWindowColumns msSpW(mss[k].spectralWindow());
228 0 : selection_p.resize(2,nSpw);
229 : //Drat...need to figure this one out....
230 : // fill in default selection
231 0 : selection_p.row(0)=0; //start
232 0 : selection_p.row(1)=msSpW.numChan().getColumn();
233 0 : blockStart[k].resize(selection_p.row(0).nelements());
234 0 : blockStart[k]=selection_p.row(0);
235 0 : blockWidth[k].resize(selection_p.row(1).nelements());
236 0 : blockWidth[k]=selection_p.row(1);
237 0 : for (uInt i=0; i<chanSelection[k].ncolumn(); i++) {
238 0 : Int spw=chanSelection[k](2,i);
239 0 : if (spw>=0 && spw<nSpw && chanSelection[k](0,i)>=0 &&
240 0 : chanSelection[k](0,i)+chanSelection[k](1,i)<=selection_p(1,spw)) {
241 : // looks like a valid selection, implement it
242 0 : selection_p(0,spw)=chanSelection[k](0,i);
243 0 : selection_p(1,spw)=chanSelection[k](1,i);
244 0 : blockStart[k][spw]=chanSelection[k](0,i);
245 0 : blockWidth[k][spw]=chanSelection[k](1,i);
246 : }
247 : }
248 : }
249 0 : iter_p=new VisIter(mss,columns,timeInterval);
250 :
251 0 : iter_p->selectChannel(blockNGroup, blockStart, blockWidth, blockIncr,
252 0 : blockSpw);
253 :
254 0 : if(addScratch){
255 0 : for (Int k=0; k < numMS ; ++k){
256 0 : addScratchCols(mss[k], compress, addModelData);
257 : }
258 : }
259 :
260 0 : }
261 :
262 :
263 0 : VisSet::VisSet(ROVisibilityIterator& vi){
264 :
265 : ///Valid for single ms iterators
266 0 : ms_p=vi.ms();
267 0 : blockOfMS_p= new Block<MeasurementSet> ();
268 0 : blockOfMS_p->resize(1);
269 0 : (*blockOfMS_p)[0]=ms_p;
270 0 : multims_p=false;
271 0 : Block<Int> columns(0);
272 0 : Double timeInterval=0.0;
273 0 : iter_p=new VisIter(ms_p,columns,timeInterval);
274 0 : Block<Vector<Int> > ngroup;
275 0 : Block<Vector<Int> > start;
276 0 : Block<Vector<Int> > width;
277 0 : Block<Vector<Int> > incr;
278 0 : Block<Vector<Int> > spws;
279 0 : vi.getChannelSelection(ngroup,start,width,incr,spws);
280 0 : iter_p->selectChannel(ngroup,start,width,incr,spws);
281 :
282 0 : }
283 :
284 0 : VisSet::VisSet(MeasurementSet& ms, const Matrix<Int>& chanSelection,
285 0 : Double timeInterval)
286 0 : :ms_p(ms)
287 : {
288 0 : LogSink logSink;
289 0 : LogMessage message(LogOrigin("VisSet","VisSet"));
290 :
291 0 : blockOfMS_p= new Block<MeasurementSet> ();
292 0 : blockOfMS_p->resize(1);
293 0 : (*blockOfMS_p)[0]=ms_p;
294 0 : multims_p=false;
295 : // sort out the channel selection
296 0 : Int nSpw=ms_p.spectralWindow().nrow();
297 0 : if(nSpw==0)
298 0 : throw(AipsError(String("There is no valid spectral windows in "+ms_p.tableName())));
299 0 : MSSpWindowColumns msSpW(ms_p.spectralWindow());
300 0 : selection_p.resize(2,nSpw);
301 : // fill in default selection
302 0 : selection_p.row(0)=0; //start
303 0 : selection_p.row(1)=msSpW.numChan().getColumn();
304 0 : for (uInt i=0; i<chanSelection.ncolumn(); i++) {
305 0 : Int spw=chanSelection(2,i);
306 0 : if (spw>=0 && spw<nSpw && chanSelection(0,i)>=0 &&
307 0 : chanSelection(0,i)+chanSelection(1,i)<=selection_p(1,spw)) {
308 : // looks like a valid selection, implement it
309 0 : selection_p(0,spw)=chanSelection(0,i);
310 0 : selection_p(1,spw)=chanSelection(1,i);
311 : }
312 : }
313 :
314 :
315 :
316 0 : Block<Int> columns(0);
317 :
318 0 : iter_p=new VisIter(ms_p,columns,timeInterval);
319 0 : for (uInt spw=0; spw<selection_p.ncolumn(); spw++) {
320 0 : iter_p->selectChannel(1,selection_p(0,spw),selection_p(1,spw),0,spw);
321 : }
322 :
323 :
324 0 : }
325 :
326 24 : VisSet::VisSet(const VisSet& vs,const Block<Int>& columns,
327 24 : Double timeInterval)
328 : {
329 24 : ms_p=vs.ms_p;
330 24 : blockOfMS_p=new Block<MeasurementSet>();
331 24 : blockOfMS_p->resize(1);
332 24 : (*blockOfMS_p)[0]=ms_p;
333 24 : multims_p=false;
334 24 : selection_p.resize(vs.selection_p.shape());
335 24 : selection_p=vs.selection_p;
336 :
337 24 : iter_p=new VisIter(ms_p,columns,timeInterval);
338 50 : for (uInt spw=0; spw<selection_p.ncolumn(); spw++) {
339 26 : iter_p->selectChannel(1,selection_p(0,spw),selection_p(1,spw),0,spw);
340 : }
341 24 : }
342 :
343 0 : VisSet& VisSet::operator=(const VisSet& other)
344 : {
345 0 : if (this == &other) return *this;
346 0 : ms_p=other.ms_p;
347 :
348 0 : blockOfMS_p= new Block<MeasurementSet> ();
349 0 : blockOfMS_p->resize(other.blockOfMS_p->nelements());
350 0 : for (uInt k=0; k < blockOfMS_p->nelements() ; ++k)
351 0 : (*blockOfMS_p)[k]= (*(other.blockOfMS_p))[k];
352 :
353 0 : selection_p.resize(other.selection_p.shape());
354 0 : selection_p=other.selection_p;
355 0 : *iter_p=*(other.iter_p);
356 0 : return *this;
357 : }
358 :
359 425 : VisSet::~VisSet() {
360 425 : ms_p.flush();
361 425 : ms_p=MeasurementSet();// breaking reference
362 :
363 425 : delete iter_p; iter_p=0;
364 : //Only one ms then most probably did not use the multi ms constructor
365 425 : if(!multims_p){
366 425 : (*blockOfMS_p)[0]=MeasurementSet();//breaking reference
367 : }
368 425 : delete blockOfMS_p;
369 425 : };
370 :
371 :
372 : void
373 94 : VisSet::resetVisIter(const Block<Int>& columns, Double timeInterval,
374 : asyncio::PrefetchColumns * prefetchColumns)
375 : {
376 94 : if (iter_p){
377 94 : delete iter_p; // Delete existing VisIter:
378 : }
379 :
380 : // Make new VisIter with existing ms_p, and new sort/interval
381 :
382 94 : iter_p = new VisIter (prefetchColumns, ms_p, columns, true, timeInterval);
383 :
384 : // Inform new VisIter of channel selection
385 :
386 279 : for (uInt spw=0; spw<selection_p.ncolumn(); spw++) {
387 185 : iter_p->selectChannel(1,selection_p(0,spw),selection_p(1,spw),0,spw);
388 : }
389 94 : }
390 :
391 11 : void VisSet::initCalSet(Int /*calSet*/)
392 : {
393 22 : LogSink logSink;
394 33 : LogMessage message(LogOrigin("VisSet","VisSet"));
395 11 : Bool doModel=iter_p->ms().tableDesc().isColumn("MODEL_DATA");
396 11 : doModel ? message.message("Initializing MODEL_DATA (to unity) and CORRECTED_DATA (to DATA)"):
397 11 : message.message("Initializing CORRECTED_DATA (to DATA)");
398 11 : logSink.post(message);
399 :
400 22 : Vector<Int> lastCorrType;
401 22 : Vector<Bool> zero;
402 11 : Int nRows(0);
403 11 : iter_p->setRowBlocking(10000);
404 190 : for (iter_p->originChunks(); iter_p->moreChunks(); iter_p->nextChunk()) {
405 :
406 : // figure out which correlations to set to 1. and 0. for the model.
407 358 : Vector<Int> corrType; iter_p->corrType(corrType);
408 179 : uInt nCorr = corrType.nelements();
409 347 : if (nCorr != lastCorrType.nelements() ||
410 168 : !allEQ(corrType, lastCorrType)) {
411 :
412 11 : lastCorrType.resize(nCorr);
413 11 : lastCorrType=corrType;
414 11 : zero.resize(nCorr);
415 :
416 49 : for (uInt i=0; i<nCorr; i++)
417 : {
418 114 : zero[i]=(corrType[i]==Stokes::RL || corrType[i]==Stokes::LR ||
419 76 : corrType[i]==Stokes::XY || corrType[i]==Stokes::YX);
420 : }
421 : }
422 358 : for (iter_p->origin(); iter_p->more(); (*iter_p)++) {
423 179 : nRows+=iter_p->nRow();
424 :
425 : // Read DATA
426 358 : Cube<Complex> data;
427 179 : iter_p->visibility(data, VisibilityIterator::Observed);
428 :
429 : // Set CORRECTED_DATA
430 179 : iter_p->setVis(data,VisibilityIterator::Corrected);
431 :
432 179 : if(doModel){
433 : // Set MODEL_DATA
434 0 : data.set(1.0);
435 0 : if (ntrue(zero)>0) {
436 0 : for (uInt i=0; i < nCorr; i++) {
437 0 : if (zero[i]) data(Slice(i), Slice(), Slice()) = Complex(0.0,0.0);
438 : }
439 : }
440 0 : iter_p->setVis(data,VisibilityIterator::Model);
441 : }
442 : }
443 : }
444 11 : flush();
445 :
446 11 : iter_p->originChunks();
447 11 : iter_p->setRowBlocking(0);
448 :
449 22 : ostringstream os;
450 11 : os << "Initialized " << nRows << " rows.";
451 11 : message.message(os.str());
452 11 : logSink.post(message);
453 :
454 11 : }
455 :
456 31 : void VisSet::flush() {
457 31 : if(iter_p->newMS()){
458 0 : Int msId = iter_p->msId();
459 :
460 0 : ms_p=(*blockOfMS_p)[msId];
461 : }
462 31 : ms_p.flush();
463 31 : };
464 :
465 1756 : VisIter& VisSet::iter() { return *iter_p; }
466 :
467 : // Set or reset the channel selection on all iterators
468 0 : void VisSet::selectChannel(Int nGroup,Int start, Int width, Int increment,
469 : Int spectralWindow) {
470 :
471 : // Delegate, with callOrigin=true:
472 0 : VisSet::selectChannel(nGroup,start,width,increment,spectralWindow,true);
473 :
474 0 : }
475 :
476 79 : void VisSet::selectChannel(Int nGroup,Int start, Int width, Int increment,
477 : Int spectralWindow,
478 : Bool callOrigin)
479 : {
480 79 : iter_p->selectChannel(nGroup,start,width,increment,spectralWindow);
481 79 : if (callOrigin) iter_p->origin();
482 :
483 79 : selection_p(0,spectralWindow) = start;
484 79 : selection_p(1,spectralWindow) = width;
485 :
486 79 : }
487 :
488 0 : void VisSet::selectChannel(const Matrix<Int>& chansel) {
489 :
490 : // This doesn't handle multiple selections per spw!
491 :
492 0 : LogSink logSink;
493 0 : LogMessage message(LogOrigin("VisSet","selectChannel"));
494 :
495 0 : for (uInt ispw=0;ispw<chansel.nrow();++ispw) {
496 0 : const Int& spw=chansel(ispw,0);
497 0 : const Int& start=chansel(ispw,1);
498 0 : const Int& end=chansel(ispw,2);
499 0 : Int nchan=end-start+1;
500 0 : const Int& step=chansel(ispw,3);
501 :
502 0 : ostringstream os;
503 0 : os << ". Spw " << spw << ":" << start << "~" << end
504 0 : << " (" << nchan << " channels, step by " << step << ")";
505 0 : message.message(os);
506 0 : logSink.post(message);
507 :
508 0 : this->selectChannel(1,start,nchan,step,spw,false);
509 : }
510 :
511 0 : }
512 :
513 0 : void VisSet::selectAllChans() {
514 :
515 0 : LogSink logSink;
516 0 : LogMessage message(LogOrigin("VisSet","selectAllChans"));
517 0 : ostringstream os;
518 0 : os << ". Selecting all channels in selected spws.";
519 0 : message.message(os);
520 0 : logSink.post(message);
521 :
522 0 : MSSpWindowColumns msSpW(ms_p.spectralWindow());
523 0 : Int nSpw=msSpW.nrow();
524 0 : if(nSpw==0)
525 0 : throw(AipsError(String("There is no valid spectral windows in "+ms_p.tableName())));
526 0 : selection_p.resize(2,nSpw);
527 : // fill in default selection
528 0 : selection_p.row(0)=0; //start
529 0 : selection_p.row(1)=msSpW.numChan().getColumn();
530 :
531 : // Pass to the VisIter...
532 0 : for (uInt spw=0; spw<selection_p.ncolumn(); ++spw) {
533 0 : iter_p->selectChannel(1,selection_p(0,spw),selection_p(1,spw),0,spw);
534 : }
535 0 : }
536 :
537 163 : Int VisSet::numberAnt()
538 : {
539 163 : if(iter_p->newMS()){
540 126 : ms_p=(*blockOfMS_p)[iter_p->msId()];
541 : }
542 :
543 163 : return ((MeasurementSet&)ms_p).antenna().nrow(); // for single (sub)array only..
544 : }
545 39 : Int VisSet::numberFld()
546 : {
547 :
548 39 : if(iter_p->newMS()){
549 2 : ms_p=(*blockOfMS_p)[iter_p->msId()];
550 : }
551 39 : return ((MeasurementSet&)ms_p).field().nrow();
552 : }
553 4892 : Int VisSet::numberSpw()
554 : {
555 :
556 4892 : if(iter_p->newMS()){
557 4825 : ms_p=(*blockOfMS_p)[iter_p->msId()];
558 : }
559 4892 : return ((MeasurementSet&)ms_p).spectralWindow().nrow();
560 : }
561 65 : Vector<Int> VisSet::startChan() const
562 : {
563 :
564 65 : return selection_p.row(0);
565 : }
566 80 : Vector<Int> VisSet::numberChan() const
567 : {
568 80 : return selection_p.row(1);
569 : }
570 0 : Int VisSet::numberCoh() const
571 : {
572 :
573 0 : Int numcoh=0;
574 0 : for (uInt k=0; k < blockOfMS_p->nelements(); ++k){
575 0 : numcoh+=(*blockOfMS_p)[k].nrow();
576 : }
577 0 : return numcoh;
578 : }
579 :
580 0 : void VisSet::addCalSet(MeasurementSet& ms, Bool compress, Bool doModelData) {
581 :
582 : // Add and initialize calibration set (comprising a set of CORRECTED_DATA
583 : // and MODEL_DATA columns) to the MeasurementSet.
584 :
585 0 : LogSink logSink;
586 0 : LogMessage message(LogOrigin("VisSet","addCalSet"));
587 :
588 0 : message.message("Adding MODEL_DATA (initialized to unity) and CORRECTED_DATA (initialized to DATA) columns");
589 0 : logSink.post(message);
590 :
591 : {
592 : // Define a column accessor to the observed data
593 : TableColumn* data;
594 0 : const bool data_is_float = ms.tableDesc().isColumn(MS::columnName(MS::FLOAT_DATA));
595 0 : if (data_is_float) {
596 0 : data = new ArrayColumn<Float> (ms, MS::columnName(MS::FLOAT_DATA));
597 : } else {
598 0 : data = new ArrayColumn<Complex> (ms, MS::columnName(MS::DATA));
599 : };
600 :
601 : // Check if the data column is tiled and, if so, the
602 : // smallest tile shape used.
603 0 : TableDesc td = ms.actualTableDesc();
604 0 : const ColumnDesc& cdesc = td[data->columnDesc().name()];
605 0 : String dataManType = cdesc.dataManagerType();
606 0 : String dataManGroup = cdesc.dataManagerGroup();
607 :
608 0 : IPosition dataTileShape;
609 0 : Bool tiled = (dataManType.contains("Tiled"));
610 0 : Bool simpleTiling = false;
611 :
612 :
613 0 : if (tiled) {
614 0 : ROTiledStManAccessor tsm(ms, dataManGroup);
615 0 : uInt nHyper = tsm.nhypercubes();
616 : // Find smallest tile shape
617 0 : Int lowestProduct=INT_MAX,highestProduct=-INT_MAX;
618 0 : Int highestId=0;
619 0 : for (uInt id=0; id < nHyper; id++) {
620 0 : Int product = tsm.getTileShape(id).product();
621 0 : if (product > 0 && (product < lowestProduct)) {
622 0 : lowestProduct = product;
623 : //lowestId = id;
624 : };
625 0 : if (product > 0 && (product > highestProduct)) {
626 0 : highestProduct = product;
627 0 : highestId = id;
628 : };
629 : };
630 :
631 : // 2010Oct07 (gmoellen) We will try using
632 : // maximum volumne intput tile shape, as this
633 : // improves things drastically for ALMA data with
634 : // an enormous range of nchan (e.g., 3840+:1), and
635 : // will have no impact on data with a single shape
636 : // dataTileShape = tsm.getTileShape(lowestId);
637 0 : dataTileShape = tsm.getTileShape(highestId);
638 0 : simpleTiling = (dataTileShape.nelements() == 3);
639 :
640 : };
641 :
642 0 : if (!tiled || !simpleTiling) {
643 : // Untiled, or tiled at a higher than expected dimensionality
644 : // Use a canonical tile shape of 1 MB size
645 :
646 0 : Int maxNchan = max (numberChan());
647 0 : Int tileSize = maxNchan/10 + 1;
648 0 : Int nCorr = data->shape(0)(0);
649 0 : dataTileShape = IPosition(3, nCorr, tileSize, 131072/nCorr/tileSize + 1);
650 : };
651 :
652 0 : if(doModelData){
653 : // Add the MODEL_DATA column
654 0 : TableDesc tdModel, tdModelComp, tdModelScale;
655 0 : CompressComplex* ccModel=NULL;
656 0 : String colModel=MS::columnName(MS::MODEL_DATA);
657 :
658 0 : tdModel.addColumn(ArrayColumnDesc<Complex>(colModel,"model data", 2));
659 0 : td.addColumn(ArrayColumnDesc<Complex>(colModel,"model data", 2));
660 0 : IPosition modelTileShape = dataTileShape;
661 0 : if (compress) {
662 0 : tdModelComp.addColumn(ArrayColumnDesc<Int>(colModel+"_COMPRESSED",
663 0 : "model data compressed",2));
664 0 : tdModelScale.addColumn(ScalarColumnDesc<Float>(colModel+"_SCALE"));
665 0 : tdModelScale.addColumn(ScalarColumnDesc<Float>(colModel+"_OFFSET"));
666 0 : ccModel = new CompressComplex(colModel, colModel+"_COMPRESSED",
667 0 : colModel+"_SCALE", colModel+"_OFFSET", true);
668 :
669 0 : StandardStMan modelScaleStMan("ModelScaleOffset");
670 0 : ms.addColumn(tdModelScale, modelScaleStMan);
671 :
672 :
673 0 : TiledShapeStMan modelCompStMan("ModelCompTiled", modelTileShape);
674 0 : ms.addColumn(tdModelComp, modelCompStMan);
675 0 : ms.addColumn(tdModel, *ccModel);
676 :
677 : } else {
678 0 : MeasurementSet::addColumnToDesc(tdModel, MeasurementSet::MODEL_DATA, 2);
679 : //MeasurementSet::addColumnToDesc(td, MeasurementSet::MODEL_DATA, 2);
680 0 : TiledShapeStMan modelStMan("ModelTiled", modelTileShape);
681 : //String hcolName=String("Tiled")+String("MODEL_DATA");
682 : //tdModel.defineHypercolumn(hcolName,3,
683 : // stringToVector(colModel));
684 : //td.defineHypercolumn(hcolName,3,
685 : // stringToVector(colModel));
686 0 : ms.addColumn(tdModel, modelStMan);
687 : };
688 0 : if (ccModel) delete ccModel;
689 : }
690 :
691 :
692 : // Add the CORRECTED_DATA column
693 0 : TableDesc tdCorr, tdCorrComp, tdCorrScale;
694 0 : CompressComplex* ccCorr=NULL;
695 0 : String colCorr=MS::columnName(MS::CORRECTED_DATA);
696 :
697 0 : tdCorr.addColumn(ArrayColumnDesc<Complex>(colCorr,"corrected data", 2));
698 0 : IPosition corrTileShape = dataTileShape;
699 0 : if (compress) {
700 0 : tdCorrComp.addColumn(ArrayColumnDesc<Int>(colCorr+"_COMPRESSED",
701 0 : "corrected data compressed",2));
702 0 : tdCorrScale.addColumn(ScalarColumnDesc<Float>(colCorr+"_SCALE"));
703 0 : tdCorrScale.addColumn(ScalarColumnDesc<Float>(colCorr+"_OFFSET"));
704 0 : ccCorr = new CompressComplex(colCorr, colCorr+"_COMPRESSED",
705 0 : colCorr+"_SCALE", colCorr+"_OFFSET", true);
706 :
707 0 : StandardStMan corrScaleStMan("CorrScaleOffset");
708 0 : ms.addColumn(tdCorrScale, corrScaleStMan);
709 :
710 0 : TiledShapeStMan corrCompStMan("CorrectedCompTiled", corrTileShape);
711 0 : ms.addColumn(tdCorrComp, corrCompStMan);
712 0 : ms.addColumn(tdCorr, *ccCorr);
713 :
714 : } else {
715 0 : TiledShapeStMan corrStMan("CorrectedTiled", corrTileShape);
716 0 : ms.addColumn(tdCorr, corrStMan);
717 : };
718 : //MeasurementSet::addColumnToDesc(td, MeasurementSet::CORRECTED_DATA, 2);
719 :
720 :
721 0 : if (ccCorr) delete ccCorr;
722 :
723 : /* Set the shapes for each row
724 : and initialize CORRECTED_DATA to DATA
725 : and MODEL_DATA to one
726 : */
727 0 : ArrayColumn<Complex> modelData;
728 0 : if(doModelData)
729 0 : modelData.attach(ms, "MODEL_DATA");
730 0 : ArrayColumn<Complex> correctedData(ms, "CORRECTED_DATA");
731 :
732 : // Get data description column
733 0 : ScalarColumn<Int> dd_col = MSMainColumns(ms).dataDescId();
734 :
735 : // Get polarization column in dd table
736 0 : ScalarColumn<Int> pol_col = MSDataDescColumns(ms.dataDescription()).polarizationId();
737 :
738 : // Get correlation column
739 0 : ArrayColumn<Int> corr_col(MSColumns(ms).polarization().corrType());
740 :
741 0 : Matrix<Complex> model_vis;
742 0 : Int last_dd_id = 0;
743 : /* Initialize last_dd_id to something that
744 : causes the model_vis to be (re-)computed
745 : for the first MS row
746 : */
747 0 : if (ms.nrow() > 0) {
748 0 : last_dd_id = dd_col(0) + 1;
749 : }
750 :
751 0 : for (uInt row=0; row < ms.nrow(); row++) {
752 :
753 0 : IPosition rowShape=data->shape(row); // shape of (FLOAT_)DATA column
754 :
755 0 : correctedData.setShape(row, rowShape);
756 0 : if(doModelData)
757 0 : modelData.setShape(row, rowShape);
758 :
759 0 : Matrix<Complex> vis(rowShape);
760 :
761 0 : if (data_is_float) {
762 : /* Convert to complex for the CORRECTED_DATA column */
763 0 : Matrix<Float> f(rowShape);
764 0 : dynamic_cast<ArrayColumn<Float>*>(data)->get(row, f);
765 :
766 0 : for (unsigned i = 0; (int)i < f.shape()(0); i++)
767 0 : for (unsigned j = 0; (int)j < f.shape()(1); j++)
768 0 : vis(i, j) = f(i, j);
769 : }
770 : else {
771 0 : dynamic_cast<ArrayColumn<Complex>*>(data)->get(row, vis);
772 : }
773 :
774 0 : correctedData.put(row, vis);
775 :
776 : // figure out which correlations to set to 1. and 0. for the model.
777 : // Only do that, if the ddid changed since the last row
778 0 : Int dd_id = dd_col(row);
779 0 : if(doModelData){
780 0 : if (dd_id != last_dd_id) {
781 :
782 0 : last_dd_id = dd_id;
783 0 : model_vis.resize(rowShape);
784 :
785 0 : Int pol_id = pol_col(dd_id);
786 0 : Vector<Int> corrType = corr_col(pol_id);
787 :
788 : //cerr << "row = " << row << ", dd = " << dd_id << ", pol = " << pol_id << ", corr = " << corrType << endl;
789 :
790 0 : Vector<Int> lastCorrType;
791 0 : Vector<Bool> zero;
792 :
793 0 : uInt nCorr = corrType.nelements();
794 0 : if (nCorr != lastCorrType.nelements() ||
795 0 : !allEQ(corrType, lastCorrType)) {
796 :
797 0 : lastCorrType.resize(nCorr);
798 0 : lastCorrType = corrType;
799 0 : zero.resize(nCorr);
800 :
801 0 : for (uInt i=0; i<nCorr; i++)
802 : {
803 0 : zero[i]=(corrType[i]==Stokes::RL || corrType[i]==Stokes::LR ||
804 0 : corrType[i]==Stokes::XY || corrType[i]==Stokes::YX);
805 : }
806 : }
807 :
808 0 : model_vis = Complex(1.0,0.0);
809 :
810 0 : for (uInt i=0; i < nCorr; i++) {
811 0 : if (zero[i]) {
812 0 : model_vis(Slice(i), Slice()) = Complex(0.0, 0.0);
813 : }
814 : }
815 : } // endif ddid changed
816 0 : modelData.put(row, model_vis);
817 : }
818 : }
819 :
820 0 : delete data;
821 : }
822 :
823 0 : ostringstream os;
824 0 : os << "Initialized " << ms.nrow() << " rows";
825 0 : message.message(os.str());
826 0 : logSink.post(message);
827 :
828 0 : return;
829 : }
830 :
831 0 : void VisSet::addCalSet2(MeasurementSet& ms, Bool compress, Bool doModelData) {
832 :
833 : // Add but DO NOT INITIALIZE calibration set (comprising a set of CORRECTED_DATA
834 : // and MODEL_DATA columns) to the MeasurementSet.
835 :
836 0 : LogSink logSink;
837 0 : LogMessage message(LogOrigin("VisSet","addCalSet"));
838 0 : doModelData ?
839 0 : message.message("Adding MODEL_DATA and CORRECTED_DATA columns"):
840 0 : message.message("Adding CORRECTED_DATA columns");
841 0 : logSink.post(message);
842 :
843 : {
844 : // Define a column accessor to the observed data
845 : TableColumn* data;
846 0 : const bool data_is_float = ms.tableDesc().isColumn(MS::columnName(MS::FLOAT_DATA));
847 0 : if (data_is_float) {
848 0 : data = new ArrayColumn<Float> (ms, MS::columnName(MS::FLOAT_DATA));
849 : } else {
850 0 : data = new ArrayColumn<Complex> (ms, MS::columnName(MS::DATA));
851 : };
852 :
853 : // Check if the data column is tiled and, if so, the
854 : // smallest tile shape used.
855 0 : TableDesc td = ms.actualTableDesc();
856 0 : const ColumnDesc& cdesc = td[data->columnDesc().name()];
857 0 : String dataManType = cdesc.dataManagerType();
858 0 : String dataManGroup = cdesc.dataManagerGroup();
859 :
860 0 : IPosition dataTileShape;
861 0 : Bool tiled = (dataManType.contains("Tiled"));
862 0 : Bool simpleTiling = false;
863 :
864 :
865 0 : if (tiled) {
866 0 : ROTiledStManAccessor tsm(ms, dataManGroup);
867 0 : uInt nHyper = tsm.nhypercubes();
868 : // Find smallest tile shape
869 0 : Int lowestProduct=INT_MAX,highestProduct=-INT_MAX;
870 0 : Int highestId=0;
871 0 : for (uInt id=0; id < nHyper; id++) {
872 0 : Int product = tsm.getTileShape(id).product();
873 0 : if (product > 0 && (product < lowestProduct)) {
874 0 : lowestProduct = product;
875 : //lowestId = id;
876 : };
877 0 : if (product > 0 && (product > highestProduct)) {
878 0 : highestProduct = product;
879 0 : highestId = id;
880 : };
881 : };
882 :
883 : // 2010Oct07 (gmoellen) We will try using
884 : // maximum volumne intput tile shape, as this
885 : // improves things drastically for ALMA data with
886 : // an enormous range of nchan (e.g., 3840+:1), and
887 : // will have no impact on data with a single shape
888 : // dataTileShape = tsm.getTileShape(lowestId);
889 0 : dataTileShape = tsm.getTileShape(highestId);
890 0 : simpleTiling = (dataTileShape.nelements() == 3);
891 :
892 : };
893 :
894 0 : if (!tiled || !simpleTiling) {
895 : // Untiled, or tiled at a higher than expected dimensionality
896 : // Use a canonical tile shape of 1 MB size
897 :
898 0 : Int maxNchan = max (numberChan());
899 0 : Int tileSize = maxNchan/10 + 1;
900 0 : Int nCorr = data->shape(0)(0);
901 0 : dataTileShape = IPosition(3, nCorr, tileSize, 131072/nCorr/tileSize + 1);
902 : };
903 :
904 0 : delete data;
905 :
906 0 : if(doModelData){
907 : // Add the MODEL_DATA column
908 0 : TableDesc tdModel, tdModelComp, tdModelScale;
909 0 : CompressComplex* ccModel=NULL;
910 0 : String colModel=MS::columnName(MS::MODEL_DATA);
911 :
912 0 : tdModel.addColumn(ArrayColumnDesc<Complex>(colModel,"model data", 2));
913 0 : td.addColumn(ArrayColumnDesc<Complex>(colModel,"model data", 2));
914 0 : IPosition modelTileShape = dataTileShape;
915 0 : if (compress) {
916 0 : tdModelComp.addColumn(ArrayColumnDesc<Int>(colModel+"_COMPRESSED",
917 0 : "model data compressed",2));
918 0 : tdModelScale.addColumn(ScalarColumnDesc<Float>(colModel+"_SCALE"));
919 0 : tdModelScale.addColumn(ScalarColumnDesc<Float>(colModel+"_OFFSET"));
920 0 : ccModel = new CompressComplex(colModel, colModel+"_COMPRESSED",
921 0 : colModel+"_SCALE", colModel+"_OFFSET", true);
922 :
923 0 : StandardStMan modelScaleStMan("ModelScaleOffset");
924 0 : ms.addColumn(tdModelScale, modelScaleStMan);
925 :
926 :
927 0 : TiledShapeStMan modelCompStMan("ModelCompTiled", modelTileShape);
928 0 : ms.addColumn(tdModelComp, modelCompStMan);
929 0 : ms.addColumn(tdModel, *ccModel);
930 :
931 : } else {
932 0 : MeasurementSet::addColumnToDesc(tdModel, MeasurementSet::MODEL_DATA, 2);
933 : //MeasurementSet::addColumnToDesc(td, MeasurementSet::MODEL_DATA, 2);
934 0 : TiledShapeStMan modelStMan("ModelTiled", modelTileShape);
935 : //String hcolName=String("Tiled")+String("MODEL_DATA");
936 : //tdModel.defineHypercolumn(hcolName,3,
937 : // stringToVector(colModel));
938 : //td.defineHypercolumn(hcolName,3,
939 : // stringToVector(colModel));
940 0 : ms.addColumn(tdModel, modelStMan);
941 : };
942 0 : if (ccModel) delete ccModel;
943 : }
944 : // Add the CORRECTED_DATA column
945 0 : TableDesc tdCorr, tdCorrComp, tdCorrScale;
946 0 : CompressComplex* ccCorr=NULL;
947 0 : String colCorr=MS::columnName(MS::CORRECTED_DATA);
948 :
949 0 : tdCorr.addColumn(ArrayColumnDesc<Complex>(colCorr,"corrected data", 2));
950 0 : IPosition corrTileShape = dataTileShape;
951 0 : if (compress) {
952 0 : tdCorrComp.addColumn(ArrayColumnDesc<Int>(colCorr+"_COMPRESSED",
953 0 : "corrected data compressed",2));
954 0 : tdCorrScale.addColumn(ScalarColumnDesc<Float>(colCorr+"_SCALE"));
955 0 : tdCorrScale.addColumn(ScalarColumnDesc<Float>(colCorr+"_OFFSET"));
956 0 : ccCorr = new CompressComplex(colCorr, colCorr+"_COMPRESSED",
957 0 : colCorr+"_SCALE", colCorr+"_OFFSET", true);
958 :
959 0 : StandardStMan corrScaleStMan("CorrScaleOffset");
960 0 : ms.addColumn(tdCorrScale, corrScaleStMan);
961 :
962 0 : TiledShapeStMan corrCompStMan("CorrectedCompTiled", corrTileShape);
963 0 : ms.addColumn(tdCorrComp, corrCompStMan);
964 0 : ms.addColumn(tdCorr, *ccCorr);
965 :
966 : } else {
967 0 : TiledShapeStMan corrStMan("CorrectedTiled", corrTileShape);
968 0 : ms.addColumn(tdCorr, corrStMan);
969 : };
970 : //MeasurementSet::addColumnToDesc(td, MeasurementSet::CORRECTED_DATA, 2);
971 :
972 :
973 0 : if (ccCorr) delete ccCorr;
974 :
975 0 : ms.flush();
976 :
977 : }
978 :
979 0 : return;
980 : }
981 :
982 :
983 :
984 :
985 0 : void VisSet::addScratchCols(MeasurementSet& ms, Bool compress, Bool doModelData){
986 :
987 0 : LogSink logSink;
988 0 : LogMessage message(LogOrigin("VisSet","VisSet"));
989 :
990 : //function to add scratchy column
991 0 : Bool init=true;
992 :
993 0 : Int nSpw=ms.spectralWindow().nrow();
994 0 : MSSpWindowColumns msSpW(ms.spectralWindow());
995 0 : selection_p.resize(2,nSpw);
996 : // fill in default selection
997 0 : selection_p.row(0)=0; //start
998 0 : selection_p.row(1)=msSpW.numChan().getColumn();
999 :
1000 : // Add scratch columns
1001 0 : if (init) {
1002 :
1003 0 : VisSetUtil::removeCalSet(ms);
1004 0 : addCalSet(ms, compress, doModelData);
1005 :
1006 :
1007 : // Force re-sort if it was sorted
1008 0 : if (ms.keywordSet().isDefined("SORT_COLUMNS"))
1009 0 : ms.rwKeywordSet().removeField("SORT_COLUMNS");
1010 0 : if (ms.keywordSet().isDefined("SORTED_TABLE"))
1011 0 : ms.rwKeywordSet().removeField("SORTED_TABLE");
1012 : }
1013 0 : }
1014 :
1015 248 : String VisSet::msName(){
1016 :
1017 248 : String name;
1018 248 : if(ms_p.isMarkedForDelete()){ // this is a temporary selected table
1019 260 : Block<String> refTables = ms_p.getPartNames(false);
1020 130 : name = refTables[0];
1021 : }
1022 : else{
1023 118 : name = ms_p.tableName();
1024 : }
1025 248 : return name;
1026 :
1027 : }
1028 :
1029 26 : String VisSet::sysCalTableName(){
1030 :
1031 26 : return ms_p.sysCalTableName();
1032 : }
1033 :
1034 0 : String VisSet::spectralWindowTableName(){
1035 :
1036 0 : return ms_p.spectralWindowTableName();
1037 : }
1038 :
1039 0 : String VisSet::fieldTableName(){
1040 :
1041 0 : return ms_p.fieldTableName();
1042 : }
1043 :
1044 0 : String VisSet::syspowerTableName(){
1045 :
1046 0 : return ms_p.rwKeywordSet().asTable("SYSPOWER").tableName();
1047 : }
1048 :
1049 0 : String VisSet::caldeviceTableName(){
1050 :
1051 0 : return ms_p.rwKeywordSet().asTable("CALDEVICE").tableName();
1052 : }
1053 :
1054 :
1055 :
1056 :
1057 : } //# NAMESPACE CASA - END
1058 :
|
