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           0 : StatWt::StatWt(
      43             :     MeasurementSet* ms,
      44             :     const StatWtColConfig* const statwtColConfig
      45           0 : ) : _ms(ms),
      46           0 :     _saf(), _statwtColConfig(statwtColConfig) {
      47           0 :     ThrowIf(! _ms, "Input MS pointer cannot be NULL");
      48           0 :     ThrowIf(
      49             :         ! _statwtColConfig,
      50             :         "Input column configuration pointer cannot be NULL"
      51             :     );
      52           0 : }
      53             : 
      54           0 : StatWt::~StatWt() {}
      55             : 
      56           0 : void StatWt::setOutputMS(const casacore::String& outname) {
      57           0 :     _outname = outname;
      58           0 : }
      59             : 
      60           0 : void StatWt::setTimeBinWidth(const casacore::Quantity& binWidth) {
      61           0 :     _timeBinWidth = vi::StatWtTVI::getTimeBinWidthInSec(binWidth);
      62           0 : }
      63             : 
      64           0 : void StatWt::setTimeBinWidth(Double binWidth) {
      65           0 :     vi::StatWtTVI::checkTimeBinWidth(binWidth);
      66           0 :     _timeBinWidth = binWidth;
      67           0 : }
      68             : 
      69           0 : void StatWt::setCombine(const String& combine) {
      70           0 :     _combine = downcase(combine);
      71           0 : }
      72             : 
      73           0 : void StatWt::setPreview(casacore::Bool preview) {
      74           0 :     _preview = preview;
      75           0 : }
      76             : 
      77           0 : void StatWt::setTVIConfig(const Record& config) {
      78           0 :     _tviConfig = config;
      79           0 : }
      80             : 
      81           0 : Record StatWt::writeWeights() {
      82           0 :     auto mustWriteWt = False;
      83           0 :     auto mustWriteWtSp = False;
      84           0 :     auto mustWriteSig = False;
      85           0 :     auto mustWriteSigSp = False;
      86           0 :     _statwtColConfig->getColWriteFlags(
      87             :         mustWriteWt, mustWriteWtSp, mustWriteSig, mustWriteSigSp
      88             :     );
      89           0 :     std::shared_ptr<vi::VisibilityIterator2> vi;
      90           0 :     std::shared_ptr<vi::StatWtTVILayerFactory> factory;
      91           0 :     _constructVi(vi, factory);
      92           0 :     vi::VisBuffer2 *vb = vi->getVisBuffer();
      93           0 :     ProgressMeter pm(0, _ms->nrow(), "StatWt Progress");
      94           0 :     uInt64 count = 0;
      95           0 :     for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
      96           0 :         for (vi->origin(); vi->more(); vi->next()) {
      97           0 :             auto nrow = vb->nRows();
      98           0 :             if (_preview) {
      99             :                 // just need to run the flags to accumulate
     100             :                 // flagging info
     101           0 :                 vb->flagCube();
     102             :             }
     103             :             else {
     104           0 :                 if (mustWriteWtSp) {
     105           0 :                     auto& x = vb->weightSpectrum();
     106           0 :                     ThrowIf(
     107             :                         x.empty(),
     108             :                         "WEIGHT_SPECTRUM is only partially initialized. "
     109             :                         "StatWt cannot deal with such an MS"
     110             :                     );
     111           0 :                     vb->setWeightSpectrum(x);
     112             :                 }
     113           0 :                 if (mustWriteSigSp) {
     114           0 :                     auto& x = vb->sigmaSpectrum();
     115           0 :                     ThrowIf(
     116             :                         x.empty(),
     117             :                         "SIGMA_SPECTRUM is only partially initialized. "
     118             :                         "StatWt2 cannot deal with such an MS"
     119             :                     );
     120           0 :                     vb->setSigmaSpectrum(x);
     121             :                 }
     122           0 :                 if (mustWriteWt) {
     123           0 :                     vb->setWeight(vb->weight());
     124             :                 }
     125           0 :                 if (mustWriteSig) {
     126           0 :                     vb->setSigma(vb->sigma());
     127             :                 }
     128           0 :                 vb->setFlagCube(vb->flagCube());
     129           0 :                 vb->setFlagRow(vb->flagRow());
     130           0 :                 vb->writeChangesBack();
     131             :             }
     132           0 :             count += nrow;
     133           0 :             pm.update(count);
     134             :         }
     135             :     }
     136           0 :     if (_preview) {
     137           0 :         LogIO log(LogOrigin("StatWt", __func__));
     138             :         log << LogIO::NORMAL
     139             :             << "RAN IN PREVIEW MODE. NO WEIGHTS NOR FLAGS WERE CHANGED."
     140           0 :             << LogIO::POST;
     141           0 :     }
     142           0 :     factory->getTVI()->summarizeFlagging();
     143             :     Double mean, variance;
     144           0 :     factory->getTVI()->summarizeStats(mean, variance);
     145           0 :     Record ret;
     146           0 :     ret.define("mean", mean);
     147           0 :     ret.define("variance", variance);
     148           0 :     return ret;
     149           0 : }
     150             : 
     151           0 : 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           0 :     std::vector<Int> scs;
     161           0 :     scs.push_back(MS::ARRAY_ID);
     162           0 :     scs.push_back(MS::DATA_DESC_ID);
     163           0 :     scs.push_back(MS::TIME);
     164           0 :     if (! _combine.contains("scan")) {
     165           0 :         scs.push_back(MS::SCAN_NUMBER);
     166             :     }
     167           0 :     if (! _combine.contains("state")) {
     168           0 :         scs.push_back(MS::STATE_ID);
     169             :     }
     170           0 :     if (! _combine.contains("field")) {
     171           0 :         scs.push_back(MS::FIELD_ID);
     172             :     }
     173           0 :     Block<int> sort(scs.size());
     174           0 :     uInt i = 0;
     175           0 :     for (const auto& col: scs) {
     176           0 :         sort[i] = col;
     177           0 :         ++i;
     178             :     }
     179           0 :     vi::SortColumns sc(sort, False);
     180           0 :     vi::IteratingParameters ipar(_timeBinWidth, sc);
     181           0 :     vi::VisIterImpl2LayerFactory data(_ms, ipar, True);
     182           0 :     std::unique_ptr<Record> config(dynamic_cast<Record*>(_tviConfig.clone()));
     183           0 :     factory.reset(new vi::StatWtTVILayerFactory(*config));
     184           0 :     Vector<vi::ViiLayerFactory*> facts(2);
     185           0 :     facts[0] = &data;
     186           0 :     facts[1] = factory.get();
     187           0 :     vi.reset(new vi::VisibilityIterator2(facts));
     188           0 : }
     189             : 
     190             : }