Line data    Source code

       1             : //#
       2             : //#  CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
       3             : //#  Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All
       4             : //#  rights reserved.
       5             : //#  Copyright (C) European Southern Observatory, 2011, All rights reserved.
       6             : //#
       7             : //#  This library is free software; you can redistribute it and/or
       8             : //#  modify it under the terms of the GNU Lesser General Public
       9             : //#  License as published by the Free software Foundation; either
      10             : //#  version 2.1 of the License, or (at your option) any later version.
      11             : //#
      12             : //#  This library is distributed in the hope that it will be useful,
      13             : //#  but WITHOUT ANY WARRANTY, without even the implied warranty of
      14             : //#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
      15             : //#  Lesser General Public License for more details.
      16             : //#
      17             : //#  You should have received a copy of the GNU Lesser General Public
      18             : //#  License along with this library; if not, write to the Free Software
      19             : //#  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
      20             : //#  MA 02111-1307  USA
      21             : 
      22             : #include <mstransform/MSTransform/StatWt.h>
      23             : 
      24             : #include <casacore/casa/Containers/ValueHolder.h>
      25             : #include <casacore/casa/Quanta/QuantumHolder.h>
      26             : #include <casacore/casa/System/ProgressMeter.h>
      27             : #include <casacore/ms/MSOper/MSMetaData.h>
      28             : #include <casacore/tables/Tables/ArrColDesc.h>
      29             : #include <casacore/tables/Tables/TableProxy.h>
      30             : #include <casacore/tables/DataMan/TiledShapeStMan.h>
      31             : 
      32             : #include <mstransform/MSTransform/StatWtColConfig.h>
      33             : #include <mstransform/TVI/StatWtTVI.h>
      34             : #include <mstransform/TVI/StatWtTVILayerFactory.h>
      35             : #include <msvis/MSVis/ViImplementation2.h>
      36             : #include <msvis/MSVis/IteratingParameters.h>
      37             : 
      38             : using namespace casacore;
      39             : 
      40             : namespace casa { 
      41             : 
      42          47 : StatWt::StatWt(
      43             :     MeasurementSet* ms,
      44             :     const StatWtColConfig* const statwtColConfig
      45          47 : ) : _ms(ms),
      46          47 :     _saf(), _statwtColConfig(statwtColConfig) {
      47          47 :     ThrowIf(! _ms, "Input MS pointer cannot be NULL");
      48          47 :     ThrowIf(
      49             :         ! _statwtColConfig,
      50             :         "Input column configuration pointer cannot be NULL"
      51             :     );
      52          47 : }
      53             : 
      54          47 : StatWt::~StatWt() {}
      55             : 
      56           0 : void StatWt::setOutputMS(const casacore::String& outname) {
      57           0 :     _outname = outname;
      58           0 : }
      59             : 
      60          12 : void StatWt::setTimeBinWidth(const casacore::Quantity& binWidth) {
      61          12 :     _timeBinWidth = vi::StatWtTVI::getTimeBinWidthInSec(binWidth);
      62          12 : }
      63             : 
      64          35 : void StatWt::setTimeBinWidth(Double binWidth) {
      65          35 :     vi::StatWtTVI::checkTimeBinWidth(binWidth);
      66          35 :     _timeBinWidth = binWidth;
      67          35 : }
      68             : 
      69          47 : void StatWt::setCombine(const String& combine) {
      70          47 :     _combine = downcase(combine);
      71          47 : }
      72             : 
      73          47 : void StatWt::setPreview(casacore::Bool preview) {
      74          47 :     _preview = preview;
      75          47 : }
      76             : 
      77          47 : void StatWt::setTVIConfig(const Record& config) {
      78          47 :     _tviConfig = config;
      79          47 : }
      80             : 
      81          47 : Record StatWt::writeWeights() {
      82          47 :     auto mustWriteWt = False;
      83          47 :     auto mustWriteWtSp = False;
      84          47 :     auto mustWriteSig = False;
      85          47 :     auto mustWriteSigSp = False;
      86          47 :     _statwtColConfig->getColWriteFlags(
      87             :         mustWriteWt, mustWriteWtSp, mustWriteSig, mustWriteSigSp
      88             :     );
      89          47 :     std::shared_ptr<vi::VisibilityIterator2> vi;
      90          47 :     std::shared_ptr<vi::StatWtTVILayerFactory> factory;
      91          47 :     _constructVi(vi, factory);
      92          46 :     vi::VisBuffer2 *vb = vi->getVisBuffer();
      93          46 :     ProgressMeter pm(0, _ms->nrow(), "StatWt Progress");
      94          46 :     uInt64 count = 0;
      95         374 :     for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
      96        3887 :         for (vi->origin(); vi->more(); vi->next()) {
      97        3559 :             auto nrow = vb->nRows();
      98        3559 :             if (_preview) {
      99             :                 // just need to run the flags to accumulate
     100             :                 // flagging info
     101          60 :                 vb->flagCube();
     102             :             }
     103             :             else {
     104        3499 :                 if (mustWriteWtSp) {
     105        2723 :                     auto& x = vb->weightSpectrum();
     106        2723 :                     ThrowIf(
     107             :                         x.empty(),
     108             :                         "WEIGHT_SPECTRUM is only partially initialized. "
     109             :                         "StatWt cannot deal with such an MS"
     110             :                     );
     111        2723 :                     vb->setWeightSpectrum(x);
     112             :                 }
     113        3499 :                 if (mustWriteSigSp) {
     114         672 :                     auto& x = vb->sigmaSpectrum();
     115         672 :                     ThrowIf(
     116             :                         x.empty(),
     117             :                         "SIGMA_SPECTRUM is only partially initialized. "
     118             :                         "StatWt2 cannot deal with such an MS"
     119             :                     );
     120         672 :                     vb->setSigmaSpectrum(x);
     121             :                 }
     122        3499 :                 if (mustWriteWt) {
     123        3379 :                     vb->setWeight(vb->weight());
     124             :                 }
     125        3499 :                 if (mustWriteSig) {
     126        1328 :                     vb->setSigma(vb->sigma());
     127             :                 }
     128        3499 :                 vb->setFlagCube(vb->flagCube());
     129        3499 :                 vb->setFlagRow(vb->flagRow());
     130        3499 :                 vb->writeChangesBack();
     131             :             }
     132        3559 :             count += nrow;
     133        3559 :             pm.update(count);
     134             :         }
     135             :     }
     136          46 :     if (_preview) {
     137           2 :         LogIO log(LogOrigin("StatWt", __func__));
     138             :         log << LogIO::NORMAL
     139             :             << "RAN IN PREVIEW MODE. NO WEIGHTS NOR FLAGS WERE CHANGED."
     140           1 :             << LogIO::POST;
     141           1 :     }
     142          46 :     factory->getTVI()->summarizeFlagging();
     143             :     Double mean, variance;
     144          46 :     factory->getTVI()->summarizeStats(mean, variance);
     145          46 :     Record ret;
     146          46 :     ret.define("mean", mean);
     147          46 :     ret.define("variance", variance);
     148          92 :     return ret;
     149          48 : }
     150             : 
     151          47 : void StatWt::_constructVi(
     152             :     std::shared_ptr<vi::VisibilityIterator2>& vi,
     153             :     std::shared_ptr<vi::StatWtTVILayerFactory>& factory
     154             : ) const {
     155             :     // default sort columns are from MSIter and are ARRAY_ID, FIELD_ID,
     156             :     // DATA_DESC_ID, and TIME. I'm adding scan and state because, according to
     157             :     // the statwt requirements, by default, scan and state changes should mark
     158             :     // boundaries in the weights computation.
     159             :     // ORDER MATTERS. Columns are sorted in the order they appear in the vector.
     160          47 :     std::vector<Int> scs;
     161          47 :     scs.push_back(MS::ARRAY_ID);
     162          47 :     scs.push_back(MS::DATA_DESC_ID);
     163          47 :     scs.push_back(MS::TIME);
     164          47 :     if (! _combine.contains("scan")) {
     165          40 :         scs.push_back(MS::SCAN_NUMBER);
     166             :     }
     167          47 :     if (! _combine.contains("state")) {
     168          43 :         scs.push_back(MS::STATE_ID);
     169             :     }
     170          47 :     if (! _combine.contains("field")) {
     171          40 :         scs.push_back(MS::FIELD_ID);
     172             :     }
     173          47 :     Block<int> sort(scs.size());
     174          47 :     uInt i = 0;
     175         311 :     for (const auto& col: scs) {
     176         264 :         sort[i] = col;
     177         264 :         ++i;
     178             :     }
     179          47 :     vi::SortColumns sc(sort, False);
     180          47 :     vi::IteratingParameters ipar(_timeBinWidth, sc);
     181          47 :     vi::VisIterImpl2LayerFactory data(_ms, ipar, True);
     182          47 :     std::unique_ptr<Record> config(dynamic_cast<Record*>(_tviConfig.clone()));
     183          47 :     factory.reset(new vi::StatWtTVILayerFactory(*config));
     184          47 :     Vector<vi::ViiLayerFactory*> facts(2);
     185          47 :     facts[0] = &data;
     186          47 :     facts[1] = factory.get();
     187          47 :     vi.reset(new vi::VisibilityIterator2(facts));
     188          53 : }
     189             : 
     190             : }