Line data Source code
1 : //# Reweighter.cc
2 : //# Copyright (C) 1996-2007
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify
6 : //# it under the terms of the GNU General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or
8 : //# (at your option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful,
11 : //# but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 : //# GNU General Public License for more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License
16 : //# along with this library; if not, write to the Free Software
17 : //# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id: $
27 :
28 : // To make Timer reports in the inner loop of the simple copy,
29 : // uncomment the following line:
30 : //#define COPYTIMER
31 :
32 : #include <msvis/MSVis/Reweighter.h>
33 : #include <casacore/ms/MSSel/MSSelection.h>
34 : //#include <ms/MSSel/MSTimeGram.h>
35 : //#include <tables/TaQL/ExprNode.h>
36 : #include <casacore/tables/Tables/RefRows.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <casacore/casa/Arrays/Matrix.h>
39 : #include <casacore/casa/Arrays/Cube.h>
40 : #include <casacore/casa/Arrays/ArrayMath.h>
41 : #include <casacore/casa/Arrays/ArrayLogical.h>
42 : #include <casacore/casa/Arrays/ArrayUtil.h>
43 : #include <casacore/casa/Arrays/IPosition.h>
44 : #include <casacore/casa/Arrays/Slice.h>
45 : #include <casacore/casa/Logging/LogIO.h>
46 : #include <casacore/casa/OS/File.h>
47 : #include <casacore/casa/OS/HostInfo.h>
48 : #include <casacore/casa/OS/Memory.h> // Can be commented out along with
49 : // // Memory:: calls.
50 :
51 : //#ifdef COPYTIMER
52 : #include <casacore/casa/OS/Timer.h>
53 : //#endif
54 :
55 : #include <casacore/casa/Containers/Record.h>
56 : #include <casacore/casa/BasicSL/String.h>
57 : #include <casacore/casa/Utilities/Assert.h>
58 : #include <casacore/casa/Utilities/GenSort.h>
59 : #include <casacore/casa/System/AppInfo.h>
60 : #include <casacore/casa/System/ProgressMeter.h>
61 : #include <casacore/casa/Quanta/QuantumHolder.h>
62 : #include <msvis/MSVis/GroupProcessor.h>
63 : //#include <msvis/MSVis/VisSet.h>
64 : #include <msvis/MSVis/StatWT.h>
65 : #include <msvis/MSVis/VisBuffer.h>
66 : #include <msvis/MSVis/VisBufferComponents.h>
67 : #include <msvis/MSVis/VisChunkAverager.h>
68 : #include <msvis/MSVis/VisIterator.h>
69 : //#include <msvis/MSVis/VisibilityIterator.h>
70 : #include <casacore/tables/DataMan/IncrementalStMan.h>
71 : #include <casacore/tables/Tables/ScalarColumn.h>
72 : #include <casacore/tables/Tables/ScaColDesc.h>
73 : #include <casacore/tables/DataMan/StandardStMan.h>
74 : #include <casacore/tables/Tables/Table.h>
75 : #include <casacore/tables/Tables/PlainTable.h>
76 : #include <casacore/tables/Tables/TableDesc.h>
77 : #include <casacore/tables/Tables/TableInfo.h>
78 : #include <casacore/tables/Tables/TableLock.h>
79 : #include <casacore/tables/Tables/TableRecord.h>
80 : #include <casacore/tables/Tables/TableCopy.h>
81 : #include <casacore/tables/Tables/TableRow.h>
82 : #include <casacore/tables/DataMan/TiledColumnStMan.h>
83 : #include <casacore/tables/DataMan/TiledShapeStMan.h>
84 : #include <casacore/tables/DataMan/TiledDataStMan.h>
85 : #include <casacore/tables/DataMan/TiledStManAccessor.h>
86 : #include <casacore/ms/MeasurementSets/MSTileLayout.h>
87 : #include <sstream>
88 : #include <iomanip>
89 : #include <functional>
90 : #include <map>
91 : #include <set>
92 : #include <casacore/scimath/Mathematics/Smooth.h>
93 : #include <casacore/casa/Quanta/MVTime.h>
94 :
95 : using namespace casacore;
96 : namespace casa {
97 :
98 : typedef ROVisibilityIterator ROVisIter;
99 : typedef VisibilityIterator VisIter;
100 :
101 0 : Reweighter::Reweighter(const String& theMS, const Bool dorms, const uInt minsamp) :
102 0 : ms_p(MeasurementSet(theMS, Table::Update)),
103 0 : mssel_p(ms_p),
104 0 : dorms_p(dorms),
105 0 : minsamp_p(minsamp),
106 0 : msc_p(NULL),
107 0 : antennaSel_p(false),
108 0 : timeBin_p(-1.0),
109 0 : scanString_p(""),
110 0 : intentString_p(""),
111 0 : obsString_p(""),
112 0 : timeRange_p(""),
113 0 : arrayExpr_p(""),
114 0 : combine_p(""),
115 0 : fitspw_p("*"),
116 0 : outspw_p("*")
117 : {
118 0 : }
119 :
120 0 : Reweighter::~Reweighter()
121 : {
122 0 : delete msc_p;
123 0 : msc_p = NULL;
124 0 : }
125 :
126 0 : Bool Reweighter::selectSpw(std::set<Int>& spwset, Vector<Int>& chanStartv,
127 : Vector<Int>& chanEndv, Vector<Int>& chanStepv,
128 : const String& spwstr)
129 : {
130 0 : LogIO os(LogOrigin("Reweighter", "selectSpw()"));
131 :
132 0 : MSSelection mssel;
133 0 : String myspwstr(spwstr == "" ? "*" : spwstr);
134 :
135 0 : mssel.setSpwExpr(myspwstr);
136 :
137 : // Each row should have spw, start, stop, step
138 : // A single width is a default, but multiple widths should be used
139 : // literally.
140 0 : Matrix<Int> chansel = mssel.getChanList(&ms_p, 1);
141 :
142 0 : if(chansel.nrow() > 0) { // Use myspwstr if it selected anything...
143 0 : Vector<Int> spwv(chansel.column(0));
144 :
145 0 : uInt nspwsel = spwv.size();
146 0 : chanStartv.resize(nspwsel);
147 0 : chanEndv.resize(nspwsel);
148 0 : chanStepv.resize(nspwsel);
149 :
150 0 : chanStartv = chansel.column(1);
151 0 : chanEndv = chansel.column(2);
152 0 : chanStepv = chansel.column(3);
153 :
154 0 : uInt nspw = chanEndv.nelements();
155 :
156 0 : for(uInt k = 0; k < nspw; ++k){
157 0 : spwset.insert(spwv[k]);
158 :
159 0 : if(chanStepv[k] == 0) // CAS-2224, triggered by spw='0:2'
160 0 : chanStepv[k] = 1; // (as opposed to '0:2~2').
161 :
162 : //nchan[k] = 1 + (chanEndv[k] - chanStartv[k]) / chanStepv[k];
163 : //if(nchan[k] < 1)
164 : // nchan[k] = 1;
165 : }
166 0 : }
167 : else{ // select everything and rely on widths.
168 0 : MSSpWindowColumns mySpwTab(ms_p.spectralWindow());
169 0 : uInt nspw = mySpwTab.nrow();
170 0 : Vector<Int> nchan(nspw);
171 :
172 0 : nchan = mySpwTab.numChan().getColumn();
173 :
174 0 : chanStartv.resize(nspw);
175 0 : chanStepv.resize(nspw);
176 0 : for(uInt k = 0; k < nspw; ++k){
177 0 : spwset.insert(k);
178 0 : chanStartv[k] = 0;
179 0 : chanEndv[k] = nchan[k] - 1;
180 0 : chanStepv[k] = 1;
181 : }
182 0 : }
183 :
184 0 : mssel.getChanSlices(chanSlices_p, &ms_p, 1);
185 0 : return true;
186 0 : }
187 :
188 0 : Bool Reweighter::selectCorrelations(const String& corrstr)
189 : {
190 0 : LogIO os(LogOrigin("Reweighter", "selectCorrelations()"));
191 0 : MSSelection mssel;
192 0 : const Bool areSelecting = corrstr != "" && corrstr != "*";
193 :
194 0 : if(areSelecting)
195 0 : mssel.setPolnExpr(corrstr);
196 0 : mssel.getCorrSlices(corrSlices_p, &ms_p);
197 : //return getCorrTypes(polIDs_p, corrTypes_p, mssel, ms_p, areSelecting);
198 0 : return areSelecting;
199 0 : }
200 :
201 0 : Bool Reweighter::getCorrTypes(Vector<Int>& polIDs,
202 : Vector<Vector<Int> >& corrTypes,
203 : const MSColumns& msc)
204 : {
205 0 : Bool cando = true;
206 :
207 0 : polIDs = msc.dataDescription().polarizationId().getColumn();
208 0 : corrTypes = msc.polarization().corrType().getColumn();
209 0 : return cando;
210 : }
211 :
212 0 : Bool Reweighter::setmsselect(const String& fitspw, const String& outspw,
213 : const String& field, const String& baseline,
214 : const String& scan, const String& subarray,
215 : const String& correlation, const String& intent,
216 : const String& obs)
217 : {
218 0 : LogIO os(LogOrigin("Reweighter", "setmsselect()"));
219 : Bool ok;
220 :
221 : // All of the requested selection functions will be tried, even if an
222 : // earlier one has indicated its failure. This allows all of the selection
223 : // strings to be tested, yielding more complete feedback for the user
224 : // (fewer retries). This is a matter of taste, though. If the selections
225 : // turn out to be slow, this function should return on the first false.
226 :
227 : // Do the output spw selection first because the channel vectors are dummies;
228 : // for output we only care about the spws.
229 0 : if(!selectSpw(outspwset_p, fitStart_p, fitEnd_p, fitStep_p, outspw)){
230 : os << LogIO::SEVERE
231 : << "No spectral windows selected for reweighting."
232 0 : << LogIO::POST;
233 0 : ok = false;
234 : }
235 0 : chanSlices_p.resize(0);
236 : // Check for : in outspw, and warn about it, at the Python level.
237 :
238 0 : if(!selectSpw(fitspwset_p, fitStart_p, fitEnd_p, fitStep_p, fitspw)){
239 : os << LogIO::SEVERE
240 : << "No channels selected for calculating the scatter."
241 0 : << LogIO::POST;
242 0 : ok = false;
243 : }
244 0 : Record selrec = ms_p.msseltoindex(outspw, field);
245 0 : ok = selectSource(selrec.asArrayInt("field"));
246 :
247 0 : if(baseline != ""){
248 0 : Vector<Int> antid(0);
249 0 : Vector<String> antstr(1,baseline);
250 0 : selectAntenna(antid, antstr);
251 0 : }
252 0 : scanString_p = scan;
253 0 : intentString_p = intent;
254 0 : obsString_p = obs;
255 :
256 0 : if(subarray != "")
257 0 : selectArray(subarray);
258 :
259 0 : if(!selectCorrelations(correlation)){
260 : //os << LogIO::SEVERE << "No correlations selected." << LogIO::POST;
261 : //For now, false here means no selection== use all available correlations
262 : //so no nead to raise an error.
263 : //ok = false;
264 : }
265 :
266 0 : return ok;
267 0 : }
268 :
269 0 : Bool Reweighter::selectSource(const Vector<Int>& fieldid)
270 : {
271 0 : LogIO os(LogOrigin("Reweighter", "selectSource()"));
272 0 : Bool cando = true;
273 :
274 0 : if(fieldid.nelements() < 1){
275 0 : fieldId_p = Vector<Int>(1, -1);
276 : }
277 0 : else if(fieldid.nelements() > ms_p.field().nrow()){
278 : os << LogIO::SEVERE
279 : << "More fields were requested than are in the input MS.\n"
280 0 : << LogIO::POST;
281 0 : cando = false;
282 : }
283 0 : else if(max(fieldid) >= static_cast<Int>(ms_p.field().nrow())){
284 : // Arriving here is very unlikely since if fieldid came from MSSelection
285 : // bad fields were presumably already quietly dropped.
286 : os << LogIO::SEVERE
287 : << "At least 1 field was requested that is not in the input MS.\n"
288 0 : << LogIO::POST;
289 0 : cando = false;
290 : }
291 : else{
292 0 : fieldId_p = fieldid;
293 : }
294 :
295 0 : if(fieldId_p.nelements() == 1 && fieldId_p[0] < 0){
296 0 : fieldId_p.resize(ms_p.field().nrow());
297 0 : indgen(fieldId_p);
298 : }
299 0 : return cando;
300 0 : }
301 :
302 0 : MS::PredefinedColumns Reweighter::dataColStrToEnum(const String& col)
303 : {
304 0 : LogIO os(LogOrigin("Reweighter", "dataColStrToEnum()"));
305 0 : String tmpName(col);
306 :
307 0 : tmpName.upcase();
308 :
309 0 : MS::PredefinedColumns result = MS::DATA;
310 :
311 0 : if(tmpName == "OBSERVED" || tmpName == "DATA" || tmpName == MS::columnName(MS::DATA))
312 0 : result = MS::DATA;
313 0 : else if(tmpName == "FLOAT" || tmpName == "FLOAT_DATA" || tmpName == MS::columnName(MS::FLOAT_DATA))
314 0 : result = MS::FLOAT_DATA;
315 0 : else if(tmpName == "LAG" || tmpName == "LAG_DATA" || tmpName == MS::columnName(MS::LAG_DATA))
316 0 : result = MS::LAG_DATA;
317 0 : else if(tmpName == "MODEL" || tmpName == "MODEL_DATA" || tmpName == MS::columnName(MS::MODEL_DATA))
318 0 : result = MS::MODEL_DATA;
319 0 : else if(tmpName == "CORRECTED" || tmpName == "CORRECTED_DATA" ||
320 0 : tmpName == MS::columnName(MS::CORRECTED_DATA))
321 0 : result = MS::CORRECTED_DATA;
322 : else
323 : os << LogIO::WARN
324 : << "Unrecognized data column " << tmpName << "...using DATA."
325 0 : << LogIO::POST;
326 0 : return result;
327 0 : }
328 :
329 0 : void Reweighter::selectTime(Double timeBin, String timerng)
330 : {
331 0 : timeBin_p = timeBin;
332 0 : timeRange_p = timerng;
333 0 : }
334 :
335 0 : Bool Reweighter::reweight(String& colname, const String& combine)
336 : {
337 : //Bool retval = true;
338 0 : LogIO os(LogOrigin("Reweighter", "reweight()"));
339 :
340 : try{
341 0 : if(!makeSelection()){
342 : os << LogIO::SEVERE
343 : << "Failed on selection: the combination of spw, field, antenna, correlation, "
344 : << "and timerange may be invalid."
345 0 : << LogIO::POST;
346 0 : ms_p=MeasurementSet();
347 0 : return false;
348 : }
349 0 : msc_p = new MSColumns(mssel_p);
350 : // Note again the parseColumnNames() a few lines back that stops setupMS()
351 : // from being called if the MS doesn't have the requested columns.
352 :
353 0 : combine_p = combine;
354 :
355 : //Detaching the selected part
356 0 : ms_p = MeasurementSet();
357 :
358 0 : Block<Int> sort;
359 0 : if(!setSortOrder(sort, "obs,scan,state", false)){
360 : os << LogIO::WARN
361 : << "Something in the combine string is not supported for reweighting."
362 0 : << LogIO::POST;
363 : }
364 :
365 : // Aaargh...everywhere else VisIter is used a timeInterval of 0 is treated as
366 : // DBL_MAX, meaning that TIME can be in sort but effectively be ignored for
367 : // major chunking. Why couldn't they just have said DBL_MAX in the first
368 : // place?
369 0 : ROVisibilityIterator viIn(mssel_p, sort, false, DBL_MIN);
370 :
371 : // Make sure it is initialized before any copies are made.
372 0 : viIn.originChunks();
373 :
374 :
375 : // NB: The following correlation selection---which doesn't work
376 : // in multi-corrshape datasets---has been rendered meaningless
377 : // within StatWT. In general, it is probably not desirable to
378 : // make the weight disposition correlation-dependent, in any
379 : // case. As such, correlation selection has been disabled
380 : // in the user interface, too. (gmoellen, 2015Aug12)
381 :
382 : // Make a list of indices of selected correlations
383 0 : vector<uInt> selcorrs;
384 0 : uInt nSelCorrs=corrSlices_p[0].nelements();
385 0 : if (nSelCorrs>0) {
386 0 : for (uInt i = 0; i < nSelCorrs; i++) {
387 0 : Slice slice=corrSlices_p[0][i];
388 0 : selcorrs.push_back(slice.start());
389 : }
390 : }
391 : else { // no selection == all
392 0 : for (Int i = 0; i < viIn.nCorr(); i++) selcorrs.push_back(static_cast<uInt>(i));
393 : }
394 :
395 0 : StatWT statwt(viIn, dataColStrToEnum(colname), fitspw_p, outspw_p, dorms_p,
396 0 : minsamp_p,selcorrs);
397 0 : GroupProcessor gp(viIn, &statwt);
398 :
399 0 : gp.go();
400 :
401 : // There should be now be statistically determined sigmas and weights for
402 : // each selected row. If smoothing is wanted, i.e. convolving by a
403 : // gaussian in time (taking the ends into account), do it now with a pass
404 : // over just WEIGHT, SIGMA, and the flags.
405 : //
406 : // In theory two passes is a little less accurate than gathering all the
407 : // data for a time window and then calculating the variance, but I show
408 : // here that the loss of accuracy is very small at worst, and in practice
409 : // two passes is much better.
410 : //
411 : // The variance of \sigma^2 is \sigma^4 (\frac{2}{n - 1} + \frac{k}{n}),
412 : // where the kurtosis k is 0 for a normal distribution. It will be dropped
413 : // from here on since it does not affect the argument.
414 : //
415 : // So for a sample of mn visibilities (gather up all the data for a time
416 : // window before calculating the variances), the variance of \sigma^2 is
417 : // 2 \sigma^4 / (mn - 1).
418 : //
419 : // If the variance of the time window is instead calculated from the
420 : // variances of m groups of size n, the variance of the variance is
421 : // 2 \sigma^4 / [m (n - 1)]
422 : //
423 : // Thus there is a difference, but a small one if n >> 1 (i.e. a reasonable
424 : // number of channels). More importantly, the m groups of n method
425 : // * is robust against the n groups having different means.
426 : // The only likely way that all the groups would have exactly the same
427 : // mean is if the mean is 0, in which case rms should be used instead of
428 : // stddev and the relevant sample size is simply mn no matter how
429 : // they're grouped. (I think; untested)
430 : // * m groups of n is a lot more flexible for programming and needs less
431 : // memory. It does mean some extra I/O for another pass, but it's only
432 : // over WEIGHT, SIGMA, FLAG_ROW, and unfortunately FLAG.
433 :
434 0 : return true;
435 0 : }
436 0 : catch(AipsError x){
437 0 : ms_p = MeasurementSet();
438 0 : throw(x);
439 0 : }
440 0 : catch(...){
441 0 : ms_p = MeasurementSet();
442 0 : throw(AipsError("Unknown exception caught"));
443 0 : }
444 0 : }
445 :
446 0 : void Reweighter::makeUnionSpw()
447 : {
448 0 : std::set<Int> unionset = fitspwset_p;
449 0 : std::set<Int>::iterator it;
450 0 : uInt i = 0;
451 :
452 0 : for(it = outspwset_p.begin(); it != outspwset_p.end(); ++it)
453 0 : unionset.insert(*it);
454 0 : unionspw_p.resize(unionset.size());
455 :
456 0 : for(it = unionset.begin(); it != unionset.end(); ++it)
457 0 : unionspw_p[i++] = *it;
458 0 : }
459 :
460 0 : Bool Reweighter::makeSelection()
461 : {
462 0 : LogIO os(LogOrigin("Reweighter", "makeSelection()"));
463 :
464 : //VisSet/MSIter will check if the SORTED exists
465 : //and resort if necessary
466 : {
467 : // Matrix<Int> noselection;
468 : // VisSet vs(ms_p, noselection);
469 0 : Block<Int> sort;
470 0 : ROVisibilityIterator(ms_p, sort);
471 0 : }
472 :
473 : const MeasurementSet *elms;
474 0 : elms = &ms_p;
475 0 : MeasurementSet sorted;
476 0 : if(ms_p.keywordSet().isDefined("SORTED_TABLE")){
477 0 : sorted = ms_p.keywordSet().asTable("SORTED_TABLE");
478 : //If ms is not writable and sort is a subselection...use original ms
479 0 : if( ms_p.nrow() == sorted.nrow())
480 0 : elms = &sorted;
481 : }
482 :
483 0 : MSSelection thisSelection;
484 0 : if(fieldId_p.nelements() > 0)
485 0 : thisSelection.setFieldExpr(MSSelection::indexExprStr(fieldId_p));
486 0 : if(unionspw_p.nelements() > 0)
487 0 : thisSelection.setSpwExpr(MSSelection::indexExprStr(unionspw_p));
488 0 : if(antennaSel_p){
489 0 : if(antennaId_p.nelements() > 0){
490 0 : thisSelection.setAntennaExpr(MSSelection::indexExprStr( antennaId_p ));
491 : }
492 0 : if(antennaSelStr_p[0] != "")
493 0 : thisSelection.setAntennaExpr(MSSelection::nameExprStr( antennaSelStr_p));
494 : }
495 0 : if(timeRange_p != "")
496 0 : thisSelection.setTimeExpr(timeRange_p);
497 :
498 0 : thisSelection.setScanExpr(scanString_p);
499 0 : thisSelection.setStateExpr(intentString_p);
500 0 : thisSelection.setObservationExpr(obsString_p);
501 0 : if(arrayExpr_p != "")
502 0 : thisSelection.setArrayExpr(arrayExpr_p);
503 0 : if(corrString_p != "")
504 0 : thisSelection.setPolnExpr(corrString_p);
505 :
506 0 : TableExprNode exprNode=thisSelection.toTableExprNode(elms);
507 0 : selTimeRanges_p = thisSelection.getTimeList();
508 :
509 : // Now remake the selected ms
510 0 : if(!(exprNode.isNull())){
511 0 : mssel_p = MeasurementSet((*elms)(exprNode));
512 : }
513 : else{
514 : // Null take all the ms ...setdata() blank means that
515 0 : mssel_p = MeasurementSet((*elms));
516 : }
517 : //mssel_p.rename(ms_p.tableName()+"/SELECTED_TABLE", Table::Scratch);
518 0 : if(mssel_p.nrow() == 0)
519 0 : return false;
520 :
521 : // Setup antNewIndex_p now that mssel_p is ready.
522 0 : if(antennaSel_p){
523 : // Watch out! getAntenna*List() and getBaselineList() return negative
524 : // numbers for negated antennas!
525 : //Vector<Int> selAnt1s(thisSelection.getAntenna1List());
526 : //Vector<Int> selAnt2s(thisSelection.getAntenna2List());
527 0 : ScalarColumn<Int> ant1c(mssel_p, MS::columnName(MS::ANTENNA1));
528 0 : ScalarColumn<Int> ant2c(mssel_p, MS::columnName(MS::ANTENNA2));
529 0 : Vector<Int> selAnts(ant1c.getColumn());
530 0 : uInt nAnts = selAnts.nelements();
531 :
532 0 : selAnts.resize(2 * nAnts, true);
533 0 : selAnts(Slice(nAnts, nAnts)) = ant2c.getColumn();
534 0 : nAnts = GenSort<Int>::sort(selAnts, Sort::Ascending, Sort::NoDuplicates);
535 0 : selAnts.resize(nAnts, true);
536 0 : Int maxAnt = max(selAnts);
537 :
538 0 : if(maxAnt < 0){
539 : os << LogIO::SEVERE
540 : << "The maximum selected antenna number, " << maxAnt
541 : << ", seems to be < 0."
542 0 : << LogIO::POST;
543 0 : return false;
544 : }
545 :
546 0 : Bool trivial = true;
547 0 : for(uInt k = 0; k < nAnts; ++k)
548 0 : trivial &= (selAnts[k] == static_cast<Int>(k)); // trivial = selAnts == indgen(nAnts)
549 : // It is possible to exclude baselines
550 0 : antennaSel_p = !trivial; // without excluding any antennas.
551 0 : } // This still gets tripped up by VLA:OUT.
552 :
553 0 : if(mssel_p.nrow() < ms_p.nrow()){
554 : os << LogIO::NORMAL
555 : << mssel_p.nrow() << " out of " << ms_p.nrow() << " rows are going to be"
556 : << " considered due to the selection criteria."
557 0 : << LogIO::POST;
558 : }
559 0 : return true;
560 0 : }
561 :
562 0 : Bool Reweighter::shouldWatch(Bool& conflict, const String& col,
563 : const String& uncombinable,
564 : const Bool verbose) const
565 : {
566 0 : Bool wantWatch = !combine_p.contains(col);
567 :
568 0 : if(!wantWatch && uncombinable.contains(col)){
569 0 : conflict = true;
570 0 : wantWatch = false;
571 :
572 0 : if(verbose){
573 0 : LogIO os(LogOrigin("Reweighter", "shouldWatch()"));
574 :
575 : os << LogIO::WARN
576 : << "Combining by " << col
577 : << " was requested, but it is not allowed by this operation and will be ignored."
578 0 : << LogIO::POST;
579 0 : }
580 : }
581 0 : return wantWatch;
582 : }
583 :
584 0 : Bool Reweighter::setSortOrder(Block<Int>& sort, const String& uncombinable,
585 : const Bool verbose) const
586 : {
587 0 : Bool conflict = false;
588 0 : uInt n_cols_to_watch = 7; // 3 + #(watchables), whether or not they are watched.
589 :
590 : // Already separated by the chunking.
591 : //const Bool watch_array(!combine_p.contains("arr")); // Pirate talk for "array".
592 :
593 0 : Bool watch_obs = shouldWatch(conflict, "obs", uncombinable, verbose);
594 0 : Bool watch_scan = shouldWatch(conflict, "scan", uncombinable, verbose);
595 0 : Bool watch_spw = shouldWatch(conflict, "spw", uncombinable, verbose);
596 0 : Bool watch_state = shouldWatch(conflict, "state", uncombinable, verbose);
597 :
598 : // if(watch_obs)
599 : // ++n_cols_to_watch;
600 : // if(watch_scan)
601 : // ++n_cols_to_watch;
602 : // if(watch_spw)
603 : // ++n_cols_to_watch;
604 : // if(watch_state)
605 : // ++n_cols_to_watch;
606 :
607 0 : uInt colnum = 1;
608 :
609 0 : sort.resize(n_cols_to_watch);
610 0 : sort[0] = MS::ARRAY_ID;
611 0 : if(watch_scan){
612 0 : sort[colnum] = MS::SCAN_NUMBER;
613 0 : ++colnum;
614 : }
615 0 : if(watch_state){
616 0 : sort[colnum] = MS::STATE_ID;
617 0 : ++colnum;
618 : }
619 0 : sort[colnum] = MS::FIELD_ID;
620 0 : ++colnum;
621 0 : if(watch_spw){
622 0 : sort[colnum] = MS::DATA_DESC_ID;
623 0 : ++colnum;
624 : }
625 0 : sort[colnum] = MS::TIME;
626 0 : ++colnum;
627 0 : if(watch_obs){
628 0 : sort[colnum] = MS::OBSERVATION_ID;
629 0 : ++colnum;
630 : }
631 :
632 : // Now all the axes that should be combined should be added, so that they end
633 : // up in the same chunk.
634 0 : if(!watch_scan){
635 0 : sort[colnum] = MS::SCAN_NUMBER;
636 0 : ++colnum;
637 : }
638 0 : if(!watch_state){
639 0 : sort[colnum] = MS::STATE_ID;
640 0 : ++colnum;
641 : }
642 0 : if(!watch_spw){
643 0 : sort[colnum] = MS::DATA_DESC_ID;
644 0 : ++colnum;
645 : }
646 0 : if(!watch_obs){
647 0 : sort[colnum] = MS::OBSERVATION_ID;
648 : //++colnum;
649 : }
650 :
651 0 : return !conflict;
652 : }
653 :
654 0 : const ArrayColumn<Complex>& Reweighter::right_column(const MSColumns *msclala,
655 : const MS::PredefinedColumns col)
656 : {
657 0 : if(col == MS::DATA)
658 0 : return msclala->data();
659 0 : else if(col == MS::MODEL_DATA)
660 0 : return msclala->modelData();
661 : // else if(col == MS::FLOAT_DATA) // Not complex.
662 : // return msclala->floatData();
663 0 : else if(col == MS::LAG_DATA)
664 0 : return msclala->lagData();
665 : else // The honored-by-time-if-nothing-else
666 0 : return msclala->correctedData(); // default.
667 : }
668 :
669 : } //#End casa namespace
|
