Line data Source code
1 : //# VisibilityIterator.cc: Step through MeasurementEquation by visibility
2 : //# Copyright (C) 1996-2012
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: VisibilityIterator.cc,v 19.15 2006/02/01 01:25:14 kgolap Exp $
27 :
28 : #include <stdcasa/UtilJ.h>
29 : #include <msvis/MSVis/VisibilityIteratorImpl.h>
30 : #include <msvis/MSVis/VisibilityIterator.h>
31 : #include <msvis/MSVis/VisBuffer.h>
32 : #include <msvis/MSVis/MSUtil.h>
33 : ////#include <synthesis/TransformMachines/VisModelData.h>
34 : #include <casacore/scimath/Mathematics/InterpolateArray1D.h>
35 : #include <casacore/ms/MeasurementSets/MSColumns.h>
36 : #include <casacore/ms/MSSel/MSSpwIndex.h>
37 : #include <casacore/tables/Tables/TableDesc.h>
38 : #include <casacore/tables/Tables/ColDescSet.h>
39 : #include <casacore/tables/Tables/TableRecord.h>
40 : #include <casacore/tables/DataMan/TiledStManAccessor.h>
41 : #include <casacore/tables/DataMan/StandardStManAccessor.h>
42 : #include <casacore/tables/DataMan/IncrStManAccessor.h>
43 : #include <casacore/casa/Arrays/ArrayLogical.h>
44 : #include <casacore/casa/System/AipsrcValue.h>
45 : #include <casacore/casa/BasicSL/Constants.h>
46 : #include <casacore/casa/Quanta/MVTime.h>
47 : #include <casacore/casa/Containers/Record.h>
48 : #include <casacore/casa/Arrays/ArrayMath.h>
49 : #include <casacore/casa/Arrays/MaskedArray.h>
50 : #include <casacore/casa/Exceptions/Error.h>
51 : #include <casacore/casa/Utilities/Assert.h>
52 : #include <casacore/casa/Utilities/Sort.h>
53 : #include <casacore/casa/OS/EnvVar.h>
54 :
55 : #include <cassert>
56 : #include <limits>
57 : #include <memory>
58 :
59 : using std::make_pair;
60 :
61 : using namespace casacore;
62 : namespace casa { //# NAMESPACE CASA - BEGIN
63 :
64 : SubChunkPair
65 0 : SubChunkPair::noMoreData ()
66 : {
67 0 : Int maxInt = std::numeric_limits<Int>::max ();
68 0 : return SubChunkPair (maxInt, maxInt);
69 : }
70 :
71 : String
72 0 : SubChunkPair::toString () const
73 : {
74 0 : return String::format ("(%d,%d)", first, second);
75 : }
76 :
77 0 : VisibilityIteratorReadImpl::VisibilityIteratorReadImpl ()
78 : : rovi_p (NULL),
79 : selRows_p (0, 0),
80 0 : tileCacheIsSet_p (0)
81 0 : {}
82 :
83 2 : VisibilityIteratorReadImpl::VisibilityIteratorReadImpl (ROVisibilityIterator * rovi,
84 : const Block<MeasurementSet> &mss,
85 : const Block<Int> & sortColumns,
86 : const Bool addDefaultSort,
87 2 : Double timeInterval)
88 : : addDefaultSort_p (addDefaultSort),
89 : curChanGroup_p (0),
90 : floatDataFound_p (false),
91 : initialized_p (false),
92 : msIterAtOrigin_p (false),
93 : msIter_p (mss, sortColumns, timeInterval, addDefaultSort, false),
94 : nChan_p (0),
95 : nRowBlocking_p (0),
96 : rovi_p (NULL),
97 : selRows_p (0, 0),
98 : sortColumns_p (sortColumns),
99 : stateOk_p (false),
100 : tileCacheIsSet_p (0),
101 2 : timeInterval_p (timeInterval)
102 : {
103 : // Make sure the pointer to the containing ROVI (rovi_p) is NULL when calling initialize
104 : // otherwise the call back to the VI can result in it trying to use an uninitialized pointer
105 : // to this object (since it is in the process or being constructed).
106 :
107 2 : initialize (mss);
108 :
109 2 : rovi_p = rovi;
110 :
111 2 : }
112 :
113 : void
114 2 : VisibilityIteratorReadImpl::initialize (const Block<MeasurementSet> &mss)
115 : {
116 :
117 2 : AipsrcValue<Bool>::find (autoTileCacheSizing_p, VisibilityIterator::getAipsRcBase () + ".AutoTileCacheSizing", false);
118 :
119 2 : asyncEnabled_p = false;
120 2 : cache_p.lastazelUT_p = -1;
121 2 : cache_p.lastazel0UT_p = -1;
122 2 : cache_p.lasthourangUT_p = -1;
123 2 : cache_p.lastfeedpaUT_p = -1;
124 2 : cache_p.lastParangUT_p = -1;
125 2 : cache_p.lastParang0UT_p = -1;
126 :
127 2 : msCounter_p = 0;
128 :
129 2 : Int numMS = mss.nelements ();
130 2 : isMultiMS_p = numMS > 1;
131 :
132 4 : Block<Vector<Int> > blockNGroup (numMS);
133 4 : Block<Vector<Int> > blockStart (numMS);
134 4 : Block<Vector<Int> > blockWidth (numMS);
135 4 : Block<Vector<Int> > blockIncr (numMS);
136 4 : Block<Vector<Int> > blockSpw (numMS);
137 :
138 2 : measurementSets_p.clear ();
139 :
140 4 : for (Int k = 0; k < numMS; ++k) {
141 :
142 4 : MSSpWindowColumns msSpW (mss[k].spectralWindow ());
143 :
144 2 : Int nspw = msSpW.nrow ();
145 :
146 2 : blockNGroup[k].resize (nspw);
147 2 : blockNGroup[k].set (1);
148 2 : blockStart[k].resize (nspw);
149 2 : blockStart[k].set (0);
150 2 : blockWidth[k].resize (nspw);
151 2 : blockWidth[k] = msSpW.numChan ().getColumn ();
152 2 : blockIncr[k].resize (nspw);
153 2 : blockIncr[k].set (1);
154 2 : blockSpw[k].resize (nspw);
155 2 : indgen (blockSpw[k]);
156 :
157 2 : measurementSets_p.push_back (mss [k]);
158 : }
159 :
160 : selectChannel (blockNGroup, blockStart, blockWidth, blockIncr,
161 2 : blockSpw);
162 2 : }
163 :
164 :
165 0 : VisibilityIteratorReadImpl::VisibilityIteratorReadImpl (const VisibilityIteratorReadImpl & other,
166 0 : ROVisibilityIterator * rovi)
167 : : rovi_p (rovi),
168 0 : selRows_p (other.selRows_p) // no default constructor so init it here
169 : {
170 0 : operator=(other);
171 0 : }
172 :
173 2 : VisibilityIteratorReadImpl::~VisibilityIteratorReadImpl ()
174 : {
175 2 : }
176 :
177 : VisibilityIteratorReadImpl &
178 0 : VisibilityIteratorReadImpl::operator=(const VisibilityIteratorReadImpl & other)
179 : {
180 0 : if (this == &other) {
181 0 : return *this;
182 : }
183 :
184 0 : cache_p = other.cache_p;
185 0 : channels_p = other.channels_p;
186 0 : columns_p = other.columns_p;
187 0 : msChannels_p = other.msChannels_p;
188 0 : velocity_p = other.velocity_p;
189 :
190 0 : addDefaultSort_p = other.addDefaultSort_p;
191 0 : channelGroupSize_p = other.channelGroupSize_p;
192 0 : curChanGroup_p = other.curChanGroup_p;
193 0 : curEndRow_p = other.curEndRow_p;
194 0 : curChanGroup_p = other.curChanGroup_p;
195 0 : curNumRow_p = other.curNumRow_p;
196 0 : curStartRow_p = other.curStartRow_p;
197 0 : curTableNumRow_p = other.curTableNumRow_p;
198 0 : floatDataFound_p = other.floatDataFound_p;
199 0 : imwgt_p = other.imwgt_p;
200 0 : initialized_p = other.initialized_p;
201 0 : isMultiMS_p = other.isMultiMS_p;
202 0 : measurementSets_p = other.measurementSets_p;
203 0 : more_p = other.more_p;
204 0 : msCounter_p = other.msCounter_p;
205 0 : msIterAtOrigin_p = other.msIterAtOrigin_p;
206 0 : msIter_p = other.msIter_p;
207 0 : msIter_p.origin();
208 0 : msd_p = other.msd_p;
209 0 : nAnt_p = other.nAnt_p;
210 0 : nChan_p = other.nChan_p;
211 0 : nPol_p = other.nPol_p;
212 0 : nRowBlocking_p = other.nRowBlocking_p;
213 0 : newChanGroup_p = other.newChanGroup_p;
214 0 : selRows_p = other.selRows_p;
215 0 : slicer_p = other.slicer_p;
216 0 : sortColumns_p = other.sortColumns_p;
217 0 : stateOk_p = other.stateOk_p;
218 0 : tileCacheIsSet_p.resize ();
219 0 : timeInterval_p = other.timeInterval_p;
220 0 : time_p.assign (other.time_p);
221 0 : useSlicer_p = other.useSlicer_p;
222 0 : weightSlicer_p = other.weightSlicer_p;
223 :
224 0 : return *this;
225 : }
226 :
227 : VisibilityIteratorReadImpl::Velocity &
228 0 : VisibilityIteratorReadImpl::Velocity::operator= (const VisibilityIteratorReadImpl::Velocity & other)
229 : {
230 0 : cFromBETA_p = other.cFromBETA_p;
231 0 : lsrFreq_p.assign (other.lsrFreq_p);
232 0 : nVelChan_p = other.nVelChan_p;
233 0 : selFreq_p = other.selFreq_p;
234 0 : vDef_p = other.vDef_p;
235 0 : vInc_p = other.vInc_p;
236 0 : vInterpolation_p = other.vInterpolation_p;
237 0 : vPrecise_p = other.vPrecise_p;
238 0 : vStart_p = other.vStart_p;
239 0 : velSelection_p = other.velSelection_p;
240 :
241 0 : return * this;
242 : }
243 :
244 2 : VisibilityIteratorReadImpl::Cache::Cache()
245 : : flagOK_p (false),
246 : floatDataCubeOK_p (false),
247 : freqCacheOK_p (false),
248 : hourang_p (0),
249 : lastParang0UT_p (-1), // set last cache update as invalid
250 : lastParangUT_p (-1), // set last cache update as invalid
251 : lastazelUT_p (-1), // set last cache update as invalid
252 : lastazel0UT_p (-1), // set last cache update as invalid
253 : lasthourangUT_p (-1), // set last cache update as invalid
254 : lastfeedpaUT_p (-1), // set last cache update as invalid
255 : msHasFC_p(false),
256 : msHasWtSp_p (false),
257 : parang0_p (0),
258 2 : weightSpOK_p (false)
259 2 : {}
260 :
261 : VisibilityIteratorReadImpl::Cache &
262 0 : VisibilityIteratorReadImpl::Cache::operator= (const VisibilityIteratorReadImpl::Cache & other)
263 : {
264 0 : azel0_p = other.azel0_p;
265 0 : azel_p.assign (other.azel_p);
266 0 : feedpa_p.assign (other.feedpa_p);
267 0 : flagCube_p.assign (other.flagCube_p);
268 0 : flagOK_p = other.flagOK_p;
269 0 : floatDataCubeOK_p = other.floatDataCubeOK_p;
270 0 : floatDataCube_p.assign (other.floatDataCube_p);
271 0 : freqCacheOK_p = other.freqCacheOK_p;
272 0 : frequency_p.assign (frequency_p);
273 0 : hourang_p = other.hourang_p;
274 0 : lastazelUT_p = other.lastazelUT_p;
275 0 : lastazel0UT_p = other.lastazel0UT_p;
276 0 : lasthourangUT_p = other.lasthourangUT_p;
277 0 : lastfeedpaUT_p = other.lastfeedpaUT_p;
278 0 : lastParangUT_p = other.lastParangUT_p;
279 0 : lastParang0UT_p = other.lastParang0UT_p;
280 0 : msHasFC_p = other.msHasFC_p;
281 0 : msHasWtSp_p = other.msHasWtSp_p;
282 0 : parang0_p = other.parang0_p;
283 0 : parang_p.assign (other.parang_p);
284 0 : rowIds_p = other.rowIds_p;
285 0 : uvwMat_p.assign (other.uvwMat_p);
286 0 : visCube_p.assign (other.visCube_p);
287 0 : visOK_p = other.visOK_p;
288 0 : weightSpOK_p = other.weightSpOK_p;
289 0 : wtSp_p.assign (other.wtSp_p);
290 :
291 0 : return * this;
292 : }
293 :
294 : VisibilityIteratorReadImpl::Columns &
295 0 : VisibilityIteratorReadImpl::Columns::operator= (const VisibilityIteratorReadImpl::Columns & other)
296 : {
297 0 : antenna1_p.reference (other.antenna1_p);
298 0 : antenna2_p.reference (other.antenna2_p);
299 0 : corrVis_p.reference (other.corrVis_p);
300 0 : exposure_p.reference (other.exposure_p);
301 0 : feed1_p.reference (other.feed1_p);
302 0 : feed2_p.reference (other.feed2_p);
303 0 : flag_p.reference (other.flag_p);
304 0 : flagCategory_p.reference (other.flagCategory_p);
305 0 : flagRow_p.reference (other.flagRow_p);
306 0 : floatVis_p.reference (other.floatVis_p);
307 0 : modelVis_p.reference (other.modelVis_p);
308 0 : observation_p.reference (other.observation_p);
309 0 : processor_p.reference (other.processor_p);
310 0 : scan_p.reference (other.scan_p);
311 0 : sigma_p.reference (other.sigma_p);
312 0 : state_p.reference (other.state_p);
313 0 : time_p.reference (other.time_p);
314 0 : timeCentroid_p.reference (other.timeCentroid_p);
315 0 : timeInterval_p.reference (other.timeInterval_p);
316 0 : uvw_p.reference (other.uvw_p);
317 0 : vis_p.reference (other.vis_p);
318 0 : weight_p.reference (other.weight_p);
319 0 : weightSpectrum_p.reference (other.weightSpectrum_p);
320 :
321 0 : return * this;
322 : }
323 :
324 2 : VisibilityIteratorReadImpl::Velocity::Velocity ()
325 : : nVelChan_p (0),
326 : velSelection_p (false),
327 2 : vPrecise_p (false)
328 2 : {}
329 :
330 :
331 : VisibilityIteratorReadImpl *
332 0 : VisibilityIteratorReadImpl::clone (ROVisibilityIterator * rovi) const
333 : {
334 0 : return new VisibilityIteratorReadImpl (* this, rovi);
335 : }
336 :
337 : void
338 0 : VisibilityIteratorReadImpl::setAsyncEnabled (Bool enabled)
339 : {
340 0 : asyncEnabled_p = enabled;
341 0 : }
342 :
343 : Bool
344 0 : VisibilityIteratorReadImpl::isAsyncEnabled () const
345 : {
346 0 : return asyncEnabled_p;
347 : }
348 :
349 :
350 : void
351 2 : VisibilityIteratorReadImpl::setRowBlocking (Int nRow)
352 : {
353 2 : nRowBlocking_p = nRow;
354 2 : }
355 :
356 : Bool
357 0 : VisibilityIteratorReadImpl::existsColumn (VisBufferComponents::EnumType id) const
358 : {
359 : Bool result;
360 0 : switch (id){
361 :
362 0 : case VisBufferComponents::Corrected:
363 : case VisBufferComponents::CorrectedCube:
364 :
365 0 : result = ! columns_p.corrVis_p.isNull();
366 0 : break;
367 :
368 0 : case VisBufferComponents::Model:
369 : case VisBufferComponents::ModelCube:
370 :
371 0 : result = ! columns_p.modelVis_p.isNull();
372 0 : break;
373 :
374 0 : case VisBufferComponents::Observed:
375 : case VisBufferComponents::ObservedCube:
376 :
377 0 : result = ! (columns_p.vis_p.isNull() && columns_p.floatVis_p.isNull());
378 0 : break;
379 :
380 : // RR: I can't tell if the other columns should checked for here or not.
381 : // It's not true that all the other columns are required.
382 : // existsFlagCategory uses caching anyway.
383 : // case VisBufferComponents::FlagCategory:
384 : // result = false;
385 : // if(!columns_p.flagCategory().isNull() &&
386 : // columns_p.flagCategory().isDefined(0)){
387 : // IPosition fcshape(columns_p.flagCategory().shape(0));
388 : // IPosition fshape(columns_p.flag().shape(0));
389 :
390 : // result = fcshape(0) == fshape(0) && fcshape(1) == fshape(1);
391 : // }
392 :
393 0 : default:
394 0 : result = true; // required columns
395 : }
396 :
397 0 : return result;
398 : }
399 :
400 :
401 : void
402 1 : VisibilityIteratorReadImpl::useImagingWeight (const VisImagingWeight & imWgt)
403 : {
404 1 : imwgt_p = imWgt;
405 1 : }
406 : void
407 22 : VisibilityIteratorReadImpl::origin ()
408 : {
409 :
410 22 : if (!initialized_p) {
411 0 : originChunks ();
412 : } else {
413 22 : curChanGroup_p = 0;
414 22 : newChanGroup_p = true;
415 22 : curStartRow_p = 0;
416 22 : cache_p.freqCacheOK_p = false;
417 22 : cache_p.flagOK_p = false;
418 22 : cache_p.weightSpOK_p = false;
419 22 : cache_p.visOK_p.resize (3);
420 22 : cache_p.visOK_p = false;
421 22 : cache_p.floatDataCubeOK_p = false;
422 22 : setSelTable ();
423 22 : attachColumnsSafe (attachTable ());
424 22 : getTopoFreqs ();
425 22 : updateSlicer ();
426 22 : more_p = curChanGroup_p < curNGroups_p;
427 : // invalidate any attached VisBuffer
428 22 : if (!vbStack_p.empty ()) {
429 11 : vbStack_p.top ()->invalidate ();
430 : }
431 :
432 22 : subchunk_p.resetSubChunk ();
433 : }
434 22 : }
435 :
436 : void
437 8 : VisibilityIteratorReadImpl::originChunks ()
438 : {
439 8 : originChunks (false);
440 8 : }
441 :
442 :
443 : void
444 8 : VisibilityIteratorReadImpl::originChunks (Bool forceRewind)
445 : {
446 8 : initialized_p = true;
447 8 : subchunk_p.resetToOrigin();
448 :
449 8 : if (forceRewind) {
450 0 : msIterAtOrigin_p = false;
451 : }
452 :
453 8 : if (!msIterAtOrigin_p) {
454 :
455 1 : msIter_p.origin ();
456 1 : msIterAtOrigin_p = true;
457 :
458 1 : while ((! isInSelectedSPW (msIter_p.spectralWindowId ())) &&
459 0 : msIter_p.more ()) {
460 :
461 0 : msIter_p++;
462 : }
463 :
464 1 : stateOk_p = false;
465 1 : msCounter_p = msId ();
466 :
467 : }
468 :
469 8 : setState ();
470 8 : origin ();
471 8 : setTileCache ();
472 8 : }
473 :
474 : Bool
475 16 : VisibilityIteratorReadImpl::isInSelectedSPW (const Int & spw)
476 : {
477 :
478 16 : for (uInt k = 0; k < msChannels_p.spw_p[msId ()].nelements () ; ++k) {
479 16 : if (spw == msChannels_p.spw_p[msId ()][k]) {
480 16 : return true;
481 : }
482 : }
483 0 : return false;
484 : }
485 :
486 : void
487 266 : VisibilityIteratorReadImpl::advance ()
488 : {
489 266 : newChanGroup_p = false;
490 266 : cache_p.flagOK_p = false;
491 266 : cache_p.visOK_p = false;
492 266 : cache_p.floatDataCubeOK_p = false;
493 266 : cache_p.weightSpOK_p = false;
494 266 : curStartRow_p = curEndRow_p + 1;
495 266 : if (curStartRow_p >= curTableNumRow_p) {
496 13 : if (++ curChanGroup_p >= curNGroups_p) {
497 13 : curChanGroup_p--;
498 13 : more_p = false;
499 : } else {
500 0 : curStartRow_p = 0;
501 0 : newChanGroup_p = true;
502 0 : cache_p.freqCacheOK_p = false;
503 0 : updateSlicer ();
504 : }
505 : }
506 266 : if (more_p) {
507 253 : subchunk_p.incrementSubChunk();
508 253 : setSelTable ();
509 253 : getTopoFreqs ();
510 : // invalidate any attached VisBuffer
511 253 : if (!vbStack_p.empty ()) {
512 253 : vbStack_p.top ()->invalidate ();
513 : }
514 : }
515 266 : }
516 :
517 : SubChunkPair
518 0 : VisibilityIteratorReadImpl::getSubchunkId () const
519 : {
520 0 : return subchunk_p;
521 : }
522 :
523 0 : const Block<Int>& VisibilityIteratorReadImpl::getSortColumns() const
524 : {
525 0 : return sortColumns_p;
526 : }
527 :
528 : VisibilityIteratorReadImpl &
529 13 : VisibilityIteratorReadImpl::nextChunk ()
530 : {
531 :
532 13 : if (msIter_p.more ()) {
533 13 : msIter_p++;
534 13 : if ((!isInSelectedSPW (msIter_p.spectralWindowId ()))) {
535 0 : while ( (!isInSelectedSPW (msIter_p.spectralWindowId ()))
536 0 : && (msIter_p.more ())) {
537 0 : msIter_p++;
538 : }
539 0 : stateOk_p = false;
540 : }
541 :
542 13 : if (msIter_p.newMS ()) {
543 0 : msCounter_p = msId ();
544 0 : doChannelSelection ();
545 : }
546 13 : msIterAtOrigin_p = false;
547 13 : stateOk_p = false;
548 : }
549 13 : if (msIter_p.more ()) {
550 10 : subchunk_p.incrementChunk();
551 10 : setState ();
552 10 : getTopoFreqs ();
553 10 : if (!vbStack_p.empty ()) {
554 5 : vbStack_p.top ()->invalidate ();
555 : }
556 : }
557 13 : more_p = msIter_p.more ();
558 13 : return *this;
559 : }
560 :
561 : void
562 288 : VisibilityIteratorReadImpl::setSelTable ()
563 : {
564 : // work out how many rows to return
565 : // for the moment we return all rows with the same value for time
566 : // unless row blocking is set, in which case we return more rows at once.
567 288 : if (nRowBlocking_p > 0) {
568 14 : curEndRow_p = curStartRow_p + nRowBlocking_p;
569 14 : if (curEndRow_p >= curTableNumRow_p) {
570 14 : curEndRow_p = curTableNumRow_p - 1;
571 : }
572 : } else {
573 195797 : for (curEndRow_p = curStartRow_p + 1; curEndRow_p < curTableNumRow_p &&
574 97895 : time_p (curEndRow_p) == time_p (curEndRow_p - 1);
575 97628 : curEndRow_p++) {
576 : ;
577 : }
578 274 : curEndRow_p--;
579 : }
580 :
581 288 : curNumRow_p = curEndRow_p - curStartRow_p + 1;
582 288 : selRows_p = RefRows (curStartRow_p, curEndRow_p);
583 288 : cache_p.rowIds_p.resize (0);
584 288 : }
585 :
586 : void
587 285 : VisibilityIteratorReadImpl::getTopoFreqs ()
588 : {
589 285 : if (velocity_p.velSelection_p) {
590 :
591 : // Convert selected velocities to TOPO frequencies.
592 : // First convert observatory vel to correct frame (for this time).
593 :
594 0 : msd_p.setEpoch (msIter_p.msColumns ().timeMeas ()(curStartRow_p));
595 0 : if (msIter_p.newMS ()) {
596 0 : msd_p.setObservatoryPosition (msIter_p.telescopePosition ());
597 : }
598 :
599 0 : MRadialVelocity obsRV = msd_p.obsVel (); // get obs velocity in required frame
600 :
601 0 : Double obsVel = velocity_p.cFromBETA_p (obsRV.toDoppler ()).getValue ().get ().getValue ();
602 : // convert to doppler in required definition and get out in m/s
603 :
604 : // Now compute corresponding TOPO freqs
605 :
606 0 : velocity_p.selFreq_p.resize (velocity_p.nVelChan_p);
607 0 : velocity_p.lsrFreq_p.resize (velocity_p.nVelChan_p);
608 0 : Double v0 = velocity_p.vStart_p.getValue ();
609 0 : Double dv = velocity_p.vInc_p.getValue ();
610 :
611 0 : if (aips_debug) {
612 0 : cout << "obsVel=" << obsVel << endl;
613 : }
614 :
615 0 : for (Int i = 0; i < velocity_p.nVelChan_p; i++) {
616 :
617 0 : Double vTopo = v0 + i * dv - obsVel;
618 0 : MDoppler dTopo (Quantity (vTopo, "m/s"), velocity_p.vDef_p);
619 0 : velocity_p.selFreq_p (i) = MFrequency::fromDoppler
620 0 : (dTopo, msIter_p.restFrequency ().getValue ()).getValue ().getValue ();
621 :
622 : // also calculate the frequencies in the requested frame for matching
623 : // up with the image planes
624 : // (they are called lsr here, but don't need to be in that frame)
625 :
626 0 : MDoppler dLSR (Quantity (v0 + i * dv, "m/s"), velocity_p.vDef_p);
627 0 : const MFrequency & restFrequency = msIter_p.restFrequency ();
628 0 : velocity_p.lsrFreq_p (i) = MFrequency::fromDoppler (dLSR, restFrequency.getValue ()).getValue ().getValue ();
629 : }
630 : }
631 285 : }
632 :
633 : void
634 0 : VisibilityIteratorReadImpl::getTopoFreqs (Vector<Double> & lsrFreq, Vector<Double> & selFreq)
635 : {
636 0 : getTopoFreqs ();
637 0 : lsrFreq.assign (velocity_p.lsrFreq_p);
638 0 : selFreq.assign (velocity_p.selFreq_p);
639 0 : }
640 :
641 :
642 :
643 : void
644 18 : VisibilityIteratorReadImpl::setState ()
645 : {
646 18 : if (stateOk_p) {
647 5 : return;
648 : }
649 :
650 13 : curTableNumRow_p = msIter_p.table ().nrow ();
651 : // get the times for this (major) iteration, so we can do (minor)
652 : // iteration by constant time (needed for VisBuffer averaging).
653 26 : ScalarColumn<Double> lcolTime (msIter_p.table (), MS::columnName (MS::TIME));
654 13 : time_p.resize (curTableNumRow_p);
655 13 : lcolTime.getColumn (time_p);
656 0 : ScalarColumn<Double> lcolTimeInterval (msIter_p.table (),
657 13 : MS::columnName (MS::INTERVAL));
658 : ///////////timeInterval_p.resize (curTableNumRow_p);
659 : ///////////lcolTimeInterval.getColumn (timeInterval_p);
660 13 : curStartRow_p = 0;
661 13 : setSelTable ();
662 13 : attachColumnsSafe (attachTable ());
663 : // If this is a new MeasurementSet then set up the antenna locations
664 13 : if (msIter_p.newMS ()) {
665 3 : nAnt_p = msd_p.setAntennas (msIter_p.msColumns ().antenna ());
666 3 : cache_p.feedpa_p.resize (nAnt_p);
667 3 : cache_p.feedpa_p.set (0);
668 3 : cache_p.lastfeedpaUT_p = -1;
669 3 : cache_p.parang_p.resize (nAnt_p);
670 3 : cache_p.parang_p.set (0);
671 3 : cache_p.lastParangUT_p = -1;
672 3 : cache_p.parang0_p = 0;
673 3 : cache_p.lastParang0UT_p = -1;
674 3 : cache_p.azel_p.resize (nAnt_p);
675 3 : cache_p.lastazelUT_p = -1;
676 3 : cache_p.lastazel0UT_p = -1;
677 3 : cache_p.lasthourangUT_p = -1;
678 :
679 : }
680 13 : if (msIter_p.newField () || msIterAtOrigin_p) {
681 13 : msd_p.setFieldCenter (msIter_p.phaseCenter ());
682 : }
683 13 : if ( msIter_p.newDataDescriptionId () || msIterAtOrigin_p) {
684 3 : Int spw = msIter_p.spectralWindowId ();
685 3 : nChan_p = msColumns ().spectralWindow ().numChan ()(spw);
686 3 : nPol_p = msColumns ().polarization ().numCorr ()(msIter_p.polarizationId ());
687 :
688 6 : if (Int (channels_p.nGroups_p.nelements ()) <= spw ||
689 3 : channels_p.nGroups_p[spw] == 0) {
690 : // no selection set yet, set default = all
691 : // for a reference MS this will normally be set appropriately in VisSet
692 0 : selectChannel (1, 0, nChan_p);
693 : }
694 3 : channelGroupSize_p = channels_p.width_p[spw];
695 3 : curNGroups_p = channels_p.nGroups_p[spw];
696 3 : cache_p.freqCacheOK_p = false;
697 : }
698 :
699 13 : stateOk_p = true;
700 : }
701 :
702 : const MSDerivedValues &
703 0 : VisibilityIteratorReadImpl::getMSD () const
704 : {
705 0 : return msd_p;
706 : }
707 :
708 :
709 : void
710 22 : VisibilityIteratorReadImpl::updateSlicer ()
711 : {
712 :
713 22 : if (msIter_p.newMS ()) {
714 12 : channels_p.nGroups_p.resize (0, true, false);
715 12 : doChannelSelection ();
716 : }
717 :
718 : // set the Slicer to get the selected part of spectrum out of the table
719 22 : Int spw = msIter_p.spectralWindowId ();
720 : //Fixed what i think was a confusion between chanWidth and chanInc
721 : // 2007/11/12
722 22 : Int start = channels_p.start_p[spw] + curChanGroup_p * channels_p.width_p[spw];
723 22 : AlwaysAssert (start >= 0 && start + channelGroupSize_p <= nChan_p, AipsError);
724 : // slicer_p=Slicer (Slice (),Slice (start,channelGroupSize_p));
725 : // above is slow, use IPositions instead.
726 44 : slicer_p = Slicer (IPosition (2, 0, start),
727 44 : IPosition (2, nPol_p, channelGroupSize_p),
728 66 : IPosition (2, 1, (channels_p.inc_p[spw] <= 0) ? 1 : channels_p.inc_p[spw] ));
729 44 : weightSlicer_p = Slicer (IPosition (1, start), IPosition (1, channelGroupSize_p),
730 66 : IPosition (1, (channels_p.inc_p[spw] <= 0) ? 1 : channels_p.inc_p[spw]));
731 22 : useSlicer_p = channelGroupSize_p < nChan_p;
732 :
733 : //if (msIter_p.newDataDescriptionId ()){
734 22 : setTileCache ();
735 : //}
736 22 : }
737 :
738 :
739 :
740 : void
741 30 : VisibilityIteratorReadImpl::setTileCache ()
742 : {
743 : // This function sets the tile cache because of a feature in
744 : // sliced data access that grows memory dramatically in some cases
745 : // if (useSlicer_p){
746 :
747 30 : if (autoTileCacheSizing_p){
748 10 : return; // rest of method does manual sizing so skip it
749 : }
750 :
751 30 : if (! (msIter_p.newDataDescriptionId () || msIter_p.newMS ()) ) {
752 10 : return;
753 : }
754 :
755 20 : const MeasurementSet & theMs = msIter_p.ms ();
756 20 : if (theMs.tableType () == Table::Memory) {
757 0 : return;
758 : }
759 :
760 20 : const ColumnDescSet & cds = theMs.tableDesc ().columnDescSet ();
761 :
762 20 : uInt startrow = msIter_p.table ().rowNumbers ()(0); // Get the first row number for this DDID.
763 :
764 20 : if (tileCacheIsSet_p.nelements () != 8) {
765 2 : tileCacheIsSet_p.resize (8);
766 2 : tileCacheIsSet_p.set (false);
767 : }
768 :
769 40 : Vector<String> columns (8);
770 20 : columns (0) = MS::columnName (MS::DATA); // complex
771 20 : columns (1) = MS::columnName (MS::CORRECTED_DATA); // complex
772 20 : columns (2) = MS::columnName (MS::MODEL_DATA); // complex
773 20 : columns (3) = MS::columnName (MS::FLAG); // boolean
774 20 : columns (4) = MS::columnName (MS::WEIGHT_SPECTRUM); // float
775 20 : columns (5) = MS::columnName (MS::WEIGHT); // float
776 20 : columns (6) = MS::columnName (MS::SIGMA); // float
777 20 : columns (7) = MS::columnName (MS::UVW); // double
778 :
779 180 : for (uInt k = 0; k < columns.nelements (); ++k) {
780 :
781 160 : if (! cds.isDefined (columns (k)) || ! usesTiledDataManager (columns[k], theMs)){
782 30 : continue;
783 : }
784 :
785 : try {
786 : //////////////////
787 : //////Temporary fix for virtual ms of multiple real ms's ...miracle of statics
788 : //////setting the cache size of hypercube at row 0 of each ms.
789 : ///will not work if each subms of a virtual ms has multi hypecube being
790 : ///accessed.
791 130 : if (theMs.tableInfo ().subType () == "CONCATENATED" &&
792 130 : msIterAtOrigin_p &&
793 0 : ! tileCacheIsSet_p[k]) {
794 :
795 0 : Block<String> refTables = theMs.getPartNames (true);
796 :
797 0 : for (uInt kk = 0; kk < refTables.nelements (); ++kk) {
798 :
799 0 : Table elms (refTables[kk]);
800 0 : ROTiledStManAccessor tacc (elms, columns[k], true);
801 :
802 : // Cleverly sense full-row cache size (in tiles)
803 0 : uInt cacheSizeInTiles(1);
804 :
805 : // Is startrow correct for each subMS?
806 : //uInt rRow=startrow;
807 0 : uInt rRow=0; // This preserves previous "single-hypercube" assumption
808 :
809 0 : IPosition hypercubeShape=tacc.hypercubeShape(rRow);
810 0 : IPosition tileShape=tacc.tileShape(rRow);
811 0 : uInt nax=hypercubeShape.size(); // how many axes
812 :
813 : // Accumulate axis factors up to--but NOT including--row (last) axis
814 : // "ceil" catches partially filled tiles...
815 0 : for (uInt iax=0;iax<nax-1;++iax)
816 0 : cacheSizeInTiles*= (uInt) ceil(hypercubeShape[iax]/ (Float)(tileShape[iax]));
817 :
818 : // Double for ALMA data
819 : // (this satisfies cases where required rows require >1 non-contiguous tiles)
820 0 : cacheSizeInTiles*=2;
821 :
822 : // Now set it
823 0 : tacc.setCacheSize (rRow, cacheSizeInTiles);
824 0 : tileCacheIsSet_p[k] = true;
825 : //cerr << "set cache on kk " << kk << " vol " << columns[k] << " " << refTables[kk] << endl;
826 : }
827 : }
828 : else {
829 :
830 260 : ROTiledStManAccessor tacc (theMs, columns[k], true);
831 :
832 :
833 : // Cleverly sense full-row cache size (in tiles)
834 130 : uInt cacheSizeInTiles(1);
835 :
836 260 : IPosition hypercubeShape=tacc.hypercubeShape(startrow);
837 260 : IPosition tileShape=tacc.tileShape(startrow);
838 130 : uInt nax=hypercubeShape.size(); // how many axes
839 :
840 : // Accumulate axis factors up to--but NOT including--row (last) axis
841 : // "ceil" catches partially filled tiles...
842 338 : for (uInt iax=0;iax<nax-1;++iax)
843 208 : cacheSizeInTiles*= (uInt) ceil(hypercubeShape[iax]/ (Float)(tileShape[iax]));
844 :
845 : // Double for ALMA data
846 : // (this satisfies cases where required rows require >1 non-contiguous tiles)
847 130 : cacheSizeInTiles*=2;
848 :
849 130 : Bool setCache = true;
850 130 : if (! useSlicer_p){ // always set cache if slicer in use
851 378 : for (uInt jj = 0 ; jj < tacc.nhypercubes (); ++jj) {
852 248 : if (tacc.getBucketSize (jj) == 0) {
853 118 : setCache = false;
854 : }
855 : }
856 : }
857 :
858 :
859 : /// If some bucketSize is 0...there is trouble in setting cache
860 : /// but if slicer is used it gushes anyways if one does not set cache
861 : /// need to fix the 0 bucket size in the filler anyways...then this is not needed
862 130 : if (setCache) {
863 : /*
864 : cout << columns[k]
865 : << " bucketSize=" << tacc.bucketSize(startrow)
866 : << " hcShape=" << hypercubeShape
867 : << " tShape=" << tileShape
868 : << " cacheSizeInTiles = " << cacheSizeInTiles << " (nhypercubes=" << tacc.nhypercubes() <<")"<< endl;
869 : */
870 12 : if (tacc.nhypercubes () == 1) {
871 12 : tacc.setCacheSize (0, cacheSizeInTiles);
872 : } else {
873 0 : tacc.setCacheSize (startrow, cacheSizeInTiles);
874 : }
875 : }
876 :
877 :
878 :
879 : }
880 : }
881 0 : catch (AipsError x) {
882 : // cerr << "Data man type " << dataManType << " " << dataManType.contains ("Tiled") << " && " << (!String (cdesc.dataManagerGroup ()).empty ()) << endl;
883 : // cerr << "Failed to set settilecache due to " << x.getMesg () << " column " << columns[k] <<endl;
884 : //It failed so leave the caching as is
885 0 : continue;
886 : }
887 : }
888 : }
889 :
890 : Bool
891 152 : VisibilityIteratorReadImpl::usesTiledDataManager (const String & columnName,
892 : const MeasurementSet & theMs) const
893 : {
894 152 : Bool noData = false;
895 :
896 : // Have to do something special about weight_spectrum as it tend to exist but
897 : // has no valid data.
898 :
899 304 : noData = noData ||
900 152 : (columnName == MS::columnName (MS::WEIGHT_SPECTRUM) && ! existsWeightSpectrum ());
901 :
902 : // Check to see if the column exist and have valid data
903 :
904 304 : noData = noData ||
905 172 : (columnName == MS::columnName (MS::DATA) &&
906 40 : (columns_p.vis_p.isNull () || ! columns_p.vis_p.isDefined (0)));
907 :
908 304 : noData = noData ||
909 172 : (columnName == MS::columnName (MS::MODEL_DATA) &&
910 40 : (columns_p.modelVis_p.isNull () || ! columns_p.modelVis_p.isDefined (0)));
911 :
912 297 : noData = noData ||
913 165 : (columnName == MS::columnName (MS::CORRECTED_DATA) &&
914 40 : (columns_p.corrVis_p.isNull () || ! columns_p.corrVis_p.isDefined (0)));
915 :
916 290 : noData = noData ||
917 158 : (columnName == MS::columnName (MS::FLAG) &&
918 40 : (columns_p.flag_p.isNull () || ! columns_p.flag_p.isDefined (0)));
919 :
920 290 : noData = noData ||
921 158 : (columnName == MS::columnName (MS::WEIGHT) &&
922 40 : (columns_p.weight_p.isNull () || ! columns_p.weight_p.isDefined (0)));
923 :
924 290 : noData = noData ||
925 158 : (columnName == MS::columnName (MS::SIGMA) &&
926 40 : (columns_p.sigma_p.isNull () || ! columns_p.sigma_p.isDefined (0)));
927 :
928 290 : noData = noData ||
929 158 : (columnName == MS::columnName (MS::UVW) &&
930 40 : (columns_p.uvw_p.isNull () || ! columns_p.uvw_p.isDefined (0)));
931 :
932 152 : Bool usesTiles = false;
933 :
934 152 : if (! noData){
935 276 : String dataManType = RODataManAccessor (theMs, columnName, true).dataManagerType ();
936 :
937 138 : usesTiles = dataManType.contains ("Tiled");
938 : }
939 :
940 152 : return usesTiles;
941 : }
942 :
943 :
944 : /*
945 : void VisibilityIteratorReadImpl::setTileCache (){
946 : // This function sets the tile cache because of a feature in
947 : // sliced data access that grows memory dramatically in some cases
948 : //if (useSlicer_p){
949 :
950 : {
951 : const MeasurementSet& thems=msIter_p.ms ();
952 : const ColumnDescSet& cds=thems.tableDesc ().columnDescSet ();
953 : ArrayColumn<Complex> columns_p.vis_p;
954 : ArrayColumn<Float> colwgt;
955 : Vector<String> columns (3);
956 : columns (0)=MS::columnName (MS::DATA);
957 : columns (1)=MS::columnName (MS::CORRECTED_DATA);
958 : columns (2)=MS::columnName (MS::MODEL_DATA);
959 : //cout << "COL " << columns << endl;
960 : for (uInt k=0; k< 3; ++k){
961 : //cout << "IN loop k " << k << endl;
962 : if (thems.tableDesc ().isColumn (columns (k)) ) {
963 :
964 : columns_p.vis_p.attach (thems,columns (k));
965 : String dataManType;
966 : dataManType = columns_p.vis_p.columnDesc ().dataManagerType ();
967 : //cout << "dataManType " << dataManType << endl;
968 : if (dataManType.contains ("Tiled")){
969 :
970 : ROTiledStManAccessor tacc (thems,
971 : columns_p.vis_p.columnDesc ().dataManagerGroup ());
972 : uInt nHyper = tacc.nhypercubes ();
973 : // Find smallest tile shape
974 : Int lowestProduct = 0;
975 : Int lowestId = 0;
976 : Bool firstFound = false;
977 : for (uInt id=0; id < nHyper; id++) {
978 : Int product = tacc.getTileShape (id).product ();
979 : if (product > 0 && (!firstFound || product < lowestProduct)) {
980 : lowestProduct = product;
981 : lowestId = id;
982 : if (!firstFound) firstFound = true;
983 : }
984 : }
985 : Int nchantile;
986 : IPosition tileshape=tacc.getTileShape (lowestId);
987 : IPosition axisPath (3,2,0,1);
988 : //nchantile=tileshape (1);
989 : tileshape (1)=channelGroupSize_p;
990 : tileshape (2)=curNumRow_p;
991 : //cout << "cursorshape " << tileshape << endl;
992 : nchantile=tacc.calcCacheSize (0, tileshape, axisPath);
993 :
994 : // if (nchantile > 0)
995 : // nchantile=channelGroupSize_p/nchantile*10;
996 : // if (nchantile<3)
997 : // nchantile=10;
998 :
999 : ///////////////
1000 : //nchantile *=8;
1001 : nchantile=1;
1002 : //tileshape (2)=tileshape (2)*8;
1003 : //////////////
1004 : //cout << tacc.cacheSize (0) << " nchantile "<< nchantile << " max cache size " << tacc.maximumCacheSize () << endl;
1005 : tacc.clearCaches ();
1006 : tacc.setCacheSize (0, 1);
1007 : //tacc.setCacheSize (0, tileshape, axisPath);
1008 : //cout << k << " " << columns (k) << " cache size " << tacc.cacheSize (0) << endl;
1009 :
1010 : }
1011 : }
1012 : }
1013 : }
1014 :
1015 : }
1016 : */
1017 :
1018 : void
1019 35 : VisibilityIteratorReadImpl::attachColumnsSafe (const Table & t)
1020 : {
1021 : // Normally, the call to attachColumns is redirected back to the ROVI class.
1022 : // This allows writable VIs to attach columns in both the read and write impls.
1023 : // However, this referral to the ROVI doesn't work during construction of this
1024 : // class (VIRI) since there is as yet no pointer to the object under construction.
1025 : // In that case, simply perform it locally.
1026 :
1027 35 : if (rovi_p == NULL){
1028 4 : attachColumns (t);
1029 : }
1030 : else{
1031 31 : rovi_p->attachColumns (t);
1032 : }
1033 35 : }
1034 :
1035 :
1036 : void
1037 35 : VisibilityIteratorReadImpl::attachColumns (const Table & t)
1038 : {
1039 35 : const ColumnDescSet & cds = t.tableDesc ().columnDescSet ();
1040 :
1041 35 : columns_p.antenna1_p.attach (t, MS::columnName (MS::ANTENNA1));
1042 35 : columns_p.antenna2_p.attach (t, MS::columnName (MS::ANTENNA2));
1043 :
1044 35 : if (cds.isDefined ("CORRECTED_DATA")) {
1045 35 : columns_p.corrVis_p.attach (t, "CORRECTED_DATA");
1046 : }
1047 :
1048 35 : columns_p.exposure_p.attach (t, MS::columnName (MS::EXPOSURE));
1049 35 : columns_p.feed1_p.attach (t, MS::columnName (MS::FEED1));
1050 35 : columns_p.feed2_p.attach (t, MS::columnName (MS::FEED2));
1051 35 : columns_p.flag_p.attach (t, MS::columnName (MS::FLAG));
1052 35 : columns_p.flagCategory_p.attach (t, MS::columnName (MS::FLAG_CATEGORY));
1053 35 : columns_p.flagRow_p.attach (t, MS::columnName (MS::FLAG_ROW));
1054 :
1055 35 : if (cds.isDefined (MS::columnName (MS::FLOAT_DATA))) {
1056 0 : columns_p.floatVis_p.attach (t, MS::columnName (MS::FLOAT_DATA));
1057 0 : floatDataFound_p = true;
1058 : } else {
1059 35 : floatDataFound_p = false;
1060 : }
1061 :
1062 35 : if (cds.isDefined ("MODEL_DATA")) {
1063 35 : columns_p.modelVis_p.attach (t, "MODEL_DATA");
1064 : }
1065 :
1066 35 : columns_p.observation_p.attach (t, MS::columnName (MS::OBSERVATION_ID));
1067 35 : columns_p.processor_p.attach (t, MS::columnName (MS::PROCESSOR_ID));
1068 35 : columns_p.scan_p.attach (t, MS::columnName (MS::SCAN_NUMBER));
1069 35 : columns_p.sigma_p.attach (t, MS::columnName (MS::SIGMA));
1070 35 : columns_p.state_p.attach (t, MS::columnName (MS::STATE_ID));
1071 35 : columns_p.time_p.attach (t, MS::columnName (MS::TIME));
1072 35 : columns_p.timeCentroid_p.attach (t, MS::columnName (MS::TIME_CENTROID));
1073 35 : columns_p.timeInterval_p.attach (t, MS::columnName (MS::INTERVAL));
1074 35 : columns_p.uvw_p.attach (t, MS::columnName (MS::UVW));
1075 :
1076 35 : if (cds.isDefined (MS::columnName (MS::DATA))) {
1077 35 : columns_p.vis_p.attach (t, MS::columnName (MS::DATA));
1078 : }
1079 :
1080 35 : columns_p.weight_p.attach (t, MS::columnName (MS::WEIGHT));
1081 :
1082 35 : if (cds.isDefined ("WEIGHT_SPECTRUM")) {
1083 29 : columns_p.weightSpectrum_p.attach (t, "WEIGHT_SPECTRUM");
1084 : }
1085 35 : }
1086 :
1087 : void
1088 0 : VisibilityIteratorReadImpl::update_rowIds () const
1089 : {
1090 0 : if (cache_p.rowIds_p.nelements () == 0) {
1091 0 : cache_p.rowIds_p = selRows_p.convert ();
1092 :
1093 0 : Vector<uInt> msIter_rowIds (msIter_p.table ().rowNumbers (msIter_p.ms ()));
1094 :
1095 0 : for (uInt i = 0; i < cache_p.rowIds_p.nelements (); i++) {
1096 0 : cache_p.rowIds_p (i) = msIter_rowIds (cache_p.rowIds_p (i));
1097 : }
1098 : }
1099 0 : return;
1100 : }
1101 :
1102 :
1103 : Int
1104 0 : VisibilityIteratorReadImpl::getDataDescriptionId () const
1105 : {
1106 0 : return msIter_p.dataDescriptionId ();
1107 : }
1108 :
1109 :
1110 : const MeasurementSet &
1111 0 : VisibilityIteratorReadImpl::getMeasurementSet () const
1112 : {
1113 0 : return msIter_p.ms ();
1114 : }
1115 :
1116 : Int
1117 0 : VisibilityIteratorReadImpl::getMeasurementSetId () const
1118 : {
1119 0 : return msIter_p.msId ();
1120 : }
1121 :
1122 :
1123 : Int
1124 0 : VisibilityIteratorReadImpl::getNAntennas () const
1125 : {
1126 0 : Int nAntennas = msIter_p.receptorAngle ().shape ()(1);
1127 :
1128 0 : return nAntennas;
1129 : }
1130 :
1131 : MEpoch
1132 0 : VisibilityIteratorReadImpl::getEpoch () const
1133 : {
1134 0 : MEpoch mEpoch = msIter_p.msColumns ().timeMeas ()(0);
1135 :
1136 0 : return mEpoch;
1137 : }
1138 :
1139 : Vector<Float>
1140 0 : VisibilityIteratorReadImpl::getReceptor0Angle ()
1141 : {
1142 0 : Int nAntennas = getNAntennas ();
1143 :
1144 0 : Vector<Float> receptor0Angle (nAntennas);
1145 :
1146 0 : for (int i = 0; i < nAntennas; i++) {
1147 0 : receptor0Angle [i] = msIter_p.receptorAngle ()(0, i);
1148 : }
1149 :
1150 0 : return receptor0Angle;
1151 : }
1152 :
1153 : Vector<rownr_t>
1154 0 : VisibilityIteratorReadImpl::getRowIds () const
1155 : {
1156 0 : update_rowIds ();
1157 :
1158 0 : return cache_p.rowIds_p;
1159 : }
1160 :
1161 :
1162 : Vector<rownr_t> &
1163 0 : VisibilityIteratorReadImpl::rowIds (Vector<rownr_t> & rowids) const
1164 : {
1165 : /* Calculate the row numbers in the original MS only when needed,
1166 : i.e. when this function is called */
1167 0 : update_rowIds ();
1168 0 : rowids.resize (cache_p.rowIds_p.nelements ());
1169 0 : rowids = cache_p.rowIds_p;
1170 0 : return rowids;
1171 : }
1172 :
1173 :
1174 : Vector<Int> &
1175 200 : VisibilityIteratorReadImpl::antenna1(Vector<Int> & ant1) const
1176 : {
1177 200 : ant1.resize (curNumRow_p);
1178 200 : getCol (columns_p.antenna1_p, ant1);
1179 200 : return ant1;
1180 : }
1181 :
1182 : Vector<Int> &
1183 200 : VisibilityIteratorReadImpl::antenna2(Vector<Int> & ant2) const
1184 : {
1185 200 : ant2.resize (curNumRow_p);
1186 200 : getCol (columns_p.antenna2_p, ant2);
1187 200 : return ant2;
1188 : }
1189 :
1190 : Vector<Int> &
1191 0 : VisibilityIteratorReadImpl::feed1(Vector<Int> & fd1) const
1192 : {
1193 0 : fd1.resize (curNumRow_p);
1194 0 : getCol (columns_p.feed1_p, fd1);
1195 0 : return fd1;
1196 : }
1197 :
1198 : Vector<Int> &
1199 0 : VisibilityIteratorReadImpl::feed2(Vector<Int> & fd2) const
1200 : {
1201 0 : fd2.resize (curNumRow_p);
1202 0 : getCol (columns_p.feed2_p, fd2);
1203 0 : return fd2;
1204 : }
1205 :
1206 : Vector<Int> &
1207 260 : VisibilityIteratorReadImpl::channel (Vector<Int> & chan) const
1208 : {
1209 260 : Int spw = msIter_p.spectralWindowId ();
1210 260 : chan.resize (channelGroupSize_p);
1211 260 : Int inc = channels_p.inc_p[spw] <= 0 ? 1 : channels_p.inc_p[spw];
1212 24040 : for (Int i = 0; i < channelGroupSize_p; i++) {
1213 23780 : chan (i) = channels_p.start_p[spw] + curChanGroup_p * channels_p.width_p[spw] + i * inc;
1214 : }
1215 260 : return chan;
1216 : }
1217 :
1218 : Vector<Int> &
1219 207 : VisibilityIteratorReadImpl::corrType (Vector<Int> & corrTypes) const
1220 : {
1221 207 : Int polId = msIter_p.polarizationId ();
1222 207 : msIter_p.msColumns ().polarization ().corrType ().get (polId, corrTypes, true);
1223 207 : return corrTypes;
1224 : }
1225 :
1226 : Cube<Bool> &
1227 260 : VisibilityIteratorReadImpl::flag (Cube<Bool> & flags) const
1228 : {
1229 260 : if (useSlicer_p) {
1230 0 : getCol (columns_p.flag_p, slicer_p, flags, true);
1231 : } else {
1232 260 : getCol (columns_p.flag_p, flags, true);
1233 : }
1234 260 : return flags;
1235 : }
1236 :
1237 : Matrix<Bool> &
1238 200 : VisibilityIteratorReadImpl::flag (Matrix<Bool> & flags) const
1239 : {
1240 200 : if (useSlicer_p) {
1241 0 : getCol (columns_p.flag_p, slicer_p, cache_p.flagCube_p, true);
1242 : } else {
1243 200 : getCol (columns_p.flag_p, cache_p.flagCube_p, true);
1244 : }
1245 :
1246 200 : flags.resize (channelGroupSize_p, curNumRow_p);
1247 : // need to optimize this...
1248 : //for (Int row=0; row<curNumRow_p; row++) {
1249 : // for (Int chn=0; chn<channelGroupSize_p; chn++) {
1250 : // flags (chn,row)=flagCube (0,chn,row);
1251 : // for (Int pol=1; pol<nPol_p; pol++) {
1252 : // flags (chn,row)|=flagCube (pol,chn,row);
1253 : // }
1254 : // }
1255 : //}
1256 : Bool deleteIt1;
1257 : Bool deleteIt2;
1258 200 : const Bool * pcube = cache_p.flagCube_p.getStorage (deleteIt1);
1259 200 : Bool * pflags = flags.getStorage (deleteIt2);
1260 70400 : for (uInt row = 0; row < curNumRow_p; row++) {
1261 7090200 : for (Int chn = 0; chn < channelGroupSize_p; chn++) {
1262 7020000 : *pflags = *pcube++;
1263 28080000 : for (Int pol = 1; pol < nPol_p; pol++, pcube++) {
1264 21060000 : *pflags = *pcube ? *pcube : *pflags;
1265 : }
1266 7020000 : pflags++;
1267 : }
1268 : }
1269 200 : cache_p.flagCube_p.freeStorage (pcube, deleteIt1);
1270 200 : flags.putStorage (pflags, deleteIt2);
1271 200 : return flags;
1272 : }
1273 :
1274 0 : Bool VisibilityIteratorReadImpl::existsFlagCategory() const
1275 : {
1276 0 : if(msIter_p.newMS()){ // Cache to avoid testing unnecessarily.
1277 : try{
1278 0 : cache_p.msHasFC_p = columns_p.flagCategory_p.hasContent();
1279 : }
1280 0 : catch (AipsError x){
1281 0 : cache_p.msHasFC_p = false;
1282 : }
1283 : }
1284 0 : return cache_p.msHasFC_p;
1285 : }
1286 :
1287 : Array<Bool> &
1288 0 : VisibilityIteratorReadImpl::flagCategory (Array<Bool> & flagCategories) const
1289 : {
1290 0 : if (columns_p.flagCategory_p.isNull () || !columns_p.flagCategory_p.isDefined (0)) { // It often is.
1291 0 : flagCategories.resize (); // Zap it.
1292 : } else {
1293 0 : if (velocity_p.velSelection_p) {
1294 0 : throw (AipsError ("velocity selection not allowed in flagCategory ()."));
1295 : } else {
1296 0 : if (useSlicer_p) {
1297 0 : getCol (columns_p.flagCategory_p, slicer_p, flagCategories, true);
1298 : } else {
1299 0 : getCol (columns_p.flagCategory_p, flagCategories, true);
1300 : }
1301 : }
1302 : }
1303 0 : return flagCategories;
1304 : }
1305 :
1306 : Vector<Bool> &
1307 200 : VisibilityIteratorReadImpl::flagRow (Vector<Bool> & rowflags) const
1308 : {
1309 200 : rowflags.resize (curNumRow_p);
1310 200 : getCol (columns_p.flagRow_p, rowflags);
1311 200 : return rowflags;
1312 : }
1313 :
1314 : Vector<Int> &
1315 0 : VisibilityIteratorReadImpl::observationId (Vector<Int> & obsIDs) const
1316 : {
1317 0 : obsIDs.resize (curNumRow_p);
1318 0 : getCol (columns_p.observation_p, obsIDs);
1319 0 : return obsIDs;
1320 : }
1321 :
1322 : Vector<Int> &
1323 0 : VisibilityIteratorReadImpl::processorId (Vector<Int> & procIDs) const
1324 : {
1325 0 : procIDs.resize (curNumRow_p);
1326 0 : getCol (columns_p.processor_p, procIDs);
1327 0 : return procIDs;
1328 : }
1329 :
1330 : Vector<Int> &
1331 0 : VisibilityIteratorReadImpl::scan (Vector<Int> & scans) const
1332 : {
1333 0 : scans.resize (curNumRow_p);
1334 0 : getCol (columns_p.scan_p, scans);
1335 0 : return scans;
1336 : }
1337 :
1338 : Vector<Int> &
1339 0 : VisibilityIteratorReadImpl::stateId (Vector<Int> & stateIds) const
1340 : {
1341 0 : stateIds.resize (curNumRow_p);
1342 0 : getCol (columns_p.state_p, stateIds);
1343 0 : return stateIds;
1344 : }
1345 :
1346 : Vector<Double> &
1347 201 : VisibilityIteratorReadImpl::frequency (Vector<Double> & freq) const
1348 : {
1349 201 : if (velocity_p.velSelection_p) {
1350 0 : freq.resize (velocity_p.nVelChan_p);
1351 0 : freq = velocity_p.selFreq_p;
1352 : } else {
1353 201 : if (! cache_p.freqCacheOK_p) {
1354 2 : cache_p.freqCacheOK_p = true;
1355 2 : Int spw = msIter_p.spectralWindowId ();
1356 2 : cache_p.frequency_p.resize (channelGroupSize_p);
1357 2 : const Vector<Double> & chanFreq = msIter_p.frequency ();
1358 2 : Int start = channels_p.start_p[spw];
1359 2 : Int inc = channels_p.inc_p[spw] <= 0 ? 1 : channels_p.inc_p[spw];
1360 202 : for (Int i = 0; i < channelGroupSize_p; i++) {
1361 200 : cache_p.frequency_p (i) = chanFreq (start + curChanGroup_p * channels_p.width_p[spw] + i * inc);
1362 : }
1363 : }
1364 201 : freq.resize (channelGroupSize_p);
1365 201 : freq = cache_p.frequency_p;
1366 : }
1367 201 : return freq;
1368 : }
1369 :
1370 :
1371 : Vector<Double> &
1372 201 : VisibilityIteratorReadImpl::time (Vector<Double> & t) const
1373 : {
1374 201 : t.resize (curNumRow_p);
1375 :
1376 201 : getCol (columns_p.time_p, t);
1377 :
1378 201 : return t;
1379 : }
1380 :
1381 : Vector<Double> &
1382 0 : VisibilityIteratorReadImpl::timeCentroid (Vector<Double> & t) const
1383 : {
1384 0 : t.resize (curNumRow_p);
1385 0 : getCol (columns_p.timeCentroid_p, t);
1386 0 : return t;
1387 : }
1388 :
1389 : Vector<Double> &
1390 0 : VisibilityIteratorReadImpl::timeInterval (Vector<Double> & t) const
1391 : {
1392 0 : t.resize (curNumRow_p);
1393 0 : getCol (columns_p.timeInterval_p, t);
1394 0 : return t;
1395 : }
1396 :
1397 : Vector<Double> &
1398 0 : VisibilityIteratorReadImpl::exposure (Vector<Double> & expo) const
1399 : {
1400 0 : expo.resize (curNumRow_p);
1401 0 : getCol (columns_p.exposure_p, expo);
1402 0 : return expo;
1403 : }
1404 :
1405 : Cube<Complex> &
1406 266 : VisibilityIteratorReadImpl::visibility (Cube<Complex> & vis, DataColumn whichOne) const
1407 : {
1408 :
1409 266 : if (useSlicer_p) {
1410 0 : getDataColumn (whichOne, slicer_p, vis);
1411 : } else {
1412 266 : getDataColumn (whichOne, vis);
1413 : }
1414 :
1415 266 : return vis;
1416 : }
1417 :
1418 :
1419 : // helper function to swap the y and z axes of a Cube
1420 : void
1421 0 : swapyz (Cube<Complex> & out, const Cube<Complex> & in)
1422 : {
1423 0 : IPosition inShape = in.shape ();
1424 0 : uInt nx = inShape (0), ny = inShape (2), nz = inShape (1);
1425 0 : out.resize (nx, ny, nz);
1426 : Bool deleteIn, deleteOut;
1427 0 : const Complex * pin = in.getStorage (deleteIn);
1428 0 : Complex * pout = out.getStorage (deleteOut);
1429 0 : uInt i = 0, zOffset = 0;
1430 0 : for (uInt iz = 0; iz < nz; iz++, zOffset += nx) {
1431 0 : Int yOffset = zOffset;
1432 0 : for (uInt iy = 0; iy < ny; iy++, yOffset += nx * nz) {
1433 0 : for (uInt ix = 0; ix < nx; ix++) {
1434 0 : pout[i++] = pin[ix + yOffset];
1435 : }
1436 : }
1437 : }
1438 0 : out.putStorage (pout, deleteOut);
1439 0 : in.freeStorage (pin, deleteIn);
1440 0 : }
1441 :
1442 : // helper function to swap the y and z axes of a Cube
1443 : void
1444 0 : swapyz (Cube<Bool> & out, const Cube<Bool> & in)
1445 : {
1446 0 : IPosition inShape = in.shape ();
1447 0 : uInt nx = inShape (0), ny = inShape (2), nz = inShape (1);
1448 0 : out.resize (nx, ny, nz);
1449 : Bool deleteIn, deleteOut;
1450 0 : const Bool * pin = in.getStorage (deleteIn);
1451 0 : Bool * pout = out.getStorage (deleteOut);
1452 0 : uInt i = 0, zOffset = 0;
1453 0 : for (uInt iz = 0; iz < nz; iz++, zOffset += nx) {
1454 0 : Int yOffset = zOffset;
1455 0 : for (uInt iy = 0; iy < ny; iy++, yOffset += nx * nz) {
1456 0 : for (uInt ix = 0; ix < nx; ix++) {
1457 0 : pout[i++] = pin[ix + yOffset];
1458 : }
1459 : }
1460 : }
1461 0 : }
1462 :
1463 : Cube<Float> &
1464 0 : VisibilityIteratorReadImpl::floatData (Cube<Float> & fcube) const
1465 : {
1466 0 : if (useSlicer_p) {
1467 0 : getFloatDataColumn (slicer_p, fcube);
1468 : } else {
1469 0 : getFloatDataColumn (fcube);
1470 : }
1471 0 : return fcube;
1472 : }
1473 :
1474 : // transpose a matrix
1475 : void
1476 0 : transpose (Matrix<Float> & out, const Matrix<Float> & in)
1477 : {
1478 0 : uInt ny = in.nrow (), nx = in.ncolumn ();
1479 0 : out.resize (nx, ny);
1480 : Bool deleteIn, deleteOut;
1481 0 : const Float * pin = in.getStorage (deleteIn);
1482 0 : Float * pout = out.getStorage (deleteOut);
1483 0 : uInt i = 0;
1484 0 : for (uInt iy = 0; iy < ny; iy++) {
1485 0 : uInt yOffset = 0;
1486 0 : for (uInt ix = 0; ix < nx; ix++, yOffset += ny) {
1487 0 : pout[i++] = pin[iy + yOffset];
1488 : }
1489 : }
1490 0 : out.putStorage (pout, deleteOut);
1491 0 : in.freeStorage (pin, deleteIn);
1492 0 : }
1493 :
1494 : void
1495 0 : VisibilityIteratorReadImpl::getDataColumn (DataColumn whichOne,
1496 : const Slicer & slicer,
1497 : Cube<Complex> & data) const
1498 : {
1499 :
1500 :
1501 : // Return the visibility (observed, model or corrected);
1502 : // deal with DATA and FLOAT_DATA seamlessly for observed data.
1503 0 : switch (whichOne) {
1504 :
1505 0 : case ROVisibilityIterator::Observed:
1506 0 : if (floatDataFound_p) {
1507 0 : Cube<Float> dataFloat;
1508 0 : getCol (columns_p.floatVis_p, slicer, dataFloat, true);
1509 0 : data.resize (dataFloat.shape ());
1510 0 : convertArray (data, dataFloat);
1511 : } else {
1512 0 : getCol (columns_p.vis_p, slicer, data, true);
1513 : }
1514 0 : break;
1515 :
1516 0 : case ROVisibilityIterator::Corrected:
1517 0 : getCol (columns_p.corrVis_p, slicer, data, true);
1518 0 : break;
1519 :
1520 0 : case ROVisibilityIterator::Model:
1521 0 : getCol (columns_p.modelVis_p, slicer, data, true);
1522 0 : break;
1523 :
1524 0 : default:
1525 0 : Assert (false);
1526 : }
1527 :
1528 0 : }
1529 :
1530 : void
1531 266 : VisibilityIteratorReadImpl::getDataColumn (DataColumn whichOne,
1532 : Cube<Complex> & data) const
1533 : {
1534 : // Return the visibility (observed, model or corrected);
1535 : // deal with DATA and FLOAT_DATA seamlessly for observed data.
1536 :
1537 266 : switch (whichOne) {
1538 :
1539 6 : case ROVisibilityIterator::Observed:
1540 6 : if (floatDataFound_p) {
1541 0 : Cube<Float> dataFloat;
1542 0 : getCol (columns_p.floatVis_p, dataFloat, true);
1543 0 : data.resize (dataFloat.shape ());
1544 0 : convertArray (data, dataFloat);
1545 : } else {
1546 6 : getCol (columns_p.vis_p, data, true);
1547 : }
1548 6 : break;
1549 :
1550 60 : case ROVisibilityIterator::Corrected:
1551 60 : getCol (columns_p.corrVis_p, data, true);
1552 60 : break;
1553 :
1554 200 : case ROVisibilityIterator::Model:
1555 200 : getCol (columns_p.modelVis_p, data, true);
1556 200 : break;
1557 :
1558 0 : default:
1559 0 : Assert (false);
1560 : }
1561 266 : }
1562 :
1563 : void
1564 0 : VisibilityIteratorReadImpl::getFloatDataColumn (const Slicer & slicer,
1565 : Cube<Float> & data) const
1566 : {
1567 : // Return FLOAT_DATA as real Floats.
1568 0 : if (floatDataFound_p) {
1569 0 : getCol (columns_p.floatVis_p, slicer, data, true);
1570 : }
1571 0 : }
1572 :
1573 : void
1574 0 : VisibilityIteratorReadImpl::getFloatDataColumn (Cube<Float> & data) const
1575 : {
1576 : // Return FLOAT_DATA as real Floats.
1577 0 : if (floatDataFound_p) {
1578 0 : getCol (columns_p.floatVis_p, data, true);
1579 : }
1580 0 : }
1581 :
1582 : Matrix<CStokesVector> &
1583 0 : VisibilityIteratorReadImpl::visibility (Matrix<CStokesVector> & vis,
1584 : DataColumn whichOne) const
1585 : {
1586 0 : if (useSlicer_p) {
1587 0 : getDataColumn (whichOne, slicer_p, cache_p.visCube_p);
1588 : } else {
1589 0 : getDataColumn (whichOne, cache_p.visCube_p);
1590 : }
1591 :
1592 0 : vis.resize (channelGroupSize_p, curNumRow_p);
1593 : Bool deleteIt;
1594 0 : Complex * pcube = cache_p.visCube_p.getStorage (deleteIt);
1595 0 : if (deleteIt) {
1596 0 : cerr << "Problem in ROVisIter::visibility - deleteIt true" << endl;
1597 : }
1598 : // Here we cope in a limited way with cases where not all 4
1599 : // polarizations are present: if only 2, assume XX,YY or RR,LL
1600 : // if only 1, assume it's an estimate of Stokes I (one of RR,LL,XX,YY)
1601 : // The cross terms are zero filled in these cases.
1602 0 : switch (nPol_p) {
1603 0 : case 4: {
1604 0 : for (uInt row = 0; row < curNumRow_p; row++) {
1605 0 : for (Int chn = 0; chn < channelGroupSize_p; chn++, pcube += 4) {
1606 0 : vis (chn, row) = pcube;
1607 : }
1608 : }
1609 0 : break;
1610 : }
1611 0 : case 2: {
1612 0 : vis.set (CStokesVector (Complex (0., 0.)));
1613 0 : for (uInt row = 0; row < curNumRow_p; row++) {
1614 0 : for (Int chn = 0; chn < channelGroupSize_p; chn++, pcube += 2) {
1615 0 : CStokesVector & v = vis (chn, row);
1616 0 : v (0) = *pcube;
1617 0 : v (3) = *(pcube + 1);
1618 : }
1619 : }
1620 0 : break;
1621 : }
1622 0 : case 1: {
1623 0 : vis.set (CStokesVector (Complex (0., 0.)));
1624 0 : for (uInt row = 0; row < curNumRow_p; row++) {
1625 0 : for (Int chn = 0; chn < channelGroupSize_p; chn++, pcube++) {
1626 0 : CStokesVector & v = vis (chn, row);
1627 0 : v (0) = v (3) = *pcube;
1628 : }
1629 : }
1630 : } //# case 1
1631 : } //# switch
1632 0 : return vis;
1633 : }
1634 :
1635 : Vector<RigidVector<Double, 3> > &
1636 200 : VisibilityIteratorReadImpl::uvw (Vector<RigidVector<Double, 3> > & uvwvec) const
1637 : {
1638 200 : uvwvec.resize (curNumRow_p);
1639 200 : getColArray<Double>(columns_p.uvw_p, cache_p.uvwMat_p, true);
1640 : // get a pointer to the raw storage for quick access
1641 : Bool deleteIt;
1642 200 : Double * pmat = cache_p.uvwMat_p.getStorage (deleteIt);
1643 70400 : for (uInt row = 0; row < curNumRow_p; row++, pmat += 3) {
1644 70200 : uvwvec (row) = pmat;
1645 : }
1646 200 : return uvwvec;
1647 : }
1648 :
1649 : Matrix<Double> &
1650 0 : VisibilityIteratorReadImpl::uvwMat (Matrix<Double> & uvwmat) const
1651 : {
1652 0 : getCol (columns_p.uvw_p, uvwmat, true);
1653 0 : return uvwmat;
1654 : }
1655 :
1656 : // Fill in parallactic angle.
1657 : Vector<Float>
1658 0 : VisibilityIteratorReadImpl::feed_pa (Double time) const
1659 : {
1660 : // LogMessage message (LogOrigin ("VisibilityIteratorReadImpl","feed_pa"));
1661 :
1662 : // Absolute UT
1663 0 : Double ut = time;
1664 :
1665 0 : if (ut != cache_p.lastfeedpaUT_p) {
1666 :
1667 : // Set up the Epoch using the absolute MJD in seconds
1668 : // get the Epoch reference from the column
1669 :
1670 0 : MEpoch mEpoch = getEpoch ();
1671 :
1672 0 : const Matrix<Double> & angles = receptorAngles ().xyPlane(0);
1673 0 : Int nAnt = angles.shape ()(1);
1674 :
1675 0 : Vector<Float> receptor0Angle (nAnt, 0);
1676 :
1677 0 : for (int i = 0; i < nAnt; i++) {
1678 0 : receptor0Angle [i] = angles (0, i);
1679 : }
1680 :
1681 0 : cache_p.feedpa_p.assign (feed_paCalculate (time, msd_p, nAnt, mEpoch, receptor0Angle));
1682 :
1683 0 : cache_p.lastfeedpaUT_p = ut;
1684 : }
1685 0 : return cache_p.feedpa_p;
1686 : }
1687 :
1688 : Vector<Float>
1689 0 : VisibilityIteratorReadImpl::feed_paCalculate (Double time, MSDerivedValues & msd,
1690 : Int nAntennas, const MEpoch & mEpoch0,
1691 : const Vector<Float> & receptor0Angle)
1692 : {
1693 0 : MEpoch mEpoch = mEpoch0;
1694 :
1695 0 : mEpoch.set (MVEpoch (Quantity (time, "s")));
1696 :
1697 0 : msd.setEpoch (mEpoch);
1698 :
1699 : // Calculate pa for all antennas.
1700 :
1701 0 : Vector<Float> feedpa (nAntennas);
1702 :
1703 0 : for (Int iant = 0; iant < nAntennas; iant++) {
1704 :
1705 0 : msd.setAntenna (iant);
1706 0 : feedpa (iant) = msd.parAngle ();
1707 :
1708 : // add angle for receptor 0
1709 :
1710 0 : feedpa (iant) += receptor0Angle (iant);
1711 :
1712 0 : if (aips_debug && iant == 0) {
1713 :
1714 0 : cout << "Antenna " << iant << " at time: " << MVTime (mEpoch.getValue ()) <<
1715 0 : " has PA = " << feedpa (iant) * 57.28 << endl;
1716 : }
1717 : }
1718 :
1719 0 : return feedpa;
1720 : }
1721 :
1722 : // Fill in parallactic angle.
1723 : const Float &
1724 0 : VisibilityIteratorReadImpl::parang0(Double time) const
1725 : {
1726 : // LogMessage message (LogOrigin ("VisibilityIteratorReadImpl","parang0"));
1727 :
1728 : // Absolute UT
1729 0 : Double ut = time;
1730 :
1731 0 : if (ut != cache_p.lastParang0UT_p) {
1732 :
1733 0 : cache_p.lastParang0UT_p = ut;
1734 :
1735 : // Set up the Epoch using the absolute MJD in seconds
1736 : // get the Epoch reference from the column
1737 0 : MEpoch mEpoch = msIter_p.msColumns ().timeMeas ()(0);
1738 0 : cache_p.parang0_p = parang0Calculate (time, msd_p, mEpoch);
1739 : }
1740 0 : return cache_p.parang0_p;
1741 : }
1742 :
1743 : Float
1744 0 : VisibilityIteratorReadImpl::parang0Calculate (Double time, MSDerivedValues & msd, const MEpoch & mEpoch0)
1745 : {
1746 0 : MEpoch mEpoch = mEpoch0;
1747 :
1748 0 : mEpoch.set (MVEpoch (Quantity (time, "s")));
1749 0 : msd.setEpoch (mEpoch);
1750 :
1751 : // Calculate pa for all antennas.
1752 0 : msd.setAntenna (-1);
1753 0 : Float parang0 = msd.parAngle ();
1754 :
1755 0 : if (aips_debug)
1756 0 : cout << "At time: " << MVTime (mEpoch.getValue ()) <<
1757 0 : " PA = " << parang0 * 57.28 << " deg" << endl;
1758 :
1759 0 : return parang0;
1760 : }
1761 :
1762 :
1763 : // Fill in parallactic angle (NO FEED PA offset!).
1764 : Vector<Float>
1765 0 : VisibilityIteratorReadImpl::parang (Double time) const
1766 : {
1767 : // LogMessage message (LogOrigin ("VisibilityIteratorReadImpl","parang"));
1768 :
1769 : // Absolute UT
1770 0 : Double ut = time;
1771 :
1772 0 : if (ut != cache_p.lastParangUT_p) {
1773 :
1774 0 : cache_p.lastParangUT_p = ut;
1775 :
1776 : // Set up the Epoch using the absolute MJD in seconds
1777 : // get the Epoch reference from the column
1778 :
1779 0 : MEpoch mEpoch = msIter_p.msColumns ().timeMeas ()(0);
1780 0 : Int nAnt = msIter_p.receptorAngle ().shape ()(1);
1781 :
1782 0 : cache_p.parang_p = parangCalculate (time, msd_p, nAnt, mEpoch);
1783 :
1784 : }
1785 0 : return cache_p.parang_p;
1786 : }
1787 :
1788 : Vector<Float>
1789 0 : VisibilityIteratorReadImpl::parangCalculate (Double time, MSDerivedValues & msd, int nAntennas, const MEpoch mEpoch0)
1790 : {
1791 0 : MEpoch mEpoch = mEpoch0;
1792 0 : mEpoch.set (MVEpoch (Quantity (time, "s")));
1793 :
1794 0 : msd.setEpoch (mEpoch);
1795 :
1796 : // Calculate pa for all antennas.
1797 :
1798 0 : Vector<Float> parang (nAntennas);
1799 :
1800 0 : for (Int iant = 0; iant < nAntennas; iant++) {
1801 :
1802 0 : msd.setAntenna (iant);
1803 0 : parang (iant) = msd.parAngle ();
1804 :
1805 0 : if (aips_debug && iant == 0) {
1806 0 : cout << "Antenna " << iant << " at time: " << MVTime (mEpoch.getValue ()) <<
1807 0 : " has PA = " << parang (iant) * 57.28 << endl;
1808 : }
1809 : }
1810 :
1811 0 : return parang;
1812 : }
1813 :
1814 : // Fill in azimuth/elevation of the antennas.
1815 : // Cloned from feed_pa, we need to check that this is all correct!
1816 : Vector<MDirection>
1817 0 : VisibilityIteratorReadImpl::azel (Double time) const
1818 : {
1819 : // LogMessage message (LogOrigin ("VisibilityIteratorReadImpl","azel"));
1820 :
1821 : // Absolute UT
1822 0 : Double ut = time;
1823 :
1824 0 : if (ut != cache_p.lastazelUT_p) {
1825 :
1826 0 : cache_p.lastazelUT_p = ut;
1827 0 : Int nAnt = msIter_p.receptorAngle ().shape ()(1);
1828 0 : MEpoch mEpoch = msIter_p.msColumns ().timeMeas ()(0);
1829 :
1830 0 : azelCalculate (ut, msd_p, cache_p.azel_p, nAnt, mEpoch);
1831 :
1832 : }
1833 0 : return cache_p.azel_p;
1834 : }
1835 :
1836 : void
1837 0 : VisibilityIteratorReadImpl::azelCalculate (Double time, MSDerivedValues & msd, Vector<MDirection> & azel,
1838 : Int nAnt, const MEpoch & mEpoch0)
1839 : {
1840 : // Refactored into a static method to allow VisBufferAsync to use
1841 :
1842 0 : MEpoch mEpoch = mEpoch0;
1843 :
1844 0 : mEpoch.set (MVEpoch (Quantity (time, "s")));
1845 :
1846 0 : msd.setEpoch (mEpoch);
1847 :
1848 : // Calculate az/el for all antennas.
1849 :
1850 0 : azel.resize (nAnt);
1851 :
1852 0 : for (Int iant = 0; iant < nAnt; iant++) {
1853 0 : msd.setAntenna (iant);
1854 0 : azel (iant) = msd.azel ();
1855 0 : if (aips_debug) {
1856 0 : if (iant == 0)
1857 0 : cout << "Antenna " << iant << " at time: " << MVTime (mEpoch.getValue ()) <<
1858 0 : " has AzEl = " << azel (iant).getAngle ("deg") << endl;
1859 : }
1860 : }
1861 0 : }
1862 :
1863 : // Fill in azimuth/elevation of the antennas.
1864 : // Cloned from feed_pa, we need to check that this is all correct!
1865 : MDirection
1866 0 : VisibilityIteratorReadImpl::azel0(Double time) const
1867 : {
1868 : // LogMessage message (LogOrigin ("VisibilityIteratorReadImpl","azel0"));
1869 :
1870 : // Absolute UT
1871 0 : Double ut = time;
1872 :
1873 0 : if (ut != cache_p.lastazel0UT_p) {
1874 :
1875 0 : cache_p.lastazel0UT_p = ut;
1876 :
1877 0 : MEpoch mEpoch = msIter_p.msColumns ().timeMeas ()(0);
1878 :
1879 0 : azel0Calculate (time, msd_p, cache_p.azel0_p, mEpoch);
1880 :
1881 : }
1882 0 : return cache_p.azel0_p;
1883 : }
1884 :
1885 : void
1886 0 : VisibilityIteratorReadImpl::azel0Calculate (Double time, MSDerivedValues & msd,
1887 : MDirection & azel0, const MEpoch & mEpoch0)
1888 : {
1889 : // Refactored into a static method to allow VisBufferAsync to use
1890 :
1891 0 : MEpoch mEpoch = mEpoch0;
1892 :
1893 0 : mEpoch.set (MVEpoch (Quantity (time, "s")));
1894 :
1895 0 : msd.setEpoch (mEpoch);
1896 :
1897 0 : msd.setAntenna (-1);
1898 :
1899 0 : azel0 = msd.azel ();
1900 :
1901 0 : if (aips_debug) {
1902 0 : cout << "At time: " << MVTime (mEpoch.getValue ()) <<
1903 0 : " AzEl = " << azel0.getAngle ("deg") << endl;
1904 : }
1905 :
1906 0 : }
1907 :
1908 : // Hour angle at specified time.
1909 : Double
1910 0 : VisibilityIteratorReadImpl::hourang (Double time) const
1911 : {
1912 : // LogMessage message (LogOrigin ("VisibilityIteratorReadImpl","azel"));
1913 :
1914 : // Absolute UT
1915 0 : Double ut = time;
1916 :
1917 0 : if (ut != cache_p.lasthourangUT_p) {
1918 :
1919 0 : cache_p.lasthourangUT_p = ut;
1920 :
1921 : // Set up the Epoch using the absolute MJD in seconds
1922 : // get the Epoch reference from the column keyword
1923 :
1924 0 : MEpoch mEpoch = msIter_p.msColumns ().timeMeas ()(0);
1925 :
1926 0 : cache_p.hourang_p = hourangCalculate (time, msd_p, mEpoch);
1927 :
1928 : }
1929 0 : return cache_p.hourang_p;
1930 : }
1931 :
1932 : Double
1933 0 : VisibilityIteratorReadImpl::hourangCalculate (Double time, MSDerivedValues & msd, const MEpoch & mEpoch0)
1934 : {
1935 0 : MEpoch mEpoch = mEpoch0;
1936 :
1937 0 : mEpoch.set (MVEpoch (Quantity (time, "s")));
1938 :
1939 0 : msd.setEpoch (mEpoch);
1940 :
1941 0 : msd.setAntenna (-1);
1942 :
1943 0 : Double hourang = msd.hourAngle ();
1944 :
1945 0 : return hourang;
1946 : }
1947 :
1948 : Vector<Float> &
1949 0 : VisibilityIteratorReadImpl::sigma (Vector<Float> & sig) const
1950 : {
1951 0 : Matrix<Float> sigmat;
1952 0 : getCol (columns_p.sigma_p, sigmat);
1953 : // Do a rough average of the parallel hand polarizations to get a single
1954 : // sigma. Should do this properly someday, or return all values
1955 0 : sig.resize (sigmat.ncolumn ());
1956 0 : sig = sigmat.row (0);
1957 0 : sig += sigmat.row (nPol_p - 1);
1958 0 : sig /= 2.0f;
1959 0 : return sig;
1960 : }
1961 :
1962 : Matrix<Float> &
1963 0 : VisibilityIteratorReadImpl::sigmaMat (Matrix<Float> & sigmat) const
1964 : {
1965 0 : sigmat.resize (nPol_p, curNumRow_p);
1966 0 : getCol (columns_p.sigma_p, sigmat);
1967 0 : return sigmat;
1968 : }
1969 :
1970 : Vector<Float> &
1971 200 : VisibilityIteratorReadImpl::weight (Vector<Float> & wt) const
1972 : {
1973 : // Take average of parallel hand polarizations for now.
1974 : // Later convert weight () to return full polarization dependence
1975 200 : Matrix<Float> polWeight;
1976 200 : getCol (columns_p.weight_p, polWeight);
1977 200 : wt.resize (polWeight.ncolumn ());
1978 200 : wt = polWeight.row (0);
1979 200 : wt += polWeight.row (nPol_p - 1);
1980 200 : wt /= 2.0f;
1981 400 : return wt;
1982 : }
1983 :
1984 : Matrix<Float> &
1985 0 : VisibilityIteratorReadImpl::weightMat (Matrix<Float> & wtmat) const
1986 : {
1987 0 : wtmat.resize (nPol_p, curNumRow_p);
1988 0 : getCol (columns_p.weight_p, wtmat);
1989 0 : return wtmat;
1990 : }
1991 :
1992 :
1993 : Bool
1994 212 : VisibilityIteratorReadImpl::existsWeightSpectrum () const
1995 : {
1996 212 : if (msIter_p.newMS ()) { // Cache to avoid testing unnecessarily.
1997 : try {
1998 212 : cache_p.msHasWtSp_p = columns_p.weightSpectrum_p.hasContent ();
1999 : // Comparing columns_p.weightSpectrum_p.shape (0) to
2000 : // IPosition (2, nPol_p, channelGroupSize ()) is too strict
2001 : // when channel averaging might have changed
2002 : // channelGroupSize () or weightSpectrum () out of sync. Unfortunately the
2003 : // right answer might not get cached soon enough.
2004 : //
2005 : // columns_p.weightSpectrum_p.shape (0).isEqual (IPosition (2, nPol_p,
2006 : // channelGroupSize ())));
2007 : // if (!msHasWtSp_p){
2008 : // cerr << "columns_p.weightSpectrum_p.shape (0): " << columns_p.weightSpectrum_p.shape (0) << endl;
2009 : // cerr << "(nPol_p, channelGroupSize ()): " << nPol_p
2010 : // << ", " << channelGroupSize () << endl;
2011 : // }
2012 0 : } catch (AipsError x) {
2013 0 : cache_p.msHasWtSp_p = false;
2014 : }
2015 : }
2016 212 : return cache_p.msHasWtSp_p;
2017 : }
2018 :
2019 : Cube<Float> &
2020 200 : VisibilityIteratorReadImpl::weightSpectrum (Cube<Float> & wtsp) const
2021 : {
2022 200 : if (existsWeightSpectrum ()) {
2023 0 : if (useSlicer_p) {
2024 0 : getCol (columns_p.weightSpectrum_p, slicer_p, wtsp, true);
2025 : } else {
2026 0 : getCol (columns_p.weightSpectrum_p, wtsp, true);
2027 : }
2028 : } else {
2029 200 : wtsp.resize (0, 0, 0);
2030 : }
2031 200 : return wtsp;
2032 : }
2033 :
2034 : const VisImagingWeight &
2035 200 : VisibilityIteratorReadImpl::getImagingWeightGenerator () const
2036 : {
2037 200 : return imwgt_p;
2038 : }
2039 :
2040 :
2041 : //Matrix<Float> &
2042 : //VisibilityIteratorReadImpl::imagingWeight (Matrix<Float> & wt) const
2043 : //{
2044 : // if (imwgt_p.getType () == "none") {
2045 : // throw (AipsError ("Programmer Error... imaging weights not set"));
2046 : // }
2047 : // Vector<Float> weightvec;
2048 : // weight (weightvec);
2049 : // Matrix<Bool> flagmat;
2050 : // flag (flagmat);
2051 : // wt.resize (flagmat.shape ());
2052 : // if (imwgt_p.getType () == "uniform") {
2053 : // Vector<Double> fvec;
2054 : // frequency (fvec);
2055 : // Matrix<Double> uvwmat;
2056 : // uvwMat (uvwmat);
2057 : // imwgt_p.weightUniform (wt, flagmat, uvwmat, fvec, weightvec, msId (), fieldId ());
2058 : // if (imwgt_p.doFilter ()) {
2059 : // imwgt_p.filter (wt, flagmat, uvwmat, fvec, weightvec);
2060 : // }
2061 : // } else if (imwgt_p.getType () == "radial") {
2062 : // Vector<Double> fvec;
2063 : // frequency (fvec);
2064 : // Matrix<Double> uvwmat;
2065 : // uvwMat (uvwmat);
2066 : // imwgt_p.weightRadial (wt, flagmat, uvwmat, fvec, weightvec);
2067 : // if (imwgt_p.doFilter ()) {
2068 : // imwgt_p.filter (wt, flagmat, uvwmat, fvec, weightvec);
2069 : // }
2070 : // } else {
2071 : // imwgt_p.weightNatural (wt, flagmat, weightvec);
2072 : // if (imwgt_p.doFilter ()) {
2073 : // Matrix<Double> uvwmat;
2074 : // uvwMat (uvwmat);
2075 : // Vector<Double> fvec;
2076 : // frequency (fvec);
2077 : // imwgt_p.filter (wt, flagmat, uvwmat, fvec, weightvec);
2078 : //
2079 : // }
2080 : // }
2081 : //
2082 : // return wt;
2083 : //}
2084 :
2085 : Int
2086 0 : VisibilityIteratorReadImpl::nSubInterval () const
2087 : {
2088 : // Return the number of sub-intervals in the current chunk,
2089 : // i.e. the number of unique time stamps
2090 : //
2091 : // Find all unique times in time_p
2092 0 : Int retval = 0;
2093 0 : uInt nTimes = time_p.nelements ();
2094 0 : if (nTimes > 0) {
2095 :
2096 0 : Vector<Double> times (time_p); /* Do not change time_p, make a copy */
2097 : Bool deleteIt;
2098 0 : Double * tp = times.getStorage (deleteIt);
2099 :
2100 0 : std::sort (tp, tp + nTimes);
2101 :
2102 : /* Count unique times */
2103 0 : retval = 1;
2104 0 : for (unsigned i = 0; i < nTimes - 1; i++) {
2105 0 : if (tp[i] < tp[i + 1]) {
2106 0 : retval += 1;
2107 : }
2108 : }
2109 : }
2110 0 : return retval;
2111 : }
2112 :
2113 : VisibilityIteratorReadImpl &
2114 0 : VisibilityIteratorReadImpl::selectVelocity (Int /*nChan*/,
2115 : const MVRadialVelocity & /*vStart*/,
2116 : const MVRadialVelocity & /*vInc*/,
2117 : MRadialVelocity::Types /*rvType*/,
2118 : MDoppler::Types /*dType*/,
2119 : Bool /*precise*/)
2120 : {
2121 0 : ThrowIf (true, "Method not implemented");
2122 :
2123 : // if (!initialized_p) {
2124 : // // initialize the base iterator only (avoid recursive call to originChunks)
2125 : // if (!msIterAtOrigin_p) {
2126 : // msIter_p.origin ();
2127 : // msIterAtOrigin_p = true;
2128 : // stateOk_p = false;
2129 : // }
2130 : // }
2131 : // velSelection_p = true;
2132 : // nVelChan_p = nChan;
2133 : // vstart_p = vStart;
2134 : // vinc_p = vInc;
2135 : // msd_p.setVelocityFrame (rvType);
2136 : // vDef_p = dType;
2137 : // cFromBETA_p.set (MDoppler (MVDoppler (Quantity (0., "m/s")),
2138 : // MDoppler::BETA), vDef_p);
2139 : // vPrecise_p = precise;
2140 : // if (precise) {
2141 : // // set up conversion engine for full conversion
2142 : // }
2143 : // // have to reset the iterator so all caches get filled
2144 : // originChunks ();
2145 0 : return *this;
2146 : }
2147 :
2148 :
2149 : VisibilityIteratorReadImpl &
2150 1 : VisibilityIteratorReadImpl::selectChannel (Int nGroup, Int start, Int width,
2151 : Int increment, Int spectralWindow)
2152 : {
2153 :
2154 1 : if (!initialized_p) {
2155 : // initialize the base iterator only (avoid recursive call to originChunks)
2156 0 : if (!msIterAtOrigin_p) {
2157 0 : msIter_p.origin ();
2158 0 : msIterAtOrigin_p = true;
2159 0 : stateOk_p = false;
2160 : }
2161 : }
2162 1 : Int spw = spectralWindow;
2163 1 : if (spw < 0) {
2164 0 : spw = msIter_p.spectralWindowId ();
2165 : }
2166 1 : Int n = channels_p.nGroups_p.nelements ();
2167 1 : if (n == 0) {
2168 0 : msChannels_p.spw_p.resize (1, true, false);
2169 0 : msChannels_p.spw_p[0].resize (1);
2170 0 : msChannels_p.spw_p[0][0] = spw;
2171 0 : msChannels_p.nGroups_p.resize (1, true, false);
2172 0 : msChannels_p.nGroups_p[0].resize (1);
2173 0 : msChannels_p.nGroups_p[0][0] = nGroup;
2174 0 : msChannels_p.start_p.resize (1, true, false);
2175 0 : msChannels_p.start_p[0].resize (1);
2176 0 : msChannels_p.start_p[0][0] = start;
2177 0 : msChannels_p.width_p.resize (1, true, false);
2178 0 : msChannels_p.width_p[0].resize (1);
2179 0 : msChannels_p.width_p[0][0] = width;
2180 0 : msChannels_p.inc_p.resize (1, true, false);
2181 0 : msChannels_p.inc_p[0].resize (1);
2182 0 : msChannels_p.inc_p[0][0] = increment;
2183 0 : msCounter_p = 0;
2184 :
2185 : } else {
2186 1 : Bool hasSpw = false;
2187 1 : Int spwIndex = -1;
2188 1 : for (uInt k = 0; k < msChannels_p.spw_p[0].nelements (); ++k) {
2189 1 : if (spw == msChannels_p.spw_p[0][k]) {
2190 1 : hasSpw = true;
2191 1 : spwIndex = k;
2192 1 : break;
2193 : }
2194 : }
2195 1 : if (!hasSpw) {
2196 0 : Int nspw = msChannels_p.spw_p[0].nelements () + 1;
2197 0 : msChannels_p.spw_p[0].resize (nspw, true);
2198 0 : msChannels_p.spw_p[0][nspw - 1] = spw;
2199 0 : msChannels_p.nGroups_p[0].resize (nspw, true);
2200 0 : msChannels_p.nGroups_p[0][nspw - 1] = nGroup;
2201 0 : msChannels_p.start_p[0].resize (nspw, true);
2202 0 : msChannels_p.start_p[0][nspw - 1] = start;
2203 0 : msChannels_p.width_p[0].resize (nspw, true);
2204 0 : msChannels_p.width_p[0][nspw - 1] = width;
2205 0 : msChannels_p.inc_p[0].resize (nspw, true);
2206 0 : msChannels_p.inc_p[0][nspw - 1] = increment;
2207 : } else {
2208 1 : msChannels_p.spw_p[0][spwIndex] = spw;
2209 1 : msChannels_p.nGroups_p[0][spwIndex] = nGroup;
2210 1 : msChannels_p.start_p[0][spwIndex] = start;
2211 1 : msChannels_p.width_p[0][spwIndex] = width;
2212 1 : msChannels_p.inc_p[0][spwIndex] = increment;
2213 : }
2214 :
2215 :
2216 : }
2217 1 : if (spw >= n) {
2218 : // we need to resize the blocks
2219 0 : Int newn = max (2, max (2 * n, spw + 1));
2220 0 : channels_p.nGroups_p.resize (newn);
2221 0 : channels_p.start_p.resize (newn);
2222 0 : channels_p.width_p.resize (newn);
2223 0 : channels_p.inc_p.resize (newn);
2224 0 : for (Int i = n; i < newn; i++) {
2225 0 : channels_p.nGroups_p[i] = 0;
2226 : }
2227 : }
2228 1 : channels_p.start_p[spw] = start;
2229 1 : channels_p.width_p[spw] = width;
2230 :
2231 1 : channels_p.inc_p[spw] = increment;
2232 1 : channels_p.nGroups_p[spw] = nGroup;
2233 : // have to reset the iterator so all caches get filled & slicer sizes
2234 : // get updated
2235 : // originChunks ();
2236 : //
2237 : // if (msIterAtOrigin_p){
2238 : // if (!isInSelectedSPW (msIter_p.spectralWindowId ())){
2239 : // while ((!isInSelectedSPW (msIter_p.spectralWindowId ()))
2240 : // && (msIter_p.more ()))
2241 : // msIter_p++;
2242 : // stateOk_p=false;
2243 : // setState ();
2244 : // }
2245 : // }
2246 :
2247 : //leave the state where msiter is pointing
2248 1 : channelGroupSize_p = channels_p.width_p[msIter_p.spectralWindowId ()];
2249 1 : curNGroups_p = channels_p.nGroups_p[msIter_p.spectralWindowId ()];
2250 :
2251 1 : return *this;
2252 : }
2253 :
2254 : VisibilityIteratorReadImpl &
2255 2 : VisibilityIteratorReadImpl::selectChannel (const Block<Vector<Int> > & blockNGroup,
2256 : const Block<Vector<Int> > & blockStart,
2257 : const Block<Vector<Int> > & blockWidth,
2258 : const Block<Vector<Int> > & blockIncr,
2259 : const Block<Vector<Int> > & blockSpw)
2260 : {
2261 : /*
2262 : No longer needed
2263 : if (!isMultiMS_p){
2264 : //Programmer error ...so should not reach here
2265 : cout << "Cannot use this function if Visiter was not constructed with multi-ms"
2266 : << endl;
2267 : }
2268 : */
2269 :
2270 2 : msChannels_p.nGroups_p.resize (0, true, false);
2271 2 : msChannels_p.nGroups_p = blockNGroup;
2272 2 : msChannels_p.start_p.resize (0, true, false);
2273 2 : msChannels_p.start_p = blockStart;
2274 2 : msChannels_p.width_p.resize (0, true, false);
2275 2 : msChannels_p.width_p = blockWidth;
2276 2 : msChannels_p.inc_p.resize (0, true, false);
2277 2 : msChannels_p.inc_p = blockIncr;
2278 2 : msChannels_p.spw_p.resize (0, true, false);
2279 2 : msChannels_p.spw_p = blockSpw;
2280 :
2281 2 : if (!initialized_p) {
2282 : // initialize the base iterator only (avoid recursive call to originChunks)
2283 2 : if (!msIterAtOrigin_p) {
2284 2 : msIter_p.origin ();
2285 2 : msIterAtOrigin_p = true;
2286 2 : stateOk_p = false;
2287 : }
2288 : }
2289 :
2290 2 : channels_p.nGroups_p.resize (0);
2291 2 : msCounter_p = 0;
2292 :
2293 2 : doChannelSelection ();
2294 : // have to reset the iterator so all caches get filled & slicer sizes
2295 : // get updated
2296 :
2297 2 : if (msIterAtOrigin_p) {
2298 2 : if (!isInSelectedSPW (msIter_p.spectralWindowId ())) {
2299 0 : while ((!isInSelectedSPW (msIter_p.spectralWindowId ()))
2300 0 : && (msIter_p.more ())) {
2301 0 : msIter_p++;
2302 : }
2303 0 : stateOk_p = false;
2304 : }
2305 :
2306 : }
2307 :
2308 2 : originChunks ();
2309 2 : return *this;
2310 : }
2311 :
2312 :
2313 : void
2314 0 : VisibilityIteratorReadImpl::getChannelSelection (Block< Vector<Int> > & blockNGroup,
2315 : Block< Vector<Int> > & blockStart,
2316 : Block< Vector<Int> > & blockWidth,
2317 : Block< Vector<Int> > & blockIncr,
2318 : Block< Vector<Int> > & blockSpw)
2319 : {
2320 :
2321 0 : blockNGroup.resize (0, true, false);
2322 0 : blockNGroup = msChannels_p.nGroups_p;
2323 0 : blockStart.resize (0, true, false);
2324 0 : blockStart = msChannels_p.start_p;
2325 0 : blockWidth.resize (0, true, false);
2326 0 : blockWidth = msChannels_p.width_p;
2327 0 : blockIncr.resize (0, true, false);
2328 0 : blockIncr = msChannels_p.inc_p;
2329 0 : blockSpw.resize (0, true, false);
2330 0 : blockSpw = msChannels_p.spw_p;
2331 0 : }
2332 : void
2333 14 : VisibilityIteratorReadImpl::doChannelSelection ()
2334 : {
2335 28 : for (uInt k = 0; k < msChannels_p.spw_p[msCounter_p].nelements (); ++k) {
2336 14 : Int spw = msChannels_p.spw_p[msCounter_p][k];
2337 14 : if (spw < 0) {
2338 0 : spw = msIter_p.spectralWindowId ();
2339 : }
2340 14 : Int n = channels_p.nGroups_p.nelements ();
2341 14 : if (spw >= n) {
2342 : // we need to resize the blocks
2343 14 : Int newn = max (2, max (2 * n, spw + 1));
2344 14 : channels_p.nGroups_p.resize (newn, true, true);
2345 14 : channels_p.start_p.resize (newn, true, true);
2346 14 : channels_p.width_p.resize (newn, true, true);
2347 14 : channels_p.inc_p.resize (newn, true, true);
2348 42 : for (Int i = n; i < newn; i++) {
2349 28 : channels_p.nGroups_p[i] = 0;
2350 : }
2351 : }
2352 :
2353 14 : channels_p.start_p[spw] = msChannels_p.start_p[msCounter_p][k];
2354 14 : channels_p.width_p[spw] = msChannels_p.width_p[msCounter_p][k];
2355 14 : channelGroupSize_p = msChannels_p.width_p[msCounter_p][k];
2356 14 : channels_p.inc_p[spw] = msChannels_p.inc_p[msCounter_p][k];
2357 14 : channels_p.nGroups_p[spw] = msChannels_p.nGroups_p[msCounter_p][k];
2358 14 : curNGroups_p = msChannels_p.nGroups_p[msCounter_p][k];
2359 :
2360 : }
2361 14 : Int spw = msIter_p.spectralWindowId ();
2362 14 : Int spIndex = -1;
2363 14 : for (uInt k = 0; k < msChannels_p.spw_p[msCounter_p].nelements (); ++k) {
2364 14 : if (spw == msChannels_p.spw_p[msCounter_p][k]) {
2365 14 : spIndex = k;
2366 14 : break;
2367 : }
2368 : }
2369 :
2370 :
2371 14 : if (spIndex < 0) {
2372 0 : spIndex = 0;
2373 : }
2374 : //leave this at the stage where msiter is pointing
2375 14 : channelGroupSize_p = msChannels_p.width_p[msCounter_p][spIndex];
2376 14 : curNGroups_p = msChannels_p.nGroups_p[msCounter_p][spIndex];
2377 :
2378 :
2379 :
2380 14 : }
2381 :
2382 : void
2383 0 : VisibilityIteratorReadImpl::slicesToMatrices (Vector<Matrix<Int> > & matv,
2384 : const Vector<Vector<Slice> > & slicesv,
2385 : const Vector<Int> & widthsv) const
2386 : {
2387 0 : uInt nspw = slicesv.nelements ();
2388 :
2389 0 : matv.resize (nspw);
2390 0 : uInt selspw = 0;
2391 0 : for (uInt spw = 0; spw < nspw; ++spw) {
2392 0 : uInt nSlices = slicesv[spw].nelements ();
2393 :
2394 : // Figure out how big to make matv[spw].
2395 0 : uInt totOutChan = 0;
2396 :
2397 0 : Int width = (nSlices > 0) ? widthsv[selspw] : 1;
2398 0 : if (width < 1) {
2399 0 : throw (AipsError ("Cannot channel average with width < 1"));
2400 : }
2401 :
2402 0 : for (uInt slicenum = 0; slicenum < nSlices; ++slicenum) {
2403 0 : const Slice & sl = slicesv[spw][slicenum];
2404 0 : Int firstchan = sl.start ();
2405 0 : Int lastchan = sl.all () ? firstchan + channels_p.width_p[spw] - 1 : sl.end ();
2406 0 : Int inc = sl.all () ? 1 : sl.inc ();
2407 :
2408 : // Even if negative increments are desirable, the for loop below has a <.
2409 0 : if (inc < 1) {
2410 0 : throw (AipsError ("The channel increment must be >= 1"));
2411 : }
2412 :
2413 : // This formula is very dependent on integer division. Don't rearrange it.
2414 0 : totOutChan += 1 + ((lastchan - firstchan) / inc) / (1 + (width - 1) / inc);
2415 : }
2416 0 : matv[spw].resize (totOutChan, 4);
2417 :
2418 : // Index of input channel in SELECTED list, i.e.
2419 : // mschan = vi.chanIds (chanids, spw)[selchanind].
2420 0 : uInt selchanind = 0;
2421 :
2422 : // Fill matv with channel boundaries.
2423 0 : uInt outChan = 0;
2424 0 : for (uInt slicenum = 0; slicenum < nSlices; ++slicenum) {
2425 0 : const Slice & sl = slicesv[spw][slicenum];
2426 0 : Int firstchan = sl.start ();
2427 0 : Int lastchan = sl.all () ? firstchan + channels_p.width_p[spw] - 1 : sl.end ();
2428 0 : Int inc = sl.all () ? 1 : sl.inc (); // Default to no skipping
2429 :
2430 : // Again, these depend on integer division. Don't rearrange them.
2431 0 : Int selspan = 1 + (width - 1) / inc;
2432 0 : Int span = inc * selspan;
2433 :
2434 0 : for (Int mschan = firstchan; mschan <= lastchan; mschan += span) {
2435 : // The start and end in MS channel #s.
2436 0 : matv[spw](outChan, 0) = mschan;
2437 0 : matv[spw](outChan, 1) = mschan + width - 1;
2438 :
2439 : // The start and end in selected reckoning.
2440 0 : matv[spw](outChan, 2) = selchanind;
2441 0 : selchanind += selspan;
2442 0 : matv[spw](outChan, 3) = selchanind - 1;
2443 0 : ++outChan;
2444 : }
2445 : }
2446 0 : if (nSlices > 0) { // spw was selected
2447 0 : ++selspw;
2448 : }
2449 : }
2450 0 : }
2451 :
2452 :
2453 0 : void VisibilityIteratorReadImpl::getFreqInSpwRange(Double& freqStart, Double& freqEnd, MFrequency::Types freqframe) const {
2454 0 : Int nMS = msIter_p.numMS ();
2455 0 : freqStart=C::dbl_max;
2456 0 : freqEnd = 0.0;
2457 :
2458 0 : for (Int msId=0; msId < nMS; ++msId){
2459 :
2460 0 : Vector<Int> spws = msChannels_p.spw_p[msId];
2461 0 : Vector<Int> starts = msChannels_p.start_p[msId];
2462 0 : Vector<Int> nchan;
2463 : // careful here as we don't want to change width_p as it is multiplied by incr below
2464 0 : nchan= msChannels_p.width_p[msId];
2465 0 : Vector<Int> incr = msChannels_p.inc_p[msId];
2466 0 : nchan=nchan*incr;
2467 0 : Vector<uInt> uniqIndx;
2468 0 : Vector<Int> fldId;
2469 0 : ScalarColumn<Int> (msIter_p.ms (msId), MS::columnName (MS::FIELD_ID)).getColumn (fldId);
2470 0 : uInt nFields = GenSort<Int>::sort (fldId, Sort::Ascending, Sort::QuickSort | Sort::NoDuplicates);
2471 0 : for (uInt indx=0; indx< nFields; ++indx){
2472 0 : Int fieldid=fldId(indx);
2473 : Double tmpFreqStart;
2474 : Double tmpFreqEnd;
2475 0 : MSUtil::getFreqRangeInSpw(tmpFreqStart, tmpFreqEnd, spws, starts, nchan, msIter_p.ms(msId), freqframe, fieldid);
2476 0 : if(freqStart > tmpFreqStart) freqStart=tmpFreqStart;
2477 0 : if(freqEnd < tmpFreqEnd) freqEnd=tmpFreqEnd;
2478 : }
2479 : }
2480 :
2481 0 : }
2482 :
2483 :
2484 :
2485 : void
2486 0 : VisibilityIteratorReadImpl::getSpwInFreqRange (Block<Vector<Int> > & spw,
2487 : Block<Vector<Int> > & start,
2488 : Block<Vector<Int> > & nchan,
2489 : Double freqStart,
2490 : Double freqEnd,
2491 : Double freqStep,
2492 : MFrequency::Types freqframe) const
2493 : {
2494 : // This functionality was relocated from MSIter in order to support this operation
2495 : // within the VI to make the VisibilityIteratorReadImplAsync implementation feasible.
2496 :
2497 0 : Int nMS = msIter_p.numMS ();
2498 :
2499 0 : spw.resize (nMS, true, false);
2500 0 : start.resize (nMS, true, false);
2501 0 : nchan.resize (nMS, true, false);
2502 :
2503 0 : for (Int k = 0; k < nMS; ++k) {
2504 0 : Vector<Double> t;
2505 0 : ScalarColumn<Double> (msIter_p.ms (k), MS::columnName (MS::TIME)).getColumn (t);
2506 0 : Vector<Int> ddId;
2507 0 : Vector<Int> fldId;
2508 0 : ScalarColumn<Int> (msIter_p.ms (k), MS::columnName (MS::DATA_DESC_ID)).getColumn (ddId);
2509 0 : ScalarColumn<Int> (msIter_p.ms (k), MS::columnName (MS::FIELD_ID)).getColumn (fldId);
2510 0 : MSFieldColumns fieldCol (msIter_p.ms (k).field ());
2511 0 : MSDataDescColumns ddCol (msIter_p.ms (k).dataDescription ());
2512 0 : MSSpWindowColumns spwCol (msIter_p.ms (k).spectralWindow ());
2513 0 : ROScalarMeasColumn<MEpoch> timeCol (msIter_p.ms (k), MS::columnName (MS::TIME));
2514 0 : Vector<uInt> uniqIndx;
2515 0 : uInt nTimes = GenSortIndirect<Double>::sort (uniqIndx, t, Sort::Ascending, Sort::QuickSort | Sort::NoDuplicates);
2516 : //now need to do the conversion to data frame from requested frame
2517 : //Get the epoch mesasures of the first row
2518 0 : MEpoch ep;
2519 0 : timeCol.get (0, ep);
2520 0 : MPosition obsPos = msIter_p.telescopePosition ();
2521 0 : Int oldDD = ddId[0];
2522 0 : Int oldFLD = fldId[0];
2523 : //For now we will assume that the field is not moving very far from polynome 0
2524 0 : MDirection dir = fieldCol.phaseDirMeas (fldId[0]);
2525 0 : MFrequency::Types obsMFreqType = (MFrequency::Types) (spwCol.measFreqRef ()(ddCol.spectralWindowId ()(ddId[0])));
2526 : //cout << "nTimes " << nTimes << endl;
2527 : //cout << " obsframe " << obsMFreqType << " reqFrame " << freqframe << endl;
2528 0 : MeasFrame frame (ep, obsPos, dir);
2529 : MFrequency::Convert toObs (freqframe,
2530 0 : MFrequency::Ref (obsMFreqType, frame));
2531 :
2532 0 : Double freqEndMax = freqEnd;
2533 0 : Double freqStartMin = freqStart;
2534 0 : if (freqframe != obsMFreqType) {
2535 0 : freqEndMax = 0.0;
2536 0 : freqStartMin = C::dbl_max;
2537 : }
2538 :
2539 0 : for (uInt j = 0; j < nTimes; ++j) {
2540 0 : timeCol.get (uniqIndx[j], ep);
2541 0 : if (oldDD != ddId[uniqIndx[j]]) {
2542 0 : oldDD = ddId[uniqIndx[j]];
2543 0 : if (spwCol.measFreqRef ()(ddCol.spectralWindowId ()(ddId[uniqIndx[j]])) != obsMFreqType) {
2544 0 : obsMFreqType = (MFrequency::Types) (spwCol.measFreqRef ()(ddCol.spectralWindowId ()(ddId[uniqIndx[j]])));
2545 0 : toObs.setOut (MFrequency::Ref (obsMFreqType, frame));
2546 : }
2547 : }
2548 0 : if (obsMFreqType != freqframe) {
2549 0 : frame.resetEpoch (ep);
2550 0 : if (oldFLD != fldId[uniqIndx[j]]) {
2551 0 : oldFLD = fldId[uniqIndx[j]];
2552 0 : frame.resetDirection (fieldCol.phaseDirMeas (fldId[uniqIndx[j]]));
2553 : }
2554 0 : Double freqTmp = toObs (Quantity (freqStart, "Hz")).get ("Hz").getValue ();
2555 0 : freqStartMin = (freqStartMin > freqTmp) ? freqTmp : freqStartMin;
2556 0 : freqTmp = toObs (Quantity (freqEnd, "Hz")).get ("Hz").getValue ();
2557 0 : freqEndMax = (freqEndMax < freqTmp) ? freqTmp : freqEndMax;
2558 : }
2559 : }
2560 :
2561 : //cout << "freqStartMin " << freqStartMin << " freqEndMax " << freqEndMax << endl;
2562 0 : MSSpwIndex spwIn (msIter_p.ms (k).spectralWindow ());
2563 :
2564 0 : spwIn.matchFrequencyRange (freqStartMin - 0.5 * freqStep, freqEndMax + 0.5 * freqStep, spw[k], start[k], nchan[k]);
2565 : }
2566 0 : }
2567 :
2568 : vector<MeasurementSet>
2569 0 : VisibilityIteratorReadImpl::getMeasurementSets () const
2570 : {
2571 0 : return measurementSets_p;
2572 : }
2573 :
2574 :
2575 : void
2576 1 : VisibilityIteratorReadImpl::allSelectedSpectralWindows (Vector<Int> & spws, Vector<Int> & nvischan)
2577 : {
2578 :
2579 1 : spws.resize ();
2580 1 : spws = msChannels_p.spw_p[msId ()];
2581 1 : nvischan.resize ();
2582 1 : nvischan.resize (max (spws) + 1);
2583 1 : nvischan.set (-1);
2584 2 : for (uInt k = 0; k < spws.nelements (); ++k) {
2585 1 : nvischan[spws[k]] = channels_p.width_p[spws[k]];
2586 : }
2587 1 : }
2588 :
2589 : Vector<Double> &
2590 0 : VisibilityIteratorReadImpl::lsrFrequency (Vector<Double> & freq) const
2591 : {
2592 0 : if (velocity_p.velSelection_p) {
2593 0 : freq.resize (velocity_p.nVelChan_p);
2594 0 : freq = velocity_p.lsrFreq_p;
2595 : } else {
2596 : // if there is no vel selection, we just return the observing freqs
2597 0 : frequency (freq);
2598 : }
2599 0 : return freq;
2600 : }
2601 :
2602 :
2603 : void
2604 1 : VisibilityIteratorReadImpl::lsrFrequency (const Int & spw, Vector<Double> & freq,
2605 : Bool & convert, const Bool ignoreconv)
2606 : {
2607 : // This method is not good for conversion between frames which are extremely
2608 : // time dependent over the course of the observation e.g topo to lsr unless
2609 : // the epoch is in the actual buffer
2610 :
2611 1 : if (velocity_p.velSelection_p) {
2612 0 : getTopoFreqs ();
2613 0 : lsrFrequency (freq);
2614 0 : return;
2615 : }
2616 :
2617 1 : if (! cache_p.freqCacheOK_p) {
2618 0 : frequency (freq);
2619 : }
2620 :
2621 : //MFrequency::Types obsMFreqType=(MFrequency::Types)(msIter_p.msColumns ().spectralWindow ().measFreqRef ()(spw));
2622 :
2623 : //chanFreq=msIter_p.msColumns ().spectralWindow ().chanFreq ()(spw);
2624 :
2625 1 : const ArrayColumn <Double> & chanFreqs = msIter_p.msColumns ().spectralWindow ().chanFreq ();
2626 1 : const ScalarColumn<Int> & obsMFreqTypes = msIter_p.msColumns ().spectralWindow ().measFreqRef ();
2627 2 : MEpoch ep;
2628 1 : ROScalarMeasColumn<MEpoch>(msIter_p.table (), MS::columnName (MS::TIME)).get (curStartRow_p, ep); // Setting epoch to iteration's first one
2629 2 : MPosition obsPos = msIter_p.telescopePosition ();
2630 2 : MDirection dir = msIter_p.phaseCenter ();
2631 :
2632 2 : lsrFrequency (spw, freq, convert, channels_p.start_p, channels_p.width_p, channels_p.inc_p,
2633 1 : channels_p.nGroups_p, chanFreqs, obsMFreqTypes, ep, obsPos, dir, ignoreconv);
2634 :
2635 : }
2636 :
2637 : void
2638 1 : VisibilityIteratorReadImpl::lsrFrequency (const Int & spw,
2639 : Vector<Double> & freq,
2640 : Bool & convert,
2641 : const Block<Int> & chanStart,
2642 : const Block<Int> & chanWidth,
2643 : const Block<Int> & chanInc,
2644 : const Block<Int> & numChanGroup,
2645 : const ArrayColumn <Double> & chanFreqs,
2646 : const ScalarColumn<Int> & obsMFreqTypes,
2647 : const MEpoch & ep,
2648 : const MPosition & obsPos,
2649 : const MDirection & dir, const Bool ignoreconv)
2650 : {
2651 :
2652 2 : Vector<Double> chanFreq (0);
2653 1 : chanFreq = chanFreqs (spw);
2654 :
2655 : //chanFreq=msIter_p.msColumns ().spectralWindow ().chanFreq ()(spw);
2656 : // Int start=channels_p.start_p[spw]-msIter_p.startChan ();
2657 : //Assuming that the spectral windows selected is not a reference ms from
2658 : //visset ...as this will have a start chan offseted may be.
2659 :
2660 1 : Int start = chanStart[spw];
2661 1 : freq.resize (chanWidth[spw]);
2662 1 : MFrequency::Types obsMFreqType = (MFrequency::Types) (obsMFreqTypes (spw));
2663 2 : MeasFrame frame (ep, obsPos, dir);
2664 : MFrequency::Convert tolsr (obsMFreqType,
2665 2 : MFrequency::Ref (MFrequency::LSRK, frame));
2666 :
2667 : // if (obsMFreqType != MFrequency::LSRK){
2668 : // convert=true;
2669 : // }
2670 :
2671 1 : convert = obsMFreqType != MFrequency::LSRK; // make this parameter write-only
2672 : // user requested no conversion
2673 1 : if(ignoreconv) convert=false;
2674 :
2675 101 : for (Int i = 0; i < chanWidth[spw]; i++) {
2676 100 : Int inc = chanInc[spw] <= 0 ? 1 : chanInc[spw] ;
2677 100 : if (convert) {
2678 0 : freq[i] = tolsr (chanFreq (start +
2679 0 : (numChanGroup[spw] - 1) * chanWidth[spw] + i * inc)).
2680 0 : getValue ().getValue ();
2681 : } else {
2682 100 : freq[i] = chanFreq (start +
2683 100 : (numChanGroup[spw] - 1) * chanWidth[spw] + i * inc);
2684 : }
2685 : }
2686 1 : }
2687 :
2688 : void
2689 0 : VisibilityIteratorReadImpl::getLsrInfo (Block<Int> & channelGroupNumber,
2690 : Block<Int> & channelIncrement,
2691 : Block<Int> & channelStart,
2692 : Block<Int> & channelWidth,
2693 : MPosition & observatoryPositon,
2694 : MDirection & phaseCenter,
2695 : Bool & velocitySelection) const
2696 : {
2697 0 : channelStart = channels_p.start_p;
2698 0 : channelWidth = channels_p.width_p;
2699 0 : channelIncrement = channels_p.inc_p;
2700 0 : channelGroupNumber = channels_p.nGroups_p;
2701 0 : observatoryPositon = msIter_p.telescopePosition ();
2702 0 : phaseCenter = msIter_p.phaseCenter ();
2703 0 : velocitySelection = velocity_p.velSelection_p;
2704 0 : }
2705 :
2706 :
2707 : void
2708 2 : VisibilityIteratorReadImpl::attachVisBuffer (VisBuffer & vb)
2709 : {
2710 2 : vbStack_p.push (& vb);
2711 2 : vb.invalidate ();
2712 2 : }
2713 :
2714 : VisBuffer *
2715 0 : VisibilityIteratorReadImpl::getVisBuffer ()
2716 : {
2717 0 : VisBuffer * result = NULL;
2718 :
2719 0 : if (! vbStack_p.empty ()) {
2720 0 : result = vbStack_p.top ();
2721 : }
2722 :
2723 0 : return result;
2724 : }
2725 :
2726 : void
2727 1 : VisibilityIteratorReadImpl::detachVisBuffer (VisBuffer & vb)
2728 : {
2729 1 : if (!vbStack_p.empty ()) {
2730 1 : if (vbStack_p.top () == & vb) {
2731 1 : vbStack_p.pop ();
2732 1 : if (!vbStack_p.empty ()) {
2733 0 : vbStack_p.top ()->invalidate ();
2734 : }
2735 : } else {
2736 : throw (AipsError ("ROVisIter::detachVisBuffer - attempt to detach "
2737 0 : "buffer that is not the last one attached"));
2738 : }
2739 : }
2740 1 : }
2741 :
2742 : Int
2743 0 : VisibilityIteratorReadImpl::numberAnt ()
2744 : {
2745 0 : return msColumns ().antenna ().nrow (); // for single (sub)array only..
2746 : }
2747 :
2748 : Int
2749 0 : VisibilityIteratorReadImpl::numberSpw ()
2750 : {
2751 0 : return msColumns ().spectralWindow ().nrow ();
2752 : }
2753 :
2754 : Int
2755 0 : VisibilityIteratorReadImpl::numberDDId ()
2756 : {
2757 0 : return msColumns ().dataDescription ().nrow ();
2758 : }
2759 :
2760 : Int
2761 0 : VisibilityIteratorReadImpl::numberPol ()
2762 : {
2763 0 : return msColumns ().polarization ().nrow ();
2764 : }
2765 :
2766 : Int
2767 0 : VisibilityIteratorReadImpl::numberCoh ()
2768 : {
2769 0 : Int numcoh = 0;
2770 0 : for (uInt k = 0; k < uInt (msIter_p.numMS ()) ; ++k) {
2771 0 : numcoh += msIter_p.ms (k).nrow ();
2772 : }
2773 0 : return numcoh;
2774 :
2775 : }
2776 :
2777 : template<class T>
2778 : void
2779 801 : VisibilityIteratorReadImpl::getColScalar (const ScalarColumn<T> &column, Vector<T> &array, Bool resize) const
2780 : {
2781 801 : column.getColumnCells (selRows_p, array, resize);
2782 801 : return;
2783 : }
2784 :
2785 : void
2786 200 : VisibilityIteratorReadImpl::getCol (const ScalarColumn<Bool> &column, Vector<Bool> &array, Bool resize) const
2787 : {
2788 200 : getColScalar<Bool>(column, array, resize);
2789 200 : }
2790 :
2791 : void
2792 400 : VisibilityIteratorReadImpl::getCol (const ScalarColumn<Int> &column, Vector<Int> &array, Bool resize) const
2793 : {
2794 400 : getColScalar<Int>(column, array, resize);
2795 400 : }
2796 :
2797 : void
2798 201 : VisibilityIteratorReadImpl::getCol (const ScalarColumn<Double> &column, Vector<Double> &array, Bool resize) const
2799 : {
2800 201 : getColScalar<Double>(column, array, resize);
2801 201 : }
2802 :
2803 : void
2804 460 : VisibilityIteratorReadImpl::getCol (const ArrayColumn<Bool> &column, Array<Bool> &array, Bool resize) const
2805 : {
2806 460 : column.getColumnCells (selRows_p, array, resize);
2807 460 : }
2808 :
2809 200 : void VisibilityIteratorReadImpl::getCol (const ArrayColumn<Float> &column, Array<Float> &array, Bool resize) const
2810 : {
2811 200 : column.getColumnCells (selRows_p, array, resize);
2812 200 : }
2813 :
2814 : template<class T>
2815 : void
2816 200 : VisibilityIteratorReadImpl::getColArray (const ArrayColumn<T> &column, Array<T> &array, Bool resize) const
2817 : {
2818 200 : column.getColumnCells (selRows_p, array, resize);
2819 200 : return;
2820 : }
2821 :
2822 : void
2823 0 : VisibilityIteratorReadImpl::getCol (const ArrayColumn<Double> &column, Array<Double> &array, Bool resize) const
2824 : {
2825 0 : column.getColumnCells (selRows_p, array, resize);
2826 0 : }
2827 :
2828 : void
2829 266 : VisibilityIteratorReadImpl::getCol (const ArrayColumn<Complex> &column, Array<Complex> &array, Bool resize) const
2830 : {
2831 266 : column.getColumnCells (selRows_p, array, resize);
2832 266 : }
2833 :
2834 : void
2835 0 : VisibilityIteratorReadImpl::getCol (const ArrayColumn<Bool> &column, const Slicer & slicer, Array<Bool> &array, Bool resize) const
2836 : {
2837 0 : column.getColumnCells (selRows_p, slicer, array, resize);
2838 0 : }
2839 :
2840 : void
2841 0 : VisibilityIteratorReadImpl::getCol (const ArrayColumn<Float> &column, const Slicer & slicer, Array<Float> &array, Bool resize) const
2842 : {
2843 0 : column.getColumnCells (selRows_p, slicer, array, resize);
2844 0 : }
2845 :
2846 : void
2847 0 : VisibilityIteratorReadImpl::getCol (const ArrayColumn<Complex> &column, const Slicer & slicer, Array<Complex> &array, Bool resize) const
2848 : {
2849 0 : column.getColumnCells (selRows_p, slicer, array, resize);
2850 0 : }
2851 :
2852 : const Table
2853 35 : VisibilityIteratorReadImpl::attachTable () const
2854 : {
2855 35 : return msIter_p.table ();
2856 : }
2857 :
2858 : void
2859 0 : VisibilityIteratorReadImpl::slurp () const
2860 : {
2861 : /* Set the table data manager (ISM and SSM) cache size to the full column size, for
2862 : the columns ANTENNA1, ANTENNA2, FEED1, FEED2, TIME, INTERVAL, FLAG_ROW, SCAN_NUMBER and UVW
2863 : */
2864 0 : Record dmInfo (msIter_p.ms ().dataManagerInfo ());
2865 :
2866 : // cout << "nfields = " << dmInfo.nfields () << endl;
2867 : // cout << "dminfo = " << dmInfo.description () << endl;
2868 0 : RecordDesc desc = dmInfo.description ();
2869 0 : for (unsigned i = 0; i < dmInfo.nfields (); i++) {
2870 : // cout << "field " << i << " isSubRecord = " << desc.isSubRecord (i) << endl;
2871 : // cout << "field " << i << " isArray = " << desc.isArray (i) << endl;
2872 0 : if (desc.isSubRecord (i)) {
2873 :
2874 0 : Record sub = dmInfo.subRecord (i);
2875 :
2876 : // cout << "sub = " << sub << endl;
2877 0 : if (sub.fieldNumber ("NAME") >= 0 &&
2878 0 : sub.fieldNumber ("TYPE") >= 0 &&
2879 0 : sub.fieldNumber ("COLUMNS") >= 0 &&
2880 0 : sub.type (sub.fieldNumber ("NAME")) == TpString &&
2881 0 : sub.type (sub.fieldNumber ("TYPE")) == TpString &&
2882 0 : sub.type (sub.fieldNumber ("COLUMNS")) == TpArrayString) {
2883 :
2884 0 : Array<String> columns;
2885 0 : dmInfo.subRecord (i).get ("COLUMNS", columns);
2886 :
2887 0 : bool match = false;
2888 0 : for (unsigned j = 0; j < columns.nelements (); j++) {
2889 0 : String column = columns (IPosition (1, j));
2890 0 : match |= (column == MS::columnName (MS::ANTENNA1) ||
2891 0 : column == MS::columnName (MS::ANTENNA2) ||
2892 0 : column == MS::columnName (MS::FEED1) ||
2893 0 : column == MS::columnName (MS::FEED2) ||
2894 0 : column == MS::columnName (MS::TIME) ||
2895 0 : column == MS::columnName (MS::INTERVAL) ||
2896 0 : column == MS::columnName (MS::FLAG_ROW) ||
2897 0 : column == MS::columnName (MS::SCAN_NUMBER) ||
2898 0 : column == MS::columnName (MS::UVW));
2899 : }
2900 : // cout << "columns = " << columns << endl;
2901 :
2902 0 : if (match) {
2903 :
2904 0 : String dm_name;
2905 0 : dmInfo.subRecord (i).get ("NAME", dm_name);
2906 : // cout << "dm_name = " << dm_name << endl;
2907 :
2908 0 : String dm_type;
2909 0 : dmInfo.subRecord (i).get ("TYPE", dm_type);
2910 : // cout << "dm_type = " << dm_type << endl;
2911 :
2912 0 : Bool can_exceed_nr_buckets = false;
2913 0 : uInt num_buckets = msIter_p.ms ().nrow ();
2914 : // One bucket is at least one row, so this is enough
2915 :
2916 0 : if (dm_type == "IncrementalStMan") {
2917 0 : ROIncrementalStManAccessor acc (msIter_p.ms (), dm_name);
2918 0 : acc.setCacheSize (num_buckets, can_exceed_nr_buckets);
2919 0 : } else if (dm_type == "StandardStMan") {
2920 0 : ROStandardStManAccessor acc (msIter_p.ms (), dm_name);
2921 0 : acc.setCacheSize (num_buckets, can_exceed_nr_buckets);
2922 : }
2923 : /* These are the only storage managers which use the BucketCache
2924 : (and therefore are slow for random access and small cache sizes)
2925 : */
2926 : } else {
2927 0 : String dm_name;
2928 0 : dmInfo.subRecord (i).get ("NAME", dm_name);
2929 : //cout << "IGNORING...." << dm_name << endl;
2930 : }
2931 : } else {
2932 0 : cerr << "Data manager info has unexpected shape! " << sub << endl;
2933 : }
2934 : }
2935 : }
2936 0 : return;
2937 : }
2938 :
2939 2 : VisibilityIteratorWriteImpl::VisibilityIteratorWriteImpl (VisibilityIterator * vi)
2940 2 : : vi_p (vi)
2941 : {
2942 2 : useCustomTileShape_p = useCustomTileShape();
2943 2 : }
2944 :
2945 0 : VisibilityIteratorWriteImpl::VisibilityIteratorWriteImpl (const VisibilityIteratorWriteImpl & other)
2946 : {
2947 0 : operator=(other);
2948 0 : useCustomTileShape_p = useCustomTileShape();
2949 0 : }
2950 :
2951 2 : VisibilityIteratorWriteImpl::~VisibilityIteratorWriteImpl ()
2952 : {
2953 2 : }
2954 :
2955 : VisibilityIteratorWriteImpl *
2956 0 : VisibilityIteratorWriteImpl::clone (VisibilityIterator * vi) const
2957 : {
2958 : // Return a clone of this object but ensure it's attached to the proper
2959 : // VI container object.
2960 :
2961 0 : VisibilityIteratorWriteImpl * viwi = new VisibilityIteratorWriteImpl (* this);
2962 :
2963 0 : viwi->vi_p = vi;
2964 :
2965 0 : return viwi;
2966 : }
2967 :
2968 : void
2969 31 : VisibilityIteratorWriteImpl::attachColumns (const Table & t)
2970 : {
2971 31 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
2972 :
2973 31 : readImpl->attachColumns (t); // Let the read impl handle it's tables
2974 :
2975 31 : const ColumnDescSet & cds = t.tableDesc ().columnDescSet ();
2976 :
2977 31 : if (cds.isDefined (MS::columnName (MS::DATA))) {
2978 31 : columns_p.vis_p.attach (t, MS::columnName (MS::DATA));
2979 : }
2980 :
2981 31 : if (cds.isDefined (MS::columnName (MS::FLOAT_DATA))) {
2982 0 : readImpl->floatDataFound_p = true;
2983 0 : columns_p.floatVis_p.attach (t, MS::columnName (MS::FLOAT_DATA));
2984 : } else {
2985 31 : readImpl->floatDataFound_p = false;
2986 : }
2987 :
2988 31 : if (cds.isDefined ("MODEL_DATA")) {
2989 31 : columns_p.modelVis_p.attach (t, "MODEL_DATA");
2990 : }
2991 :
2992 31 : if (cds.isDefined ("CORRECTED_DATA")) {
2993 31 : columns_p.corrVis_p.attach (t, "CORRECTED_DATA");
2994 : }
2995 :
2996 31 : columns_p.weight_p.attach (t, MS::columnName (MS::WEIGHT));
2997 :
2998 31 : if (cds.isDefined ("WEIGHT_SPECTRUM")) {
2999 27 : columns_p.weightSpectrum_p.attach (t, "WEIGHT_SPECTRUM");
3000 : }
3001 :
3002 31 : columns_p.sigma_p.attach (t, MS::columnName (MS::SIGMA));
3003 :
3004 31 : columns_p.flag_p.attach (t, MS::columnName (MS::FLAG));
3005 :
3006 31 : columns_p.flagRow_p.attach (t, MS::columnName (MS::FLAG_ROW));
3007 :
3008 31 : if (cds.isDefined ("FLAG_CATEGORY")) {
3009 31 : columns_p.flagCategory_p.attach (t, MS::columnName (MS::FLAG_CATEGORY));
3010 : }
3011 31 : }
3012 :
3013 : VisibilityIteratorReadImpl *
3014 307 : VisibilityIteratorWriteImpl::getReadImpl ()
3015 : {
3016 307 : return vi_p->getReadImpl();
3017 : }
3018 :
3019 : void
3020 0 : VisibilityIteratorWriteImpl::setFlag (const Matrix<Bool> & flag)
3021 : {
3022 : // use same value for all polarizations
3023 :
3024 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3025 :
3026 0 : readImpl->cache_p.flagCube_p.resize (readImpl->nPol_p, readImpl->channelGroupSize_p,
3027 0 : readImpl->curNumRow_p);
3028 :
3029 : Bool deleteIt;
3030 0 : Bool * p = readImpl->cache_p.flagCube_p.getStorage (deleteIt);
3031 0 : const Bool * pflag = flag.getStorage (deleteIt);
3032 0 : if (Int (flag.nrow ()) != readImpl->channelGroupSize_p) {
3033 0 : throw (AipsError ("VisIter::setFlag (flag) - inconsistent number of channels"));
3034 : }
3035 :
3036 0 : for (uInt row = 0; row < readImpl->curNumRow_p; row++) {
3037 0 : for (Int chn = 0; chn < readImpl->channelGroupSize_p; chn++) {
3038 0 : for (Int pol = 0; pol < readImpl->nPol_p; pol++) {
3039 0 : *p++ = *pflag;
3040 : }
3041 0 : pflag++;
3042 : }
3043 : }
3044 :
3045 0 : if (readImpl->useSlicer_p) {
3046 0 : putCol (columns_p.flag_p, readImpl->slicer_p, readImpl->cache_p.flagCube_p);
3047 : } else {
3048 0 : putCol (columns_p.flag_p, readImpl->cache_p.flagCube_p);
3049 : }
3050 0 : }
3051 :
3052 : void
3053 0 : VisibilityIteratorWriteImpl::setFlag (const Cube<Bool> & flags)
3054 : {
3055 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3056 :
3057 0 : if (readImpl->useSlicer_p) {
3058 0 : putCol (columns_p.flag_p, readImpl->slicer_p, flags);
3059 : } else {
3060 0 : putCol (columns_p.flag_p, flags);
3061 : }
3062 0 : }
3063 :
3064 : void
3065 0 : VisibilityIteratorWriteImpl::setFlagCategory(const Array<Bool>& flagCategory)
3066 : {
3067 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3068 :
3069 0 : if (readImpl->useSlicer_p){
3070 0 : putCol(columns_p.flagCategory_p, readImpl->slicer_p, flagCategory);
3071 : }
3072 : else{
3073 0 : putCol(columns_p.flagCategory_p, flagCategory);
3074 : }
3075 0 : }
3076 :
3077 :
3078 : void
3079 0 : VisibilityIteratorWriteImpl::setFlagRow (const Vector<Bool> & rowflags)
3080 : {
3081 0 : putCol (columns_p.flagRow_p, rowflags);
3082 0 : }
3083 :
3084 : void
3085 0 : VisibilityIteratorWriteImpl::setVis (const Matrix<CStokesVector> & vis,
3086 : DataColumn whichOne)
3087 : {
3088 : // two problems: 1. channel selection -> we can only write to reference
3089 : // MS with 'processed' channels
3090 : // 2. polarization: there could be 1, 2 or 4 in the
3091 : // original data, predict () always gives us 4. We save what was there
3092 : // originally.
3093 :
3094 : // if (!preselected_p) {
3095 : // throw (AipsError ("VisIter::setVis (vis) - cannot change original data"));
3096 : // }
3097 :
3098 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3099 :
3100 0 : if (Int (vis.nrow ()) != readImpl->channelGroupSize_p) {
3101 0 : throw (AipsError ("VisIter::setVis (vis) - inconsistent number of channels"));
3102 : }
3103 : // we need to reform the vis matrix to a cube before we can use
3104 : // putColumn to a Matrix column
3105 0 : readImpl->cache_p.visCube_p.resize (readImpl->nPol_p, readImpl->channelGroupSize_p, readImpl->curNumRow_p);
3106 : Bool deleteIt;
3107 0 : Complex * p = readImpl->cache_p.visCube_p.getStorage (deleteIt);
3108 0 : for (uInt row = 0; row < readImpl->curNumRow_p; row++) {
3109 0 : for (Int chn = 0; chn < readImpl->channelGroupSize_p; chn++) {
3110 0 : const CStokesVector & v = vis (chn, row);
3111 0 : switch (readImpl->nPol_p) {
3112 0 : case 4:
3113 0 : *p++ = v (0);
3114 0 : *p++ = v (1);
3115 0 : *p++ = v (2);
3116 0 : *p++ = v (3);
3117 0 : break;
3118 0 : case 2:
3119 0 : *p++ = v (0);
3120 0 : *p++ = v (3);
3121 0 : break;
3122 0 : case 1:
3123 0 : *p++ = (v (0) + v (3)) / 2;
3124 0 : break;
3125 : }
3126 : }
3127 : }
3128 0 : if (readImpl->useSlicer_p) {
3129 0 : putDataColumn (whichOne, readImpl->slicer_p, readImpl->cache_p.visCube_p);
3130 : } else {
3131 0 : putDataColumn (whichOne, readImpl->cache_p.visCube_p);
3132 : }
3133 0 : }
3134 :
3135 : void
3136 60 : VisibilityIteratorWriteImpl::setVisAndFlag (const Cube<Complex> & vis,
3137 : const Cube<Bool> & flag,
3138 : DataColumn whichOne)
3139 : {
3140 60 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3141 :
3142 60 : if (readImpl->useSlicer_p) {
3143 0 : putDataColumn (whichOne, readImpl->slicer_p, vis);
3144 : } else {
3145 60 : putDataColumn (whichOne, vis);
3146 : }
3147 60 : if (readImpl->useSlicer_p) {
3148 0 : putCol (columns_p.flag_p, readImpl->slicer_p, flag);
3149 : } else {
3150 60 : putCol (columns_p.flag_p, flag);
3151 : }
3152 :
3153 60 : }
3154 :
3155 : void
3156 12 : VisibilityIteratorWriteImpl::setVis (const Cube<Complex> & vis, DataColumn whichOne)
3157 : {
3158 12 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3159 :
3160 12 : if (readImpl->useSlicer_p) {
3161 0 : putDataColumn (whichOne, readImpl->slicer_p, vis);
3162 : } else {
3163 12 : putDataColumn (whichOne, vis);
3164 : }
3165 12 : }
3166 :
3167 : void
3168 0 : VisibilityIteratorWriteImpl::setWeight (const Vector<Float> & weight)
3169 : {
3170 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3171 :
3172 : // No polarization dependence for now
3173 0 : Matrix<Float> polWeight;
3174 0 : readImpl->getCol (readImpl->columns_p.weight_p, polWeight);
3175 0 : for (Int i = 0; i < readImpl->nPol_p; i++) {
3176 0 : Vector<Float> r = polWeight.row (i);
3177 0 : r = weight;
3178 : }
3179 0 : putCol (columns_p.weight_p, polWeight);
3180 0 : }
3181 :
3182 : void
3183 0 : VisibilityIteratorWriteImpl::setWeightMat (const Matrix<Float> & weightMat)
3184 : {
3185 0 : putCol (columns_p.weight_p, weightMat);
3186 0 : }
3187 :
3188 : void
3189 0 : VisibilityIteratorWriteImpl::setWeightSpectrum (const Cube<Float> & weightSpectrum)
3190 : {
3191 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3192 :
3193 0 : if (! readImpl->columns_p.weightSpectrum_p.isNull ()) {
3194 0 : putCol (columns_p.weightSpectrum_p, weightSpectrum);
3195 : }
3196 0 : }
3197 :
3198 : void
3199 0 : VisibilityIteratorWriteImpl::setSigma (const Vector<Float> & sigma)
3200 : {
3201 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3202 :
3203 0 : Matrix<Float> sigmat;
3204 0 : readImpl->getCol (readImpl->columns_p.sigma_p, sigmat);
3205 0 : for (Int i = 0; i < readImpl->nPol_p; i++) {
3206 0 : Vector<Float> r = sigmat.row (i);
3207 0 : r = sigma;
3208 : }
3209 0 : putCol (columns_p.sigma_p, sigmat);
3210 0 : }
3211 :
3212 : void
3213 0 : VisibilityIteratorWriteImpl::setSigmaMat (const Matrix<Float> & sigMat)
3214 : {
3215 0 : putCol (columns_p.sigma_p, sigMat);
3216 0 : }
3217 :
3218 :
3219 : void
3220 0 : VisibilityIteratorWriteImpl::putDataColumn (DataColumn whichOne,
3221 : const Slicer & slicer,
3222 : const Cube<Complex> & data)
3223 : {
3224 : // Set the visibility (observed, model or corrected);
3225 : // deal with DATA and FLOAT_DATA seamlessly for observed data.
3226 :
3227 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3228 :
3229 0 : switch (whichOne) {
3230 :
3231 0 : case ROVisibilityIterator::Observed:
3232 0 : if (readImpl->floatDataFound_p) {
3233 0 : Cube<Float> dataFloat = real (data);
3234 0 : putCol (columns_p.floatVis_p, slicer, dataFloat);
3235 : } else {
3236 0 : putCol (columns_p.vis_p, slicer, data);
3237 : }
3238 0 : break;
3239 :
3240 0 : case ROVisibilityIterator::Corrected:
3241 0 : putCol (columns_p.corrVis_p, slicer, data);
3242 0 : break;
3243 :
3244 0 : case ROVisibilityIterator::Model:
3245 0 : putCol (columns_p.modelVis_p, slicer, data);
3246 0 : break;
3247 :
3248 0 : default:
3249 0 : Assert (false);
3250 : }
3251 0 : }
3252 :
3253 : void
3254 72 : VisibilityIteratorWriteImpl::putDataColumn (DataColumn whichOne,
3255 : const Cube<Complex> & data)
3256 : {
3257 : // Set the visibility (observed, model or corrected);
3258 : // deal with DATA and FLOAT_DATA seamlessly for observed data.
3259 :
3260 72 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3261 :
3262 72 : switch (whichOne) {
3263 :
3264 0 : case ROVisibilityIterator::Observed:
3265 0 : if (readImpl->floatDataFound_p) {
3266 0 : Cube<Float> dataFloat = real (data);
3267 0 : if (useCustomTileShape_p) setTileShape(readImpl->selRows_p,columns_p.floatVis_p,dataFloat.shape());
3268 0 : putCol (columns_p.floatVis_p, dataFloat);
3269 : } else {
3270 0 : if (useCustomTileShape_p) setTileShape(readImpl->selRows_p,columns_p.vis_p,data.shape());
3271 0 : putCol (columns_p.vis_p, data);
3272 : }
3273 0 : break;
3274 :
3275 66 : case ROVisibilityIterator::Corrected:
3276 66 : if (useCustomTileShape_p) setTileShape(readImpl->selRows_p,columns_p.corrVis_p,data.shape());
3277 66 : putCol (columns_p.corrVis_p, data);
3278 66 : break;
3279 :
3280 6 : case ROVisibilityIterator::Model:
3281 6 : if (useCustomTileShape_p) setTileShape(readImpl->selRows_p,columns_p.modelVis_p,data.shape());
3282 6 : putCol (columns_p.modelVis_p, data);
3283 6 : break;
3284 :
3285 0 : default:
3286 0 : Assert (false);
3287 : }
3288 72 : }
3289 :
3290 : void
3291 0 : VisibilityIteratorWriteImpl::putColScalar (ScalarColumn<Bool> &column, const Vector<Bool> &array)
3292 : {
3293 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3294 :
3295 0 : column.putColumnCells (readImpl->selRows_p, array);
3296 0 : return;
3297 : }
3298 :
3299 : void
3300 0 : VisibilityIteratorWriteImpl::putCol (ScalarColumn<Bool> &column, const Vector<Bool> &array)
3301 : {
3302 0 : putColScalar (column, array);
3303 0 : }
3304 :
3305 : void
3306 60 : VisibilityIteratorWriteImpl::putCol (ArrayColumn<Bool> &column, const Array<Bool> &array)
3307 : {
3308 60 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3309 :
3310 60 : column.putColumnCells (readImpl->selRows_p, array);
3311 60 : }
3312 :
3313 : void
3314 0 : VisibilityIteratorWriteImpl::putCol (ArrayColumn<Float> &column, const Array<Float> &array)
3315 : {
3316 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3317 :
3318 0 : column.putColumnCells (readImpl->selRows_p, array);
3319 0 : }
3320 :
3321 : void
3322 72 : VisibilityIteratorWriteImpl::putCol (ArrayColumn<Complex> &column, const Array<Complex> &array)
3323 : {
3324 72 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3325 :
3326 72 : column.putColumnCells (readImpl->selRows_p, array);
3327 72 : }
3328 :
3329 : void
3330 0 : VisibilityIteratorWriteImpl::putCol (ArrayColumn<Bool> &column,
3331 : const Slicer & slicer,
3332 : const Array<Bool> &array)
3333 : {
3334 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3335 :
3336 0 : column.putColumnCells (readImpl->selRows_p, slicer, array);
3337 0 : }
3338 :
3339 : void
3340 0 : VisibilityIteratorWriteImpl::putCol (ArrayColumn<Float> &column,
3341 : const Slicer & slicer,
3342 : const Array<Float> &array)
3343 : {
3344 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3345 :
3346 0 : column.putColumnCells (readImpl->selRows_p, slicer, array);
3347 0 : }
3348 :
3349 : void
3350 0 : VisibilityIteratorWriteImpl::putCol (ArrayColumn<Complex> &column,
3351 : const Slicer & slicer,
3352 : const Array<Complex> &array)
3353 : {
3354 0 : VisibilityIteratorReadImpl * readImpl = getReadImpl();
3355 :
3356 0 : column.putColumnCells (readImpl->selRows_p, slicer, array);
3357 0 : }
3358 :
3359 0 : template <class T> void VisibilityIteratorWriteImpl::setTileShape( RefRows &rowRef,
3360 : ArrayColumn<T> &outputDataCol,
3361 : const IPosition &arrayShape)
3362 : {
3363 0 : size_t nCorr = arrayShape(0);
3364 0 : size_t nChan = arrayShape(1);
3365 0 : ssize_t nRows = 1048576 / (sizeof(T)*nCorr*nChan);
3366 0 : IPosition outputPlaneShape(2,nCorr,nChan);
3367 0 : IPosition tileShape(3,nCorr,nChan,nRows);
3368 0 : outputDataCol.setShape(rowRef.firstRow(),outputPlaneShape,tileShape);
3369 :
3370 0 : return;
3371 : }
3372 :
3373 2 : Bool VisibilityIteratorWriteImpl::useCustomTileShape()
3374 : {
3375 2 : Bool ret = false;
3376 2 : String rank = EnvironmentVariable::get("OMPI_COMM_WORLD_RANK");
3377 2 : if (!rank.empty())
3378 : {
3379 0 : ret = true;
3380 : }
3381 :
3382 4 : return ret;
3383 : }
3384 :
3385 : void
3386 0 : VisibilityIteratorWriteImpl::putModel(const RecordInterface& rec, Bool iscomponentlist, Bool incremental)
3387 : {
3388 :
3389 : /*
3390 : Vector<Int> fields = getReadImpl()->msColumns().fieldId().getColumn();
3391 : const Int option = Sort::HeapSort | Sort::NoDuplicates;
3392 : const Sort::Order order = Sort::Ascending;
3393 :
3394 : Int nfields = GenSort<Int>::sort (fields, order, option);
3395 :
3396 : // Make sure we have the right size
3397 :
3398 : fields.resize(nfields, true);
3399 : */
3400 0 : Matrix<Int> combiIndex;
3401 0 : MSUtil::getIndexCombination(getReadImpl()->msColumns(), combiIndex);
3402 0 : Int msid = getReadImpl()->msId();
3403 :
3404 0 : Vector<Int> spws = getReadImpl()->msChannels_p.spw_p[msid];
3405 0 : Vector<Int> starts = getReadImpl()->msChannels_p.start_p[msid];
3406 0 : Vector<Int> nchan = getReadImpl()->msChannels_p.width_p[msid];
3407 0 : Vector<Int> incr = getReadImpl()->msChannels_p.inc_p[msid];
3408 0 : Matrix<Int> chansel(spws.nelements(),4);
3409 0 : chansel.column(0)=spws;
3410 0 : chansel.column(1)=starts;
3411 0 : chansel.column(2)=nchan;
3412 0 : chansel.column(3)=incr;
3413 :
3414 0 : CountedPtr<VisModelDataI> visModelData = VisModelDataI::create();
3415 :
3416 : //visModelData->putModelI (getReadImpl()->ms (), rec, fields, spws, starts, nchan, incr,
3417 : // iscomponentlist, incremental);
3418 :
3419 : //version 2.0 to keep track of scan and state_id selection
3420 0 : visModelData->putModelI (getReadImpl()->ms (), rec, combiIndex, chansel, iscomponentlist, incremental);
3421 :
3422 0 : }
3423 :
3424 :
3425 : void
3426 0 : VisibilityIteratorWriteImpl::writeBack (VisBuffer * vb)
3427 : {
3428 0 : if (backWriters_p.empty ()) {
3429 0 : initializeBackWriters ();
3430 : }
3431 :
3432 0 : VbDirtyComponents dirtyComponents = vb->dirtyComponentsGet ();
3433 :
3434 0 : for (VbDirtyComponents::const_iterator dirtyComponent = dirtyComponents.begin ();
3435 0 : dirtyComponent != dirtyComponents.end ();
3436 0 : dirtyComponent ++) {
3437 :
3438 0 : ThrowIf (backWriters_p.find (* dirtyComponent) == backWriters_p.end (),
3439 : String::format ("No writer defined for VisBuffer component %d", * dirtyComponent));
3440 0 : BackWriter * backWriter = backWriters_p [ * dirtyComponent];
3441 :
3442 : try {
3443 0 : (* backWriter) (this, vb);
3444 0 : } catch (AipsError & e) {
3445 0 : Rethrow (e, String::format ("Error while writing back VisBuffer component %d", * dirtyComponent));
3446 : }
3447 : }
3448 0 : }
3449 :
3450 : void
3451 0 : VisibilityIteratorWriteImpl::initializeBackWriters ()
3452 : {
3453 0 : backWriters_p [VisBufferComponents::Flag] =
3454 0 : makeBackWriter (& VisibilityIteratorWriteImpl::setFlag, & VisBuffer::flag);
3455 0 : backWriters_p [VisBufferComponents::FlagCube] =
3456 0 : makeBackWriter (& VisibilityIteratorWriteImpl::setFlag, & VisBuffer::flagCube);
3457 0 : backWriters_p [VisBufferComponents::FlagRow] =
3458 0 : makeBackWriter (& VisibilityIteratorWriteImpl::setFlagRow, & VisBuffer::flagRow);
3459 0 : backWriters_p [VisBufferComponents::FlagCategory] =
3460 0 : makeBackWriter (& VisibilityIteratorWriteImpl::setFlagCategory, & VisBuffer::flagCategory);
3461 0 : backWriters_p [VisBufferComponents::Sigma] =
3462 0 : makeBackWriter (& VisibilityIteratorWriteImpl::setSigma, & VisBuffer::sigma);
3463 0 : backWriters_p [VisBufferComponents::SigmaMat] =
3464 0 : makeBackWriter (& VisibilityIteratorWriteImpl::setSigmaMat, & VisBuffer::sigmaMat);
3465 0 : backWriters_p [VisBufferComponents::Weight] =
3466 0 : makeBackWriter (& VisibilityIteratorWriteImpl::setWeight, & VisBuffer::weight);
3467 0 : backWriters_p [VisBufferComponents::WeightMat] =
3468 0 : makeBackWriter (& VisibilityIteratorWriteImpl::setWeightMat, & VisBuffer::weightMat);
3469 :
3470 : // Now do the visibilities. These are slightly different since the setter requires an
3471 : // enum value.
3472 :
3473 0 : backWriters_p [VisBufferComponents::Observed] =
3474 0 : makeBackWriter2 (& VisibilityIteratorWriteImpl::setVis, & VisBuffer::visibility,
3475 : ROVisibilityIterator::Observed);
3476 0 : backWriters_p [VisBufferComponents::Corrected] =
3477 0 : makeBackWriter2 (& VisibilityIteratorWriteImpl::setVis, & VisBuffer::correctedVisibility,
3478 : ROVisibilityIterator::Corrected);
3479 0 : backWriters_p [VisBufferComponents::Model] =
3480 0 : makeBackWriter2 (& VisibilityIteratorWriteImpl::setVis, & VisBuffer::modelVisibility,
3481 : ROVisibilityIterator::Model);
3482 :
3483 0 : backWriters_p [VisBufferComponents::ObservedCube] =
3484 0 : makeBackWriter2 (& VisibilityIteratorWriteImpl::setVis, & VisBuffer::visCube,
3485 : ROVisibilityIterator::Observed);
3486 0 : backWriters_p [VisBufferComponents::CorrectedCube] =
3487 0 : makeBackWriter2 (& VisibilityIteratorWriteImpl::setVis, & VisBuffer::correctedVisCube,
3488 : ROVisibilityIterator::Corrected);
3489 0 : backWriters_p [VisBufferComponents::ModelCube] =
3490 0 : makeBackWriter2 (& VisibilityIteratorWriteImpl::setVis, & VisBuffer::modelVisCube,
3491 : ROVisibilityIterator::Model);
3492 :
3493 0 : }
3494 :
3495 : VisibilityIteratorWriteImpl::Columns &
3496 0 : VisibilityIteratorWriteImpl::Columns::operator= (const VisibilityIteratorWriteImpl::Columns & other)
3497 : {
3498 0 : flag_p.reference (other.flag_p);
3499 0 : flagCategory_p.reference (other.flagCategory_p);
3500 0 : flagRow_p.reference (other.flagRow_p);
3501 0 : vis_p.reference (other.vis_p);
3502 0 : floatVis_p.reference (other.floatVis_p);
3503 0 : modelVis_p.reference (other.modelVis_p);
3504 0 : corrVis_p.reference (other.corrVis_p);
3505 0 : weight_p.reference (other.weight_p);
3506 0 : weightSpectrum_p.reference (other.weightSpectrum_p);
3507 0 : sigma_p.reference (other.sigma_p);
3508 :
3509 0 : return * this;
3510 : }
3511 :
3512 :
3513 : } //# NAMESPACE CASA - END
3514 :
3515 :
3516 :
|
