Line data Source code
1 : //# StatWtTVI.cc: This file contains the implementation of the StatWtTVI class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All
5 : //# rights reserved.
6 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
7 : //#
8 : //# This library is free software; you can redistribute it and/or
9 : //# modify it under the terms of the GNU Lesser General Public
10 : //# License as published by the Free software Foundation; either
11 : //# version 2.1 of the License, or (at your option) any later version.
12 : //#
13 : //# This library is distributed in the hope that it will be useful,
14 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
15 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 : //# Lesser General Public License for more details.
17 : //#
18 : //# You should have received a copy of the GNU Lesser General Public
19 : //# License along with this library; if not, write to the Free Software
20 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
21 : //# MA 02111-1307 USA
22 :
23 : #include <mstransform/TVI/StatWtTVI.h>
24 :
25 : #include <casacore/casa/Arrays/ArrayLogical.h>
26 : #include <casacore/casa/Quanta/QuantumHolder.h>
27 : #include <casacore/ms/MSOper/MSMetaData.h>
28 : #include <casacore/tables/Tables/ArrColDesc.h>
29 :
30 : #include <mstransform/TVI/StatWtClassicalDataAggregator.h>
31 : #include <mstransform/TVI/StatWtFloatingWindowDataAggregator.h>
32 : #include <mstransform/TVI/StatWtVarianceAndWeightCalculator.h>
33 :
34 : #ifdef _OPENMP
35 : #include <omp.h>
36 : #endif
37 :
38 : #include <iomanip>
39 :
40 : using namespace casacore;
41 : using namespace casac;
42 :
43 : namespace casa {
44 : namespace vi {
45 :
46 : const String StatWtTVI::CHANBIN = "stchanbin";
47 :
48 0 : StatWtTVI::StatWtTVI(ViImplementation2 * inputVii, const Record &configuration)
49 0 : : TransformingVi2 (inputVii) {
50 : // Parse and check configuration parameters
51 : // Note: if a constructor finishes by throwing an exception, the memory
52 : // associated with the object itself is cleaned up there is no memory leak.
53 0 : ThrowIf(
54 : ! _parseConfiguration(configuration),
55 : "Error parsing StatWtTVI configuration"
56 : );
57 0 : LogIO log(LogOrigin("StatWtTVI", __func__));
58 0 : log << LogIO::NORMAL << "Using " << StatWtTypes::asString(_column)
59 0 : << " to compute weights" << LogIO::POST;
60 : // FIXME when the TVI framework has methods to
61 : // check for metadata, like the existence of
62 : // columns, remove references to the original MS
63 0 : const auto& origMS = ms();
64 : // FIXME uses original MS explicitly
65 0 : ThrowIf(
66 : (_column == StatWtTypes::CORRECTED || _column == StatWtTypes::RESIDUAL)
67 : && ! origMS.isColumn(MSMainEnums::CORRECTED_DATA),
68 : "StatWtTVI requires the MS to have a CORRECTED_DATA column. This MS "
69 : "does not"
70 : );
71 : // FIXME uses original MS explicitly
72 0 : ThrowIf(
73 : (_column == StatWtTypes::DATA || _column == StatWtTypes::RESIDUAL_DATA)
74 : && ! origMS.isColumn(MSMainEnums::DATA),
75 : "StatWtTVI requires the MS to have a DATA column. This MS does not"
76 : );
77 0 : _mustComputeSigma = (
78 0 : _column == StatWtTypes::DATA || _column == StatWtTypes::RESIDUAL_DATA
79 : );
80 : // FIXME uses original MS explicitly
81 0 : _updateWeight = ! _mustComputeSigma
82 0 : || (_mustComputeSigma && ! origMS.isColumn(MSMainEnums::CORRECTED_DATA));
83 0 : _noModel = (
84 0 : _column == StatWtTypes::RESIDUAL || _column == StatWtTypes::RESIDUAL_DATA
85 0 : ) && ! origMS.isColumn(MSMainEnums::MODEL_DATA)
86 0 : && ! origMS.source().isColumn(MSSourceEnums::SOURCE_MODEL);
87 : // Initialize attached VisBuffer
88 0 : setVisBuffer(createAttachedVisBuffer(VbRekeyable));
89 0 : }
90 :
91 0 : StatWtTVI::~StatWtTVI() {}
92 :
93 0 : Bool StatWtTVI::_parseConfiguration(const Record& config) {
94 0 : String field = CHANBIN;
95 0 : if (config.isDefined(field)) {
96 : // channel binning
97 0 : auto fieldNum = config.fieldNumber(field);
98 0 : switch (config.type(fieldNum)) {
99 0 : case DataType::TpArrayBool:
100 : // because this is the actual default variant type, no matter
101 : // what is specified in the xml
102 0 : ThrowIf(
103 : ! config.asArrayBool(field).empty(),
104 : "Unsupported data type for " + field
105 : );
106 0 : _setDefaultChanBinMap();
107 0 : break;
108 0 : case DataType::TpInt:
109 : Int binWidth;
110 0 : config.get(CHANBIN, binWidth);
111 0 : _setChanBinMap(binWidth);
112 0 : break;
113 0 : case DataType::TpString:
114 : {
115 0 : auto chanbin = config.asString(field);
116 0 : if (chanbin == "spw") {
117 : // bin using entire spws
118 0 : _setDefaultChanBinMap();
119 0 : break;
120 : }
121 : else {
122 0 : QuantumHolder qh(casaQuantity(chanbin));
123 0 : _setChanBinMap(qh.asQuantity());
124 0 : }
125 0 : break;
126 0 : }
127 0 : default:
128 0 : ThrowCc("Unsupported data type for " + field);
129 : }
130 : }
131 : else {
132 0 : _setDefaultChanBinMap();
133 : }
134 0 : field = "minsamp";
135 0 : if (config.isDefined(field)) {
136 0 : config.get(field, _minSamp);
137 0 : ThrowIf(_minSamp < 2, "Minimum size of sample must be >= 2.");
138 : }
139 0 : field = "combine";
140 0 : if (config.isDefined(field)) {
141 0 : ThrowIf(
142 : config.type(config.fieldNumber(field)) != TpString,
143 : "Unsupported data type for combine"
144 : );
145 0 : _combineCorr = config.asString(field).contains("corr");
146 : }
147 0 : field = "wtrange";
148 0 : if (config.isDefined(field)) {
149 0 : ThrowIf(
150 : config.type(config.fieldNumber(field)) != TpArrayDouble,
151 : "Unsupported type for field '" + field + "'"
152 : );
153 0 : auto myrange = config.asArrayDouble(field);
154 0 : if (! myrange.empty()) {
155 0 : ThrowIf(
156 : myrange.size() != 2,
157 : "Array specified in '" + field
158 : + "' must have exactly two values"
159 : );
160 0 : ThrowIf(
161 : casacore::anyLT(myrange, 0.0),
162 : "Both values specified in '" + field
163 : + "' array must be non-negative"
164 : );
165 : // order the values in case they aren't ordered
166 0 : std::set<Double> rangeset(myrange.begin(), myrange.end());
167 0 : ThrowIf(
168 : rangeset.size() == 1, "Values specified in '" + field
169 : + "' array must be unique"
170 : );
171 0 : auto iter = rangeset.begin();
172 0 : auto first = *iter;
173 0 : auto second = *(++iter);
174 0 : _wtrange.reset(new std::pair<Double, Double>(first, second));
175 0 : }
176 0 : }
177 0 : auto excludeChans = False;
178 0 : field = "excludechans";
179 0 : if (config.isDefined(field)) {
180 0 : ThrowIf(
181 : config.type(config.fieldNumber(field)) != TpBool,
182 : "Unsupported type for field '" + field + "'"
183 : );
184 0 : excludeChans = config.asBool(field);
185 : }
186 0 : field = "fitspw";
187 0 : if (config.isDefined(field)) {
188 0 : ThrowIf(
189 : config.type(config.fieldNumber(field)) != TpString,
190 : "Unsupported type for field '" + field + "'"
191 : );
192 0 : auto val = config.asString(field);
193 0 : if (! val.empty()) {
194 : // FIXME references underlying MS
195 0 : const auto& myms = ms();
196 0 : MSSelection sel(myms);
197 0 : sel.setSpwExpr(val);
198 0 : auto chans = sel.getChanList();
199 0 : auto nrows = chans.nrow();
200 0 : MSMetaData md(&myms, 50);
201 0 : auto nchans = md.nChans();
202 0 : IPosition start(3, 0);
203 0 : IPosition stop(3, 0);
204 0 : IPosition step(3, 1);
205 0 : for (size_t i=0; i<nrows; ++i) {
206 0 : auto row = chans.row(i);
207 0 : const auto& spw = row[0];
208 0 : if (_chanSelFlags.find(spw) == _chanSelFlags.end()) {
209 0 : _chanSelFlags[spw]
210 0 : = Cube<Bool>(1, nchans[spw], 1, ! excludeChans);
211 : }
212 0 : start[1] = row[1];
213 0 : ThrowIf(
214 : start[1] < 0, "Invalid channel selection in spw "
215 : + String::toString(spw))
216 : ;
217 0 : stop[1] = row[2];
218 0 : step[1] = row[3];
219 0 : Slicer slice(start, stop, step, Slicer::endIsLast);
220 0 : _chanSelFlags[spw](slice) = excludeChans;
221 0 : }
222 0 : }
223 0 : }
224 0 : field = "datacolumn";
225 0 : if (config.isDefined(field)) {
226 0 : ThrowIf(
227 : config.type(config.fieldNumber(field)) != TpString,
228 : "Unsupported type for field '" + field + "'"
229 : );
230 0 : auto val = config.asString(field);
231 0 : if (! val.empty()) {
232 0 : val.downcase();
233 0 : ThrowIf (
234 : ! (
235 : val.startsWith("c") || val.startsWith("d")
236 : || val.startsWith("residual") || val.startsWith("residual_")
237 : ),
238 : "Unsupported value for " + field + ": " + val
239 : );
240 0 : _column = val.startsWith("c") ? StatWtTypes::CORRECTED
241 0 : : val.startsWith("d") ? StatWtTypes::DATA
242 0 : : val.startsWith("residual_") ? StatWtTypes::RESIDUAL_DATA
243 : : StatWtTypes::RESIDUAL;
244 :
245 : }
246 0 : }
247 0 : field = "slidetimebin";
248 0 : ThrowIf(
249 : ! config.isDefined(field), "Config param " + field + " must be defined"
250 : );
251 0 : ThrowIf(
252 : config.type(config.fieldNumber(field)) != TpBool,
253 : "Unsupported type for field '" + field + "'"
254 : );
255 0 : _timeBlockProcessing = ! config.asBool(field);
256 0 : field = "timebin";
257 0 : ThrowIf(
258 : ! config.isDefined(field), "Config param " + field + " must be defined"
259 : );
260 0 : auto mytype = config.type(config.fieldNumber(field));
261 0 : ThrowIf(
262 : ! (
263 : mytype == TpString || mytype == TpDouble
264 : || mytype == TpInt
265 : ),
266 : "Unsupported type for field '" + field + "'"
267 : );
268 0 : switch(mytype) {
269 0 : case TpDouble: {
270 0 : _binWidthInSeconds.reset(new Double(config.asDouble(field)));
271 0 : break;
272 : }
273 0 : case TpInt:
274 0 : _nTimeStampsInBin.reset(new Int(config.asInt(field)));
275 0 : ThrowIf(
276 : *_nTimeStampsInBin <= 0,
277 : "Logic Error: nTimeStamps must be positive"
278 : );
279 0 : break;
280 0 : case TpString: {
281 0 : QuantumHolder qh(casaQuantity(config.asString(field)));
282 0 : _binWidthInSeconds.reset(
283 0 : new Double(getTimeBinWidthInSec(qh.asQuantity()))
284 : );
285 0 : break;
286 0 : }
287 0 : default:
288 0 : ThrowCc("Logic Error: Unhandled type for timebin");
289 :
290 : }
291 0 : _doClassicVIVB = _binWidthInSeconds && _timeBlockProcessing;
292 0 : _configureStatAlg(config);
293 0 : if (_doClassicVIVB) {
294 0 : _dataAggregator.reset(
295 : new StatWtClassicalDataAggregator(
296 0 : getVii(), _chanBins, _samples, _column, _noModel, _chanSelFlags,
297 0 : _wtStats, _wtrange, _combineCorr, _statAlg, _minSamp
298 0 : )
299 : );
300 : }
301 : else {
302 0 : _dataAggregator.reset(
303 : new StatWtFloatingWindowDataAggregator(
304 0 : getVii(), _chanBins, _samples, _column, _noModel, _chanSelFlags,
305 0 : _combineCorr, _wtStats, _wtrange, _binWidthInSeconds,
306 0 : _nTimeStampsInBin, _timeBlockProcessing, _statAlg, _minSamp
307 0 : )
308 : );
309 : }
310 0 : _dataAggregator->setMustComputeWtSp(_mustComputeWtSp);
311 0 : return True;
312 0 : }
313 :
314 0 : void StatWtTVI::_configureStatAlg(const Record& config) {
315 0 : String field = "statalg";
316 0 : if (config.isDefined(field)) {
317 0 : ThrowIf(
318 : config.type(config.fieldNumber(field)) != TpString,
319 : "Unsupported type for field '" + field + "'"
320 : );
321 0 : auto alg = config.asString(field);
322 0 : alg.downcase();
323 0 : if (alg.startsWith("cl")) {
324 0 : _statAlg.reset(
325 : new ClassicalStatistics<
326 : Double, Array<Float>::const_iterator,
327 : Array<Bool>::const_iterator, Array<Double>::const_iterator
328 0 : >()
329 : );
330 : }
331 : else {
332 : casacore::StatisticsAlgorithmFactory<
333 : Double, Array<Float>::const_iterator,
334 : Array<Bool>::const_iterator, Array<Double>::const_iterator
335 0 : > saf;
336 0 : if (alg.startsWith("ch")) {
337 0 : Int maxiter = -1;
338 0 : field = "maxiter";
339 0 : if (config.isDefined(field)) {
340 0 : ThrowIf(
341 : config.type(config.fieldNumber(field)) != TpInt,
342 : "Unsupported type for field '" + field + "'"
343 : );
344 0 : maxiter = config.asInt(field);
345 : }
346 0 : Double zscore = -1;
347 0 : field = "zscore";
348 0 : if (config.isDefined(field)) {
349 0 : ThrowIf(
350 : config.type(config.fieldNumber(field)) != TpDouble,
351 : "Unsupported type for field '" + field + "'"
352 : );
353 0 : zscore = config.asDouble(field);
354 : }
355 0 : saf.configureChauvenet(zscore, maxiter);
356 : }
357 0 : else if (alg.startsWith("f")) {
358 0 : auto center = FitToHalfStatisticsData::CMEAN;
359 0 : field = "center";
360 0 : if (config.isDefined(field)) {
361 0 : ThrowIf(
362 : config.type(config.fieldNumber(field)) != TpString,
363 : "Unsupported type for field '" + field + "'"
364 : );
365 0 : auto cs = config.asString(field);
366 0 : cs.downcase();
367 0 : if (cs == "mean") {
368 0 : center = FitToHalfStatisticsData::CMEAN;
369 : }
370 0 : else if (cs == "median") {
371 0 : center = FitToHalfStatisticsData::CMEDIAN;
372 : }
373 0 : else if (cs == "zero") {
374 0 : center = FitToHalfStatisticsData::CVALUE;
375 : }
376 : else {
377 0 : ThrowCc("Unsupported value for '" + field + "'");
378 : }
379 0 : }
380 0 : field = "lside";
381 0 : auto ud = FitToHalfStatisticsData::LE_CENTER;
382 0 : if (config.isDefined(field)) {
383 0 : ThrowIf(
384 : config.type(config.fieldNumber(field)) != TpBool,
385 : "Unsupported type for field '" + field + "'"
386 : );
387 0 : ud = config.asBool(field)
388 0 : ? FitToHalfStatisticsData::LE_CENTER
389 : : FitToHalfStatisticsData::GE_CENTER;
390 : }
391 0 : saf.configureFitToHalf(center, ud, 0);
392 : }
393 0 : else if (alg.startsWith("h")) {
394 0 : Double fence = -1;
395 0 : field = "fence";
396 0 : if (config.isDefined(field)) {
397 0 : ThrowIf(
398 : config.type(config.fieldNumber(field)) != TpDouble,
399 : "Unsupported type for field '" + field + "'"
400 : );
401 0 : fence = config.asDouble(field);
402 : }
403 0 : saf.configureHingesFences(fence);
404 : }
405 : else {
406 0 : ThrowCc("Unsupported value for 'statalg'");
407 : }
408 : // clone needed for CountedPtr -> shared_ptr hand off
409 0 : _statAlg.reset(saf.createStatsAlgorithm()->clone());
410 0 : }
411 0 : }
412 : else {
413 0 : _statAlg.reset(
414 : new ClassicalStatistics<
415 : Double, Array<Float>::const_iterator,
416 : Array<Bool>::const_iterator, Array<Double>::const_iterator
417 0 : >());
418 : }
419 0 : std::set<StatisticsData::STATS> stats {StatisticsData::VARIANCE};
420 0 : _statAlg->setStatsToCalculate(stats);
421 : // also configure the _wtStats object here
422 : // FIXME? Does not include exposure weighting
423 0 : _wtStats.reset(
424 : new ClassicalStatistics<
425 : Double, Array<Float>::const_iterator,
426 : Array<Bool>::const_iterator
427 0 : >()
428 : );
429 0 : stats.insert(StatisticsData::MEAN);
430 0 : _wtStats->setStatsToCalculate(stats);
431 0 : _wtStats->setCalculateAsAdded(True);
432 0 : }
433 :
434 0 : void StatWtTVI::_logUsedChannels() const {
435 : // FIXME uses underlying MS
436 0 : MSMetaData msmd(&ms(), 100.0);
437 0 : const auto nchan = msmd.nChans();
438 0 : LogIO log(LogOrigin("StatWtTVI", __func__));
439 0 : log << LogIO::NORMAL << "Weights are being computed using ";
440 0 : const auto cend = _chanSelFlags.cend();
441 0 : const auto nspw = _samples->size();
442 0 : uInt spwCount = 0;
443 0 : for (const auto& kv: *_samples) {
444 0 : const auto spw = kv.first;
445 0 : log << "SPW " << spw << ", channels ";
446 0 : const auto flagCube = _chanSelFlags.find(spw);
447 0 : if (flagCube == cend) {
448 0 : log << "0~" << (nchan[spw] - 1);
449 : }
450 : else {
451 0 : vector<pair<uInt, uInt>> startEnd;
452 0 : const auto flags = flagCube->second.tovector();
453 0 : bool started = false;
454 0 : std::unique_ptr<pair<uInt, uInt>> curPair;
455 0 : for (uInt j=0; j<nchan[spw]; ++j) {
456 0 : if (started) {
457 0 : if (flags[j]) {
458 : // found a bad channel, end current range
459 0 : startEnd.push_back(*curPair);
460 0 : started = false;
461 : }
462 : else {
463 : // found a "good" channel, update end of current range
464 0 : curPair->second = j;
465 : }
466 : }
467 0 : else if (! flags[j]) {
468 : // found a good channel, start new range
469 0 : started = true;
470 0 : curPair.reset(new pair<uInt, uInt>(j, j));
471 : }
472 : }
473 0 : if (curPair) {
474 0 : if (started) {
475 : // The last pair won't get added inside the previous loop,
476 : // so add it here
477 0 : startEnd.push_back(*curPair);
478 : }
479 0 : auto nPairs = startEnd.size();
480 0 : for (uInt i=0; i<nPairs; ++i) {
481 0 : log << startEnd[i].first << "~" << startEnd[i].second;
482 0 : if (i < nPairs - 1) {
483 0 : log << ", ";
484 : }
485 : }
486 : }
487 : else {
488 : // if the pointer never got set, all the channels are bad
489 0 : log << "no channels";
490 : }
491 0 : }
492 0 : if (spwCount < (nspw - 1)) {
493 0 : log << ";";
494 : }
495 0 : ++spwCount;
496 : }
497 0 : log << LogIO::POST;
498 0 : }
499 :
500 0 : void StatWtTVI::_setChanBinMap(const casacore::Quantity& binWidth) {
501 0 : if (! binWidth.isConform(Unit("Hz"))) {
502 0 : ostringstream oss;
503 : oss << "If specified as a quantity, channel bin width must have "
504 0 : << "frequency units. " << binWidth << " does not.";
505 0 : ThrowCc(oss.str());
506 0 : }
507 0 : ThrowIf(binWidth.getValue() <= 0, "channel bin width must be positive");
508 0 : MSMetaData msmd(&ms(), 100.0);
509 0 : auto chanFreqs = msmd.getChanFreqs();
510 0 : auto nspw = chanFreqs.size();
511 0 : auto binWidthHz = binWidth.getValue("Hz");
512 0 : for (uInt i=0; i<nspw; ++i) {
513 0 : auto cfs = chanFreqs[i].getValue("Hz");
514 0 : auto citer = cfs.begin();
515 0 : auto cend = cfs.end();
516 0 : StatWtTypes::ChanBin bin;
517 0 : bin.start = 0;
518 0 : bin.end = 0;
519 0 : uInt chanNum = 0;
520 0 : auto startFreq = *citer;
521 0 : auto nchan = cfs.size();
522 0 : for (; citer!=cend; ++citer, ++chanNum) {
523 : // both could be true, in which case both conditionals
524 : // must be executed
525 0 : if (abs(*citer - startFreq) > binWidthHz) {
526 : // add bin to list
527 0 : bin.end = chanNum - 1;
528 0 : _chanBins[i].push_back(bin);
529 0 : bin.start = chanNum;
530 0 : startFreq = *citer;
531 : }
532 0 : if (chanNum + 1 == nchan) {
533 : // add last bin
534 0 : bin.end = chanNum;
535 0 : _chanBins[i].push_back(bin);
536 : }
537 : }
538 0 : }
539 : // weight spectrum must be computed
540 0 : _mustComputeWtSp.reset(new Bool(True));
541 0 : }
542 :
543 0 : void StatWtTVI::_setChanBinMap(Int binWidth) {
544 0 : ThrowIf(binWidth < 1, "Channel bin width must be positive");
545 0 : MSMetaData msmd(&ms(), 100.0);
546 0 : auto nchans = msmd.nChans();
547 0 : auto nspw = nchans.size();
548 0 : StatWtTypes::ChanBin bin;
549 0 : for (uInt i=0; i<nspw; ++i) {
550 0 : auto lastChan = nchans[i]-1;
551 0 : for (uInt j=0; j<nchans[i]; j += binWidth) {
552 0 : bin.start = j;
553 0 : bin.end = min(j+binWidth-1, lastChan);
554 0 : _chanBins[i].push_back(bin);
555 : }
556 : }
557 : // weight spectrum must be computed
558 0 : _mustComputeWtSp.reset(new Bool(True));
559 0 : }
560 :
561 0 : void StatWtTVI::_setDefaultChanBinMap() {
562 : // FIXME uses underlying MS
563 0 : MSMetaData msmd(&ms(), 0.0);
564 0 : auto nchans = msmd.nChans();
565 0 : auto niter = nchans.begin();
566 0 : auto nend = nchans.end();
567 0 : Int i = 0;
568 0 : StatWtTypes::ChanBin bin;
569 0 : bin.start = 0;
570 0 : for (; niter!=nend; ++niter, ++i) {
571 0 : bin.end = *niter - 1;
572 0 : _chanBins[i].push_back(bin);
573 : }
574 0 : }
575 :
576 0 : Double StatWtTVI::getTimeBinWidthInSec(const casacore::Quantity& binWidth) {
577 0 : ThrowIf(
578 : ! binWidth.isConform(Unit("s")),
579 : "Time bin width unit must be a unit of time"
580 : );
581 0 : auto v = binWidth.getValue("s");
582 0 : checkTimeBinWidth(v);
583 0 : return v;
584 : }
585 :
586 0 : void StatWtTVI::checkTimeBinWidth(Double binWidth) {
587 0 : ThrowIf(binWidth <= 0, "time bin width must be positive");
588 0 : }
589 :
590 0 : void StatWtTVI::sigmaSpectrum(Cube<Float>& sigmaSp) const {
591 0 : if (_mustComputeSigma) {
592 : {
593 0 : Cube<Float> wtsp;
594 : // this computes _newWtsp, ignore wtsp
595 0 : weightSpectrum(wtsp);
596 0 : }
597 0 : sigmaSp = Float(1.0)/sqrt(_newWtSp);
598 0 : if (anyEQ(_newWtSp, Float(0))) {
599 0 : auto iter = sigmaSp.begin();
600 0 : auto end = sigmaSp.end();
601 0 : auto witer = _newWtSp.cbegin();
602 0 : for ( ; iter != end; ++iter, ++witer) {
603 0 : if (*witer == 0) {
604 0 : *iter = -1;
605 : }
606 : }
607 0 : }
608 : }
609 : else {
610 0 : TransformingVi2::sigmaSpectrum(sigmaSp);
611 : }
612 0 : }
613 :
614 0 : void StatWtTVI::weightSpectrum(Cube<Float>& newWtsp) const {
615 0 : ThrowIf(! _weightsComputed, "Weights have not been computed yet");
616 0 : if (! _dataAggregator->mustComputeWtSp()) {
617 0 : newWtsp.resize(IPosition(3, 0));
618 0 : return;
619 : }
620 0 : if (! _newWtSp.empty()) {
621 : // already calculated
622 0 : if (_updateWeight) {
623 0 : newWtsp = _newWtSp.copy();
624 : }
625 : else {
626 0 : TransformingVi2::weightSpectrum(newWtsp);
627 : }
628 0 : return;
629 : }
630 0 : _computeWeightSpectrumAndFlags();
631 0 : if (_updateWeight) {
632 0 : newWtsp = _newWtSp.copy();
633 : }
634 : else {
635 0 : TransformingVi2::weightSpectrum(newWtsp);
636 : }
637 : }
638 :
639 0 : void StatWtTVI::_computeWeightSpectrumAndFlags() const {
640 : size_t nOrigFlagged;
641 0 : auto mypair = _getLowerLayerWtSpFlags(nOrigFlagged);
642 0 : auto& wtsp = mypair.first;
643 0 : auto& flagCube = mypair.second;
644 0 : if (_dataAggregator->mustComputeWtSp() && wtsp.empty()) {
645 : // This can happen in preview mode if
646 : // WEIGHT_SPECTRUM doesn't exist or is empty
647 0 : wtsp.resize(flagCube.shape());
648 : }
649 0 : auto checkFlags = False;
650 0 : Vector<Int> ant1, ant2, spws;
651 0 : antenna1(ant1);
652 0 : antenna2(ant2);
653 0 : spectralWindows(spws);
654 0 : Vector<rownr_t> rowIDs;
655 0 : getRowIds(rowIDs);
656 0 : Vector<Double> exposures;
657 0 : exposure(exposures);
658 0 : _dataAggregator->weightSpectrumFlags(
659 : wtsp, flagCube, checkFlags, ant1, ant2, spws, exposures, rowIDs
660 : );
661 0 : if (checkFlags) {
662 0 : _nNewFlaggedPts += ntrue(flagCube) - nOrigFlagged;
663 : }
664 0 : _newWtSp = wtsp;
665 0 : _newFlag = flagCube;
666 0 : }
667 :
668 0 : std::pair<Cube<Float>, Cube<Bool>> StatWtTVI::_getLowerLayerWtSpFlags(
669 : size_t& nOrigFlagged
670 : ) const {
671 0 : auto mypair = std::make_pair(Cube<Float>(), Cube<Bool>());
672 0 : if (_dataAggregator->mustComputeWtSp()) {
673 0 : getVii()->weightSpectrum(mypair.first);
674 : }
675 0 : getVii()->flag(mypair.second);
676 0 : _nTotalPts += mypair.second.size();
677 0 : nOrigFlagged = ntrue(mypair.second);
678 0 : _nOrigFlaggedPts += nOrigFlagged;
679 0 : return mypair;
680 0 : }
681 :
682 0 : void StatWtTVI::sigma(Matrix<Float>& sigmaMat) const {
683 0 : if (_mustComputeSigma) {
684 0 : if (_newWt.empty()) {
685 0 : Matrix<Float> wtmat;
686 0 : weight(wtmat);
687 0 : }
688 0 : sigmaMat = Float(1.0)/sqrt(_newWt);
689 0 : if (anyEQ(_newWt, Float(0))) {
690 0 : Matrix<Float>::iterator iter = sigmaMat.begin();
691 0 : Matrix<Float>::iterator end = sigmaMat.end();
692 0 : Matrix<Float>::iterator witer = _newWt.begin();
693 0 : for ( ; iter != end; ++iter, ++witer) {
694 0 : if (*witer == 0) {
695 0 : *iter = -1;
696 : }
697 : }
698 0 : }
699 : }
700 : else {
701 0 : TransformingVi2::sigma(sigmaMat);
702 : }
703 0 : }
704 :
705 0 : void StatWtTVI::weight(Matrix<Float> & wtmat) const {
706 0 : ThrowIf(! _weightsComputed, "Weights have not been computed yet");
707 0 : if (! _newWt.empty()) {
708 0 : if (_updateWeight) {
709 0 : wtmat = _newWt.copy();
710 : }
711 : else {
712 0 : TransformingVi2::weight(wtmat);
713 : }
714 0 : return;
715 : }
716 0 : auto nrows = nRows();
717 0 : getVii()->weight(wtmat);
718 0 : if (_dataAggregator->mustComputeWtSp()) {
719 : // always use classical algorithm to get median for weights
720 : ClassicalStatistics<
721 : Double, Array<Float>::const_iterator, Array<Bool>::const_iterator
722 0 : > cs;
723 0 : Cube<Float> wtsp;
724 0 : Cube<Bool> flagCube;
725 : // this computes _newWtsP which is what we will use, so
726 : // just ignore wtsp
727 0 : weightSpectrum(wtsp);
728 0 : flag(flagCube);
729 0 : IPosition blc(3, 0);
730 0 : IPosition trc = _newWtSp.shape() - 1;
731 0 : const auto ncorr = _newWtSp.shape()[0];
732 0 : for (rownr_t i=0; i<nrows; ++i) {
733 0 : blc[2] = i;
734 0 : trc[2] = i;
735 0 : if (_combineCorr) {
736 0 : auto flags = flagCube(blc, trc);
737 0 : if (allTrue(flags)) {
738 0 : wtmat.column(i) = 0;
739 : }
740 : else {
741 0 : auto weights = _newWtSp(blc, trc);
742 0 : auto mask = ! flags;
743 0 : cs.setData(weights.begin(), mask.begin(), weights.size());
744 0 : wtmat.column(i) = cs.getMedian();
745 0 : }
746 0 : }
747 : else {
748 0 : for (uInt corr=0; corr<ncorr; ++corr) {
749 0 : blc[0] = corr;
750 0 : trc[0] = corr;
751 0 : auto weights = _newWtSp(blc, trc);
752 0 : auto flags = flagCube(blc, trc);
753 0 : if (allTrue(flags)) {
754 0 : wtmat(corr, i) = 0;
755 : }
756 : else {
757 0 : auto mask = ! flags;
758 0 : cs.setData(
759 0 : weights.begin(), mask.begin(), weights.size()
760 : );
761 0 : wtmat(corr, i) = cs.getMedian();
762 0 : }
763 0 : }
764 : }
765 : }
766 0 : }
767 : else {
768 : // the only way this can happen is if there is a single channel bin
769 : // for each baseline/spw pair
770 0 : _dataAggregator->weightSingleChanBin(wtmat, nrows);
771 : }
772 0 : _newWt = wtmat.copy();
773 0 : if (! _updateWeight) {
774 0 : wtmat = Matrix<Float>(wtmat.shape());
775 0 : TransformingVi2::weight(wtmat);
776 : }
777 : }
778 :
779 0 : void StatWtTVI::flag(Cube<Bool>& flagCube) const {
780 0 : ThrowIf(! _weightsComputed, "Weights have not been computed yet");
781 0 : if (! _newFlag.empty()) {
782 0 : flagCube = _newFlag.copy();
783 0 : return;
784 : }
785 0 : _computeWeightSpectrumAndFlags();
786 0 : flagCube = _newFlag.copy();
787 : }
788 :
789 0 : void StatWtTVI::flagRow(Vector<Bool>& flagRow) const {
790 0 : ThrowIf(! _weightsComputed, "Weights have not been computed yet");
791 0 : if (! _newFlagRow.empty()) {
792 0 : flagRow = _newFlagRow.copy();
793 0 : return;
794 : }
795 0 : Cube<Bool> flags;
796 0 : flag(flags);
797 0 : getVii()->flagRow(flagRow);
798 0 : auto nrows = nRows();
799 0 : for (rownr_t i=0; i<nrows; ++i) {
800 0 : flagRow[i] = allTrue(flags.xyPlane(i));
801 : }
802 0 : _newFlagRow = flagRow.copy();
803 0 : }
804 :
805 0 : void StatWtTVI::originChunks(Bool forceRewind) {
806 : // Drive next lower layer
807 0 : getVii()->originChunks(forceRewind);
808 0 : _weightsComputed = False;
809 0 : _dataAggregator->aggregate();
810 0 : _weightsComputed = True;
811 0 : _clearCache();
812 : // re-origin this chunk in next layer
813 : // (ensures wider scopes see start of the this chunk)
814 0 : getVii()->origin();
815 0 : }
816 :
817 0 : void StatWtTVI::nextChunk() {
818 : // Drive next lower layer
819 0 : getVii()->nextChunk();
820 0 : _weightsComputed = False;
821 0 : _dataAggregator->aggregate();
822 0 : _weightsComputed = True;
823 0 : _clearCache();
824 : // re-origin this chunk next layer
825 : // (ensures wider scopes see start of the this chunk)
826 0 : getVii()->origin();
827 0 : }
828 :
829 0 : void StatWtTVI::_clearCache() {
830 0 : _newWtSp.resize(0, 0, 0);
831 0 : _newWt.resize(0, 0);
832 0 : _newFlag.resize(0, 0, 0);
833 0 : _newFlagRow.resize(0);
834 0 : }
835 :
836 0 : Cube<Bool> StatWtTVI::_getResultantFlags(
837 : Cube<Bool>& chanSelFlagTemplate, Cube<Bool>& chanSelFlags,
838 : Bool& initTemplate, Int spw, const Cube<Bool>& flagCube
839 : ) const {
840 0 : if (_chanSelFlags.find(spw) == _chanSelFlags.cend()) {
841 : // no selection of channels to ignore
842 0 : return flagCube;
843 : }
844 0 : if (initTemplate) {
845 : // this can be done just once per chunk because all the rows
846 : // in the chunk are guaranteed to have the same spw
847 : // because each subchunk is guaranteed to have a single
848 : // data description ID.
849 0 : chanSelFlagTemplate = _chanSelFlags.find(spw)->second;
850 0 : initTemplate = False;
851 : }
852 0 : auto dataShape = flagCube.shape();
853 0 : chanSelFlags.resize(dataShape, False);
854 0 : auto ncorr = dataShape[0];
855 0 : auto nrows = dataShape[2];
856 0 : IPosition start(3, 0);
857 0 : IPosition end = dataShape - 1;
858 0 : Slicer sl(start, end, Slicer::endIsLast);
859 0 : for (uInt corr=0; corr<ncorr; ++corr) {
860 0 : start[0] = corr;
861 0 : end[0] = corr;
862 0 : for (Int row=0; row<nrows; ++row) {
863 0 : start[2] = row;
864 0 : end[2] = row;
865 0 : sl.setStart(start);
866 0 : sl.setEnd(end);
867 0 : chanSelFlags(sl) = chanSelFlagTemplate;
868 : }
869 : }
870 0 : return flagCube || chanSelFlags;
871 0 : }
872 :
873 0 : void StatWtTVI::initWeightSpectrum (const Cube<Float>& wtspec) {
874 : // Pass to next layer down
875 0 : getVii()->initWeightSpectrum(wtspec);
876 0 : }
877 :
878 0 : void StatWtTVI::initSigmaSpectrum (const Cube<Float>& sigspec) {
879 : // Pass to next layer down
880 0 : getVii()->initSigmaSpectrum(sigspec);
881 0 : }
882 :
883 :
884 0 : void StatWtTVI::writeBackChanges(VisBuffer2 *vb) {
885 : // Pass to next layer down
886 0 : getVii()->writeBackChanges(vb);
887 0 : }
888 :
889 0 : void StatWtTVI::summarizeFlagging() const {
890 0 : auto orig = (Double)_nOrigFlaggedPts/(Double)_nTotalPts*100;
891 0 : auto stwt = (Double)_nNewFlaggedPts/(Double)_nTotalPts*100;
892 0 : auto total = orig + stwt;
893 0 : LogIO log(LogOrigin("StatWtTVI", __func__));
894 : log << LogIO::NORMAL << "Originally, " << orig
895 : << "% of the data were flagged. StatWtTVI flagged an "
896 0 : << "additional " << stwt << "%." << LogIO::POST;
897 : log << LogIO::NORMAL << "TOTAL FLAGGED DATA AFTER RUNNING STATWT: "
898 0 : << total << "%" << LogIO::POST;
899 0 : log << LogIO::NORMAL << std::endl << LogIO::POST;
900 0 : if (_nOrigFlaggedPts == _nTotalPts) {
901 : log << LogIO::WARN << "IT APPEARS THAT ALL THE DATA IN THE INPUT "
902 : << "MS/SELECTION WERE FLAGGED PRIOR TO RUNNING STATWT"
903 0 : << LogIO::POST;
904 0 : log << LogIO::NORMAL << std::endl << LogIO::POST;
905 : }
906 0 : else if (_nOrigFlaggedPts + _nNewFlaggedPts == _nTotalPts) {
907 : log << LogIO::WARN << "IT APPEARS THAT STATWT FLAGGED ALL THE DATA "
908 : "IN THE REQUESTED SELECTION THAT WASN'T ORIGINALLY FLAGGED"
909 0 : << LogIO::POST;
910 0 : log << LogIO::NORMAL << std::endl << LogIO::POST;
911 : }
912 0 : String col0 = "SPECTRAL_WINDOW";
913 0 : String col1 = "SAMPLES_WITH_NON-ZERO_VARIANCE";
914 : String col2 = "SAMPLES_WHERE_REAL_PART_VARIANCE_DIFFERS_BY_>50%_FROM_"
915 0 : "IMAGINARY_PART";
916 0 : log << LogIO::NORMAL << col0 << " " << col1 << " " << col2 << LogIO::POST;
917 0 : auto n0 = col0.size();
918 0 : auto n1 = col1.size();
919 0 : auto n2 = col2.size();
920 0 : for (const auto& sample: *_samples) {
921 0 : ostringstream oss;
922 0 : oss << std::setw(n0) << sample.first << " " << std::setw(n1)
923 0 : << sample.second.first << " " << std::setw(n2)
924 0 : << sample.second.second;
925 0 : log << LogIO::NORMAL << oss.str() << LogIO::POST;
926 0 : }
927 0 : }
928 :
929 0 : void StatWtTVI::summarizeStats(Double& mean, Double& variance) const {
930 0 : LogIO log(LogOrigin("StatWtTVI", __func__));
931 0 : _logUsedChannels();
932 : try {
933 0 : mean = _wtStats->getStatistic(StatisticsData::MEAN);
934 0 : variance = _wtStats->getStatistic(StatisticsData::VARIANCE);
935 : log << LogIO::NORMAL << "The mean of the computed weights is "
936 0 : << mean << LogIO::POST;
937 : log << LogIO::NORMAL << "The variance of the computed weights is "
938 0 : << variance << LogIO::POST;
939 : log << LogIO::NORMAL << "Weights which had corresponding flags of True "
940 : << "prior to running this application were not used to compute these "
941 0 : << "stats." << LogIO::POST;
942 : }
943 0 : catch (const AipsError& x) {
944 : log << LogIO::WARN << "There was a problem calculating the mean and "
945 : << "variance of the weights computed by this application. Perhaps there "
946 : << "was something amiss with the input MS and/or the selection criteria. "
947 : << "Examples of such issues are that all the data were originally flagged "
948 : << "or that the sample size was consistently too small for computations "
949 0 : << "of variances" << LogIO::POST;
950 0 : setNaN(mean);
951 0 : setNaN(variance);
952 0 : }
953 0 : }
954 :
955 0 : void StatWtTVI::origin() {
956 : // Drive underlying ViImplementation2
957 0 : getVii()->origin();
958 : // Synchronize own VisBuffer
959 0 : configureNewSubchunk();
960 0 : _clearCache();
961 0 : }
962 :
963 0 : void StatWtTVI::next() {
964 : // Drive underlying ViImplementation2
965 0 : getVii()->next();
966 : // Synchronize own VisBuffer
967 0 : configureNewSubchunk();
968 0 : _clearCache();
969 0 : }
970 :
971 : }
972 :
973 : }
|
