Line data Source code
1 : //# MSIter2.cc: Implementation of MSIter2.h
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 : #include <memory>
28 :
29 : #include <msvis/MSVis/MSIter2.h>
30 : #include <casacore/casa/Arrays/ArrayMath.h>
31 : #include <casacore/casa/Arrays/ArrayLogical.h>
32 : #include <casacore/casa/Arrays/Slice.h>
33 : #include <casacore/casa/Exceptions/Error.h>
34 : #include <casacore/tables/Tables/TableIter.h>
35 : #include <casacore/casa/Utilities/Assert.h>
36 : #include <casacore/casa/Utilities/BinarySearch.h>
37 : #include <casacore/casa/Utilities/GenSort.h>
38 : #include <casacore/casa/Arrays/Slicer.h>
39 : #include <casacore/ms/MeasurementSets/MSColumns.h>
40 : #include <casacore/ms/MSSel/MSSpwIndex.h>
41 : #include <casacore/measures/Measures.h>
42 : #include <casacore/measures/Measures/MeasTable.h>
43 : #include <casacore/measures/Measures/MPosition.h>
44 : #include <casacore/measures/Measures/MEpoch.h>
45 : #include <casacore/casa/Quanta/MVTime.h>
46 : #include <casacore/measures/Measures/Stokes.h>
47 : #include <casacore/tables/Tables/TableRecord.h>
48 : #include <casacore/casa/Logging/LogIO.h>
49 : #include <casacore/casa/iostream.h>
50 :
51 : using namespace casacore;
52 :
53 : namespace casa { //# NAMESPACE CASA - BEGIN
54 : namespace vi { //# NAMESPACE VI - BEGIN
55 :
56 168 : MSSmartInterval::MSSmartInterval(Double interval,
57 168 : Vector<Double>& timebounds) :
58 : MSInterval(interval),
59 168 : interval2_(interval),
60 168 : offset2_(0.0),
61 168 : timeBounds_(timebounds),
62 168 : nBounds_(timebounds.nelements()),
63 168 : zeroInterval_(abs(interval)<(2.0*DBL_MIN)), // NB: DBL_MIN is smallest positive
64 168 : found_(false),
65 168 : bidx_(0)
66 : {
67 :
68 : // If timebounds unspecified, set one element to zero
69 : // In this case, we are gridding time globally
70 168 : if (timeBounds_.nelements()==0) {
71 0 : timeBounds_.resize(1);
72 0 : timeBounds_(0)=0.0;
73 : }
74 : // Initialize offset2_ to the first timeBounds_
75 168 : offset2_=timeBounds_(0);
76 168 : }
77 :
78 342 : void MSSmartInterval::setOffset(Double offset)
79 : {
80 : // Used only for absolute reset of the offset w.r.t.
81 : // the bounts list
82 342 : if (offset==0.0) {
83 342 : bidx_=0;
84 342 : offset2_=timeBounds_(0);
85 : }
86 342 : }
87 :
88 179425475 : int MSSmartInterval::comp(const void * obj1, const void * obj2) const
89 : {
90 179425475 : double v1 = *(const Double*)obj1;
91 179425475 : double v2 = *(const Double*)obj2;
92 :
93 : // cout.precision(6);
94 : // cout << " " << Float(v1-3600.0*floor(v1/3600.0)) << " ";
95 :
96 : // Shortcut if values are equal.
97 : // In synthesis MSs, this captures most row comparisons, luckily!
98 179425475 : if (v1 == v2) return 0;
99 :
100 : // Shortcut if interval is trivially small, or supplied times
101 : // differ by more than the interval
102 : // TBD: move zeroInterval_ context to own function
103 : // and point to it (handle v1==v2 also)?
104 134672795 : if (zeroInterval_ || abs(v1-v2)>interval2_)
105 8538500 : return v1 < v2 ? -1 : 1;
106 :
107 : // Reaching here, v1 and v2 differ by less than the inverval_,
108 : // so work harder to discern if they are in the same interval,
109 : // with attention to where the specified timeBounds_ are.
110 :
111 : // Find the reference time for the 2nd specified value
112 : // This costs, at worst, nBounds_*log(nBounds_) _each_ execution,
113 : // which could be problematic....
114 126134295 : if (nBounds_>1) {
115 91338237 : if (v2<timeBounds_[bidx_]) {
116 : // search from beginning,
117 20554 : bidx_=binarySearch(found_,timeBounds_,v2,nBounds_,0);
118 20554 : if (!found_) --bidx_; // handle exact match
119 : }
120 176169946 : else if (((bidx_+1)<nBounds_) && // a forward bound available
121 84852263 : (v2>timeBounds_[bidx_+1])) { // value is beyond it
122 : // search from current position
123 25826 : bidx_=binarySearch(found_,timeBounds_,v2,nBounds_-bidx_,bidx_);
124 25826 : if (!found_) --bidx_; // handle exact match
125 : }
126 : // else current bidx_ is correct
127 :
128 : // Offset boundary identified
129 91338237 : offset2_=timeBounds_[bidx_];
130 :
131 : // If v1 is prior to the detected boundary...
132 91338237 : if (v1 < offset2_) return -1;
133 :
134 : // If v1 later than the next higher boundary...
135 160261086 : if ( ((bidx_+1)<nBounds_) && // a forward bound available
136 84829290 : (v1 > timeBounds_[bidx_+1]) ) return 1;
137 : }
138 : else {
139 : // this should match initialization...
140 34796058 : bidx_=0;
141 34796058 : offset2_=timeBounds_[0];
142 : }
143 :
144 : // Reaching here, both values are in the same timeBounds_ interval;
145 : // As last resort, determine if v1 and v2 are in different gridded
146 : // intervals within the current timeBounds_ interval
147 :
148 : // The times are binned in bins with a width of interval_.
149 110227854 : t1_ = floor((v1 - offset2_) / interval2_);
150 110227854 : t2_ = floor((v2 - offset2_) / interval2_);
151 :
152 110227854 : return (t1_==t2_ ? 0 : (t1_<t2_ ? -1 : 1));
153 : }
154 :
155 :
156 0 : MSIter2::MSIter2():MSIter() {}
157 :
158 0 : MSIter2::MSIter2(const MeasurementSet& ms,
159 : const Block<Int>& sortColumns,
160 : Double timeInterval,
161 : Bool addDefaultSortColumns,
162 0 : Bool storeSorted)
163 0 : : MSIter(ms,sortColumns,timeInterval,addDefaultSortColumns,storeSorted)
164 0 : {}
165 :
166 :
167 168 : MSIter2::MSIter2(const Block<MeasurementSet>& mss,
168 : const Block<Int>& sortColumns,
169 : Double timeInterval,
170 : Bool addDefaultSortColumns,
171 168 : Bool storeSorted)
172 168 : : MSIter()
173 : {
174 168 : bms_p=mss;
175 168 : curMS_p=0;
176 168 : lastMS_p=-1;
177 168 : storeSorted_p=storeSorted;
178 168 : interval_p=timeInterval;
179 168 : prevFirstTimeStamp_p = -1.0;
180 168 : this->construct2(sortColumns,addDefaultSortColumns);
181 :
182 168 : }
183 :
184 342 : void MSIter2::origin()
185 : {
186 : // Reset time comparer, if present
187 : // (Ensures repeatability!)
188 342 : if (timeComp_p) timeComp_p->setOffset(0.0);
189 :
190 : // Call conventional
191 342 : MSIter::origin();
192 342 : }
193 :
194 :
195 :
196 168 : void MSIter2::construct2(const Block<Int>& sortColumns,
197 : Bool addDefaultSortColumns)
198 : {
199 168 : This = (MSIter2*)this;
200 168 : nMS_p=bms_p.nelements();
201 168 : if (nMS_p==0) throw(AipsError("MSIter::construct - No input MeasurementSets"));
202 336 : for (size_t i=0; i<nMS_p; i++) {
203 168 : if (bms_p[i].nrow()==0) {
204 0 : throw(AipsError("MSIter::construct - Input MeasurementSet.has zero selected rows"));
205 : }
206 : }
207 168 : tabIter_p.resize(nMS_p);
208 168 : tabIterAtStart_p.resize(nMS_p);
209 : // 'sort out' the sort orders
210 : // We normally require the table to be sorted on ARRAY_ID and FIELD_ID,
211 : // DATA_DESC_ID and TIME for the correct operation of the
212 : // VisibilityIterator (it needs to know when any of these changes to
213 : // be able to give the correct coordinates with the data)
214 : // If these columns are not explicitly sorted on, they will be added
215 : // BEFORE any others, unless addDefaultSortColumns=False
216 :
217 168 : Block<Int> cols;
218 : // try to reuse the existing sorted table if we didn't specify
219 : // any sortColumns
220 168 : if (sortColumns.nelements()==0 &&
221 168 : bms_p[0].keywordSet().isDefined("SORT_COLUMNS")) {
222 : // note that we use the order of the first MS for all MS's
223 0 : Vector<String> colNames = bms_p[0].keywordSet().asArrayString("SORT_COLUMNS");
224 0 : uInt n=colNames.nelements();
225 0 : cols.resize(n);
226 0 : for (uInt i=0; i<n; i++) cols[i]=MS::columnType(colNames(i));
227 0 : } else {
228 168 : cols=sortColumns;
229 : }
230 :
231 168 : timeInSort_p=False, arrayInSort_p=False, ddInSort_p=False, fieldInSort_p=False;
232 168 : bool scanSeen = false;
233 168 : Int nCol=0;
234 1089 : for (uInt i=0; i<cols.nelements(); i++) {
235 1842 : if (cols[i]>0 &&
236 921 : cols[i]<MS::NUMBER_PREDEFINED_COLUMNS) {
237 921 : if (cols[i]==MS::ARRAY_ID && !arrayInSort_p) { arrayInSort_p=True; nCol++; }
238 921 : if (cols[i]==MS::FIELD_ID && !fieldInSort_p) { fieldInSort_p=True; nCol++; }
239 921 : if (cols[i]==MS::DATA_DESC_ID && !ddInSort_p) { ddInSort_p=True; nCol++; }
240 921 : if (cols[i]==MS::TIME && !timeInSort_p) { timeInSort_p=True; nCol++; }
241 921 : if (cols[i]==MS::SCAN_NUMBER && !scanSeen) { scanSeen=True; }
242 : } else {
243 0 : throw(AipsError("MSIter() - invalid sort column"));
244 : }
245 : }
246 168 : Block<String> columns;
247 :
248 168 : Int iCol=0;
249 168 : if (addDefaultSortColumns) {
250 168 : columns.resize(cols.nelements()+4-nCol);
251 168 : if (!arrayInSort_p) {
252 : // add array if it's not there
253 0 : columns[iCol++]=MS::columnName(MS::ARRAY_ID);
254 : }
255 168 : if (!fieldInSort_p) {
256 : // add field if it's not there
257 0 : columns[iCol++]=MS::columnName(MS::FIELD_ID);
258 : }
259 168 : if (!ddInSort_p) {
260 : // add dd if it's not there
261 0 : columns[iCol++]=MS::columnName(MS::DATA_DESC_ID);
262 : }
263 168 : if (!timeInSort_p) {
264 : // add time if it's not there
265 0 : columns[iCol++]=MS::columnName(MS::TIME);
266 0 : timeInSort_p = True;
267 : }
268 : } else {
269 0 : columns.resize(cols.nelements());
270 : }
271 1089 : for (uInt i=0; i<cols.nelements(); i++) {
272 921 : columns[iCol++]=MS::columnName(MS::PredefinedColumns(cols[i]));
273 : }
274 :
275 168 : if (interval_p==0.0) {
276 0 : interval_p=DBL_MAX; // semi infinite
277 : } else {
278 : // assume that we want to sort on time if interval is set
279 168 : if (!timeInSort_p) {
280 0 : columns.resize(columns.nelements()+1);
281 0 : columns[iCol++]=MS::columnName(MS::TIME);
282 : }
283 : }
284 :
285 : // now find the time column and set the compare function
286 168 : Block<std::shared_ptr<BaseCompare> > objComp(columns.nelements());
287 168 : Block<Int> sortOrd(columns.nelements());
288 168 : timeComp_p = 0;
289 168 : Bool fieldBounded(false);
290 168 : Bool scanBounded=scanSeen;
291 1089 : for (uInt i=0; i<columns.nelements(); i++) {
292 :
293 : // Meaningful if this occurs before TIME
294 921 : if (columns[i]==MS::columnName(MS::FIELD_ID)) fieldBounded=true;
295 :
296 921 : if (columns[i]==MS::columnName(MS::TIME)) {
297 :
298 168 : Vector<Double> timebounds(0);
299 : //this->discernEnforcedTimeBounds(timebounds,scanBounded);
300 : //this->discernEnforcedTimeBounds(timebounds,scanBounded,fieldBounded);
301 168 : this->discernEnforcedTimeBounds(timebounds,scanBounded,fieldBounded,interval_p);
302 :
303 168 : timeComp_p = std::make_shared<MSSmartInterval>(interval_p,timebounds);
304 : //timeComp_p = new MSInterval(interval_p);
305 168 : objComp[i] = timeComp_p;
306 168 : }
307 921 : sortOrd[i]=Sort::Ascending;
308 : }
309 168 : Block<Int> orders(columns.nelements(),TableIterator::Ascending);
310 :
311 : // Store the sorted table for future access if possible,
312 : // reuse it if already there
313 336 : for (size_t i=0; i<nMS_p; i++) {
314 168 : Bool useIn=False, store=False, useSorted=False;
315 168 : Table sorted;
316 : // check if we already have a sorted table consistent with the requested
317 : // sort order
318 336 : if (!bms_p[i].keywordSet().isDefined("SORT_COLUMNS") ||
319 176 : !bms_p[i].keywordSet().isDefined("SORTED_TABLE") ||
320 176 : bms_p[i].keywordSet().asArrayString("SORT_COLUMNS").nelements()!=
321 352 : columns.nelements() ||
322 176 : !allEQ(bms_p[i].keywordSet().asArrayString("SORT_COLUMNS"),
323 176 : Vector<String>(columns.begin(),columns.end()))) {
324 : // if not, sort and store it (if possible)
325 168 : store=(bms_p[i].isWritable() && (bms_p[i].tableType() != Table::Memory));
326 : } else {
327 0 : sorted = bms_p[i].keywordSet().asTable("SORTED_TABLE");
328 : // if sorted table is smaller it can't be useful, remake it
329 0 : if (sorted.nrow() < bms_p[i].nrow()) store = bms_p[i].isWritable();
330 : else {
331 : // if input is a sorted subset of the stored sorted table
332 : // we can use the input in the iterator
333 0 : if (isSubSet(bms_p[i].rowNumbers(),sorted.rowNumbers())) {
334 0 : useIn=True;
335 : } else {
336 : // check if #rows in input table is the same as the base table
337 : // i.e., this is the entire table, if so, use sorted version instead
338 0 : String anttab = bms_p[i].antenna().tableName(); // see comments below
339 0 : Table base (anttab.erase(anttab.length()-8));
340 0 : if (base.nrow()==bms_p[i].nrow()) {
341 0 : useSorted = True;
342 : } else {
343 0 : store=bms_p[i].isWritable();
344 : }
345 0 : }
346 : }
347 : }
348 :
349 168 : if (!useIn && !useSorted) {
350 : // we have to resort the input
351 168 : if (aips_debug) cout << "MSIter::construct - resorting table"<<endl;
352 168 : sorted = bms_p[i].sort(columns, objComp, sortOrd, Sort::QuickSort);
353 : }
354 :
355 : // Only store if globally requested _and_ locally decided
356 168 : if (storeSorted_p && store) {
357 : // We need to get the name of the base table to add a persistent
358 : // subtable (the ms used here might be a reference table)
359 : // There is no table function to get this, so we use the name of
360 : // the antenna subtable to get at it.
361 0 : String anttab = bms_p[i].antenna().tableName();
362 0 : sorted.rename(anttab.erase(anttab.length()-7)+"SORTED_TABLE",Table::New);
363 0 : sorted.flush();
364 0 : bms_p[i].rwKeywordSet().defineTable("SORTED_TABLE",sorted);
365 0 : bms_p[i].rwKeywordSet().define("SORT_COLUMNS", Vector<String>(columns.begin( ),columns.end( )));
366 0 : }
367 :
368 : // create the iterator for each MS
369 : // at this stage either the input is sorted already or we are using
370 : // the sorted table, so the iterator can avoid sorting.
371 168 : if (useIn) {
372 0 : tabIter_p[i] = new TableIterator(bms_p[i],columns,objComp,orders,
373 0 : TableIterator::NoSort);
374 : } else {
375 168 : tabIter_p[i] = new TableIterator(sorted,columns,objComp,orders,
376 168 : TableIterator::NoSort);
377 : }
378 168 : tabIterAtStart_p[i]=True;
379 168 : }
380 168 : setMSInfo();
381 :
382 168 : }
383 :
384 0 : MSIter2::MSIter2(const MSIter2& other): MSIter(other)
385 : {
386 0 : operator=(other);
387 0 : }
388 :
389 336 : MSIter2::~MSIter2()
390 336 : {}
391 :
392 : MSIter2&
393 0 : MSIter2::operator=(const MSIter2& other)
394 : {
395 0 : if (this==&other) return *this;
396 0 : MSIter::operator=(other);
397 0 : return *this;
398 : }
399 :
400 198374 : void MSIter2::setState()
401 : {
402 198374 : setMSInfo();
403 198374 : if(newMS_p)
404 0 : checkFeed_p=True;
405 198374 : curTable_p=tabIter_p[curMS_p]->table();
406 198374 : colArray_p.attach(curTable_p,MS::columnName(MS::ARRAY_ID));
407 : // msc_p is already defined here (it is set in setMSInfo)
408 198374 : if(newMS_p)
409 0 : msc_p->antenna().mount().getColumn(antennaMounts_p,True);
410 :
411 198374 : if(!ddInSort_p)
412 : {
413 : // If Data Description is not in the sorting columns, then the DD, SPW, pol
414 : // can change between elements of the same iteration, so the safest
415 : // is to signal that it changes.
416 0 : newDataDescId_p = true;
417 0 : newSpectralWindowId_p = true;
418 0 : newPolarizationId_p = true;
419 0 : freqCacheOK_p= false;
420 :
421 : // This indicates that the current DD, SPW, Pol Ids are not computed.
422 : // Note that the last* variables are not set, since the new* variables
423 : // are unconditionally set to true.
424 : // These cur* vars wiil be computed lazily if it is needed, together
425 : // with some more vars set in cacheExtraDDInfo().
426 0 : curDataDescIdFirst_p = -1;
427 0 : curSpectralWindowIdFirst_p = -1;
428 0 : curPolarizationIdFirst_p = -1;
429 : }
430 : else
431 : {
432 : // First we cache the current DD, SPW, Pol since we know it changed
433 198374 : cacheCurrentDDInfo();
434 :
435 : // In this case we know that the last* variables were computed and
436 : // we can know whether there was a changed in these keywords by
437 : // comparing the two.
438 198374 : newDataDescId_p=(lastDataDescId_p!=curDataDescIdFirst_p);
439 198374 : newSpectralWindowId_p=(lastSpectralWindowId_p!=curSpectralWindowIdFirst_p);
440 198374 : newPolarizationId_p=(lastPolarizationId_p!=curPolarizationIdFirst_p);
441 :
442 198374 : lastDataDescId_p=curDataDescIdFirst_p;
443 198374 : lastSpectralWindowId_p = curSpectralWindowIdFirst_p;
444 198374 : lastPolarizationId_p = curPolarizationIdFirst_p;
445 :
446 : // Some extra information depends on the new* keywords, so compute
447 : // it now that new* have been set.
448 198374 : cacheExtraDDInfo();
449 : }
450 :
451 198374 : setArrayInfo();
452 198374 : feedInfoCached_p = false;
453 198374 : curFieldIdFirst_p=-1;
454 : //If field is not in the sorting columns, then the field id
455 : //can change between elements of the same iteration, so the safest
456 : //is to signal that it changes.
457 198374 : if(!fieldInSort_p)
458 0 : newFieldId_p = true;
459 : else
460 : {
461 198374 : setFieldInfo();
462 198374 : newFieldId_p=(lastFieldId_p!=curFieldIdFirst_p);
463 198374 : lastFieldId_p = curFieldIdFirst_p;
464 : }
465 :
466 :
467 : // If time binning, update the MSInterval's offset to account for glitches.
468 : // For example, if averaging to 5s and the input is
469 : // TIME STATE_ID INTERVAL
470 : // 0 0 1
471 : // 1 0 1
472 : // 2 1 1
473 : // 3 1 1
474 : // 4 1 1
475 : // 5 1 1
476 : // 6 1 1
477 : // 7 0 1
478 : // 8 0 1
479 : // 9 0 1
480 : // 10 0 1
481 : // 11 0 1
482 : // we want the output to be
483 : // TIME STATE_ID INTERVAL
484 : // 0.5 0 2
485 : // 4 1 5
486 : // 9 0 5
487 : // not what we'd get without the glitch fix:
488 : // TIME STATE_ID INTERVAL
489 : // 0.5 0 2
490 : // 3 1 3
491 : // 5.5 1 2
492 : // 8 0 3
493 : // 10.5 0 2
494 : //
495 : // Resetting the offset with each advance() might be too often, i.e. we might
496 : // need different spws to share the same offset. But in testing resetting
497 : // with each advance produces results more consistent with expectations than
498 : // either not resetting at all or resetting only
499 : // if(colTime_p(0) - 0.02 > timeComp_p->getOffset()).
500 : //
501 : // Don't need to do this with the MSSmartInterval?
502 : // if(timeComp_p && keyChange()=="SCAN_NUMBER"){
503 : // timeComp_p->setOffset(0.0);
504 : // }
505 198374 : }
506 :
507 :
508 0 : void MSIter2::discernEnforcedTimeBounds(Vector<Double>& timebounds,
509 : Bool scanBounded) {
510 :
511 0 : ScalarColumn<Double> timecol;
512 0 : std::map<Int,Double> scanmap;
513 0 : Int nscan(0);
514 0 : for (size_t iMS=0; iMS<nMS_p; iMS++) {
515 0 : TableIterator ti(bms_p[iMS],String("SCAN_NUMBER"));
516 :
517 0 : if (scanBounded) {
518 0 : while (!ti.pastEnd()) {
519 0 : timecol.attach(ti.table(),MS::columnName(MS::TIME));
520 0 : scanmap[nscan++]=timecol(0)-0.001;
521 0 : ti.next();
522 : }
523 : }
524 : else {
525 : // first time in the scan
526 0 : timecol.attach(ti.table(),MS::columnName(MS::TIME));
527 0 : scanmap[nscan++]=timecol(0)-0.001;
528 : }
529 0 : }
530 :
531 0 : timebounds.resize(nscan);
532 0 : for (Int iscan=0;iscan<nscan;++iscan)
533 0 : timebounds[iscan]=scanmap[iscan];
534 :
535 : //cout << "timebounds = " << timebounds-86400.0*floor(timebounds(0)/86400.0) << endl;
536 0 : }
537 :
538 :
539 0 : void MSIter2::discernEnforcedTimeBounds(Vector<Double>& timebounds,
540 : Bool scanBounded,
541 : Bool fieldBounded) {
542 :
543 0 : Int nsort(1+(fieldBounded?1:0));
544 0 : Block<String> cols(nsort);
545 0 : cols[0]="SCAN_NUMBER";
546 0 : if (fieldBounded) cols[1]="FIELD_ID";
547 0 : ScalarColumn<Double> timecol,expcol;
548 0 : ScalarColumn<Int> fieldcol,scancol;
549 0 : std::map<Int,Double> timemap;
550 0 : Int nscan(0);
551 0 : for (size_t iMS=0; iMS<nMS_p; iMS++) {
552 0 : TableIterator ti(bms_p[iMS],cols);
553 0 : Int thisScan(-1),thisField(-1),lastScan(-1),lastField(-1);
554 0 : while (!ti.pastEnd()) {
555 0 : timecol.attach(ti.table(),MS::columnName(MS::TIME));
556 0 : expcol.attach(ti.table(),MS::columnName(MS::EXPOSURE));
557 0 : scancol.attach(ti.table(),MS::columnName(MS::SCAN_NUMBER));
558 0 : fieldcol.attach(ti.table(),MS::columnName(MS::FIELD_ID));
559 0 : thisScan=scancol(0);
560 0 : thisField=fieldcol(0);
561 :
562 0 : if (nscan==0 || // first iteration
563 0 : (scanBounded && (thisScan!=lastScan)) || // per-scan and scan change
564 0 : (fieldBounded && ( (thisField!=lastField) || // per-field and field change
565 0 : ((thisScan-lastScan)>1) ) ) // or scan discont
566 : ) {
567 :
568 : /*
569 : cout << nscan << " "
570 : << thisScan << " "
571 : << thisField << " "
572 : << MVTime((timecol(0)-expcol(0)/2.0)/C::day).string(MVTime::YMD,7)
573 : << endl;
574 : */
575 :
576 0 : timemap[nscan++]=timecol(0)-0.001;
577 :
578 : }
579 0 : lastScan=thisScan;
580 0 : lastField=thisField;
581 0 : ti.next();
582 : }
583 0 : }
584 :
585 0 : timebounds.resize(nscan);
586 0 : for (Int iscan=0;iscan<nscan;++iscan)
587 0 : timebounds[iscan]=timemap[iscan];
588 :
589 0 : }
590 :
591 :
592 168 : void MSIter2::discernEnforcedTimeBounds(Vector<Double>& timebounds,
593 : Bool scanBounded,
594 : Bool fieldBounded,
595 : Double dt) {
596 :
597 168 : Int nsort(1+(fieldBounded?1:0));
598 168 : Block<String> cols(nsort);
599 168 : cols[0]="SCAN_NUMBER";
600 168 : if (fieldBounded) cols[1]="FIELD_ID";
601 168 : ScalarColumn<Double> timecol,expcol;
602 168 : ScalarColumn<Int> fieldcol,scancol;
603 168 : std::map<Int,Double> timemap;
604 168 : Int nscan(0);
605 168 : cout.precision(12);
606 336 : for (size_t iMS=0; iMS<nMS_p; iMS++) {
607 168 : TableIterator ti(bms_p[iMS],cols);
608 168 : Int thisScan(-1),thisField(-1),lastScan(-1),lastField(-1);
609 1342 : while (!ti.pastEnd()) {
610 1174 : timecol.attach(ti.table(),MS::columnName(MS::TIME));
611 1174 : expcol.attach(ti.table(),MS::columnName(MS::EXPOSURE));
612 1174 : scancol.attach(ti.table(),MS::columnName(MS::SCAN_NUMBER));
613 1174 : fieldcol.attach(ti.table(),MS::columnName(MS::FIELD_ID));
614 1174 : thisScan=scancol(0);
615 1174 : thisField=fieldcol(0);
616 :
617 1006 : if (nscan==0 || // first iteration
618 1006 : (scanBounded && (thisScan!=lastScan)) || // per-scan and scan change
619 2512 : (fieldBounded && (thisField!=lastField)) || // per-field and field change
620 348 : ( (thisScan-lastScan)>1 && // scan discontinuity _and_
621 1518 : (timecol(0)-timemap[nscan-1])>dt ) // previous interval expired
622 : ) {
623 :
624 : /*
625 : cout << nscan << " "
626 : << thisScan << " "
627 : << thisField << " "
628 : << MVTime((timecol(0)-expcol(0)/2.0)/C::day).string(MVTime::YMD,7)
629 : << endl;
630 : */
631 :
632 826 : timemap[nscan++]=(timecol(0)-0.01);
633 :
634 : }
635 1174 : lastScan=thisScan;
636 1174 : lastField=thisField;
637 1174 : ti.next();
638 : }
639 168 : }
640 :
641 168 : timebounds.resize(nscan);
642 994 : for (Int iscan=0;iscan<nscan;++iscan)
643 826 : timebounds[iscan]=timemap[iscan];
644 :
645 168 : }
646 :
647 :
648 : } //# NAMESPACE VI - END
649 : } //# NAMESPACE CASA - END
650 :
|
