casa5
Public
Commits
| 66 66 | auto q = qh.asQuantity(); |
| 67 67 | _setChanBinMap(qh.asQuantity()); |
| 68 68 | break; |
| 69 69 | } |
| 70 70 | default: |
| 71 71 | ThrowCc("Unsupported data type for chanbin"); |
| 72 72 | } |
| 73 73 | } |
| 74 74 | else { |
| 75 75 | _setDefaultChanBinMap(); |
| 76 + | } |
| 77 + | if (config.isDefined("minsamp")) { |
| 78 + | config.get("minsamp", _minSamp); |
| 76 79 | } |
| 77 80 | return True; |
| 78 81 | } |
| 79 82 | |
| 80 83 | void StatWtTVI::_setChanBinMap(const Quantity& binWidth) { |
| 81 84 | if (! binWidth.isConform(Unit("Hz"))) { |
| 82 85 | ostringstream oss; |
| 83 86 | oss << "If specified as a quantity, chanbin must have frequency units. " |
| 84 87 | << binWidth << " does not."; |
| 85 88 | ThrowCc(oss.str()); |
| 259 262 | blcb.spw = spw; |
| 260 263 | auto bins = _chanBins.find(spw)->second; |
| 261 264 | auto biter = bins.begin(); |
| 262 265 | auto bend = bins.end(); |
| 263 266 | blc[2] = i; |
| 264 267 | trc[2] = i; |
| 265 268 | for (; biter!=bend; ++biter) { |
| 266 269 | blc[1] = biter->start; |
| 267 270 | trc[1] = biter->end; |
| 268 271 | blcb.chanBin = *biter; |
| 269 - | flagCube(blc, trc) = _weights.find(blcb)->second == 0; |
| 272 + | if (_weights.find(blcb)->second == 0) { |
| 273 + | flagCube(blc, trc) = True; |
| 274 + | } |
| 270 275 | } |
| 271 276 | } |
| 272 277 | _newFlag = flagCube.copy(); |
| 273 278 | } |
| 274 279 | |
| 275 280 | void StatWtTVI::flagRow (Vector<Bool>& flagRow) const { |
| 276 281 | ThrowIf(! _weightsComputed, "Weights have not been computed yet"); |
| 277 282 | if (! _newFlagRow.empty()) { |
| 278 283 | flagRow = _newFlagRow.copy(); |
| 279 284 | return; |
| 378 383 | data[blcb](blc, trc) = dataCube(dataCubeBLC, dataCubeTRC); |
| 379 384 | flags[blcb](blc, trc) = flagCube(dataCubeBLC, dataCubeTRC); |
| 380 385 | } |
| 381 386 | } |
| 382 387 | } |
| 383 388 | } |
| 384 389 | // data has been gathered, now compute weights |
| 385 390 | _computeWeights(data, flags); |
| 386 391 | } |
| 387 392 | |
| 388 - | void StatWtTVI::writeBackChanges(VisBuffer2 * vb) { |
| 393 + | void StatWtTVI::writeBackChanges(VisBuffer2 *vb) { |
| 389 394 | // Pass to next layer down |
| 390 395 | getVii()->writeBackChanges(vb); |
| 391 396 | } |
| 392 397 | |
| 393 398 | StatWtTVI::Baseline StatWtTVI::_baseline(uInt ant1, uInt ant2) { |
| 394 399 | Baseline baseline; |
| 395 400 | if (ant1 < ant2) { |
| 396 401 | // this may always be the case, but I'm not certain, |
| 397 402 | baseline.first = ant1; |
| 398 403 | baseline.second = ant2; |
| 413 418 | stats.insert(StatisticsData::VARIANCE); |
| 414 419 | csReal.setStatsToCalculate(stats); |
| 415 420 | csImag.setStatsToCalculate(stats); |
| 416 421 | auto diter = data.begin(); |
| 417 422 | auto dend = data.end(); |
| 418 423 | auto fiter = flags.begin(); |
| 419 424 | for (; diter!=dend; ++diter, ++fiter) { |
| 420 425 | auto blcb = diter->first; |
| 421 426 | auto dataForBLCB = diter->second; |
| 422 427 | const auto npts = dataForBLCB.size(); |
| 423 - | if (npts == 1) { |
| 424 - | // one data point, trivial |
| 428 + | if (npts < _minSamp) { |
| 429 + | // not enough points, trivial |
| 425 430 | _weights[blcb] = 0; |
| 426 431 | } |
| 427 432 | else { |
| 428 433 | auto flagsForBLCB = fiter->second; |
| 429 - | if (allTrue(flagsForBLCB)) { |
| 430 - | // all data flagged, trivial |
| 434 + | if (nfalse(flagsForBLCB) < _minSamp) { |
| 435 + | // not enough points, trivial |
| 431 436 | _weights[blcb] = 0; |
| 432 437 | } |
| 433 438 | else { |
| 434 439 | // some data not flagged |
| 435 440 | const auto realPart = real(dataForBLCB); |
| 436 441 | const auto imagPart = imag(dataForBLCB); |
| 437 442 | const auto mask = ! flagsForBLCB; |
| 438 443 | const auto riter = realPart.begin(); |
| 439 444 | const auto iiter = imagPart.begin(); |
| 440 445 | const auto miter = mask.begin(); |