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