casa / casa5 / bba3a34c967

Commits

Takeshi Nakazato authored bba3a34c967 Merge 18 May 2018

Merge pull request #517 in CASA/casa from bugfix/CAS-10530 to master

* commit '89759ec0108e36a3604148bd64dc25f07a065977':
  only initialize spectrum column if chanbin requires it if that column is uninitialized to start with
  change implementation so that _SPECTRUM columns are fully initialized even if an MS selection is used
  add ViiType() implementation to StatWtTVI
  refactoring
  enable initialization of sigma_spectrum if necessary, refactoring
  debug info added, committing so I can switch branches
  CAS-10530, yet more requirements implemented
  fix of no WEIGHT_SPECTRUM bugs

code/mstransform/MSTransform/StatWt.cc

Modified

-11
+//#  This library is distributed in the hope that it will be useful,
-12
+//#  but WITHOUT ANY WARRANTY, without even the implied warranty of
-13
+//#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-14
+//#  Lesser General Public License for more details.
-15
+//#
-16
+//#  You should have received a copy of the GNU Lesser General Public
-17
+//#  License along with this library; if not, write to the Free Software
-18
+//#  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
-19
+//#  MA 02111-1307  USA
-20
-+
+#include <mstransform/MSTransform/StatWt.h>
-+
-23
+#include <casacore/casa/Containers/ValueHolder.h>
-24
+#include <casacore/casa/Quanta/QuantumHolder.h>
-25
+#include <casacore/casa/System/ProgressMeter.h>
-26
+#include <casacore/ms/MSOper/MSMetaData.h>
-27
+#include <casacore/tables/Tables/ArrColDesc.h>
-28
+#include <casacore/tables/Tables/TableProxy.h>
-29
+#include <casacore/tables/DataMan/TiledShapeStMan.h>
-30
--
+#include <mstransform/MSTransform/StatWt.h>
-+
+#include <mstransform/MSTransform/StatWtColConfig.h>
-32
+#include <mstransform/TVI/StatWtTVI.h>
-33
+#include <mstransform/TVI/StatWtTVILayerFactory.h>
-34
+#include <msvis/MSVis/ViImplementation2.h>
-35
+#include <msvis/MSVis/IteratingParameters.h>
--
+#include <msvis/MSVis/LayeredVi2Factory.h>
-36
-37
+using namespace casacore;
-38
-39
+namespace casa {
-40
--
+StatWt::StatWt(MeasurementSet* ms) : _ms(ms), _saf() {
-+
+StatWt::StatWt(
-+
+    MeasurementSet* ms,
-+
+    const StatWtColConfig* const statwtColConfig
-+
+) : _ms(ms),
-+
+    _saf(), _statwtColConfig(statwtColConfig) {
-46
+    ThrowIf(! _ms, "Input MS pointer cannot be NULL");
-+
+    ThrowIf(
-+
+        ! _statwtColConfig,
-+
+        "Input column configuration pointer cannot be NULL"
-+
+    );
-51
+}
-52
-53
+StatWt::~StatWt() {}
-54
-55
+void StatWt::setOutputMS(const casacore::String& outname) {
-56
+    _outname = outname;
-57
+}
-58
-59
+void StatWt::setTimeBinWidth(const casacore::Quantity& binWidth) {
-60
+    _timeBinWidth = vi::StatWtTVI::getTimeBinWidthInSec(binWidth);
-78
+}
-79
-80
+void StatWt::setPreview(casacore::Bool preview) {
-81
+    _preview = preview;
-82
+}
-83
-84
+void StatWt::setTVIConfig(const Record& config) {
-85
+    _tviConfig = config;
-86
+}
-87
--
+Record StatWt::writeWeights() const {
--
+    auto hasWtSp = _ms->isColumn(MSMainEnums::WEIGHT_SPECTRUM);
--
+    auto mustWriteWtSp = ! _preview
--
+        && _tviConfig.isDefined(vi::StatWtTVI::CHANBIN);
--
+    if (mustWriteWtSp) {
--
+        auto type = _tviConfig.type(_tviConfig.fieldNumber(vi::StatWtTVI::CHANBIN));
--
+        if (type == TpArrayBool) {
--
+            // default variant type
--
+            mustWriteWtSp = False;
--
+        }
--
+        else if (type == TpString) {
--
+            auto val = _tviConfig.asString(vi::StatWtTVI::CHANBIN);
--
+            val.downcase();
--
+            if (val == "spw") {
--
+                mustWriteWtSp = False;
-+
+Record StatWt::writeWeights() {
-+
+    auto mustWriteWt = False;
-+
+    auto mustWriteWtSp = False;
-+
+    auto mustWriteSig = False;
-+
+    auto mustWriteSigSp = False;
-+
+    _statwtColConfig->getColWriteFlags(
-+
+        mustWriteWt, mustWriteWtSp, mustWriteSig, mustWriteSigSp
-+
+    );
-+
+    shared_ptr<vi::VisibilityIterator2> vi;
-+
+    std::shared_ptr<vi::StatWtTVILayerFactory> factory;
-+
+    _constructVi(vi, factory);
-+
+    vi::VisBuffer2 *vb = vi->getVisBuffer();
-+
+    ProgressMeter pm(0, _ms->nrow(), "StatWt Progress");
-+
+    uInt64 count = 0;
-+
+    for (vi->originChunks(); vi->moreChunks(); vi->nextChunk()) {
-+
+        for (vi->origin(); vi->more(); vi->next()) {
-+
+            auto nrow = vb->nRows();
-+
+            if (_preview) {
-+
+                // just need to run the flags to accumulate
-+
+                // flagging info
-+
+                vb->flagCube();
-109
+            }
--
+        }
--
+    }
--
+    auto mustInitWtSp = False;
--
+    if (! hasWtSp && mustWriteWtSp) {
--
+        // we must create WEIGHT_SPECTRUM
--
+        hasWtSp = True;
--
+        mustInitWtSp = True;
--
+        // from Calibrater.cc
--
+        // Nominal default tile shape
--
+        IPosition dts(3, 4, 32, 1024);
--
+        // Discern DATA's default tile shape and use it
--
+        const auto dminfo = _ms->dataManagerInfo();
--
+        for (uInt i=0; i<dminfo.nfields(); ++i) {
--
+            Record col = dminfo.asRecord(i);
--
+            if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
--
+                dts = IPosition(col.asRecord("SPEC").asArrayInt("DEFAULTTILESHAPE"));
--
+                break;
-+
+            else {
-+
+                if (mustWriteWtSp) {
-+
+                    auto& x = vb->weightSpectrum();
-+
+                    ThrowIf(
-+
+                        x.empty(),
-+
+                        "WEIGHT_SPECTRUM is only partially initialized. "
-+
+                        "StatWt2 cannot deal with such an MS"
-+
+                    );
-+
+                    vb->setWeightSpectrum(x);
-+
+                }
-+
+                if (mustWriteSigSp) {
-+
+                    auto& x = vb->sigmaSpectrum();
-+
+                    ThrowIf(
-+
+                        x.empty(),
-+
+                        "SIGMA_SPECTRUM is only partially initialized. "
-+
+                        "StatWt2 cannot deal with such an MS"
-+
+                    );
-+
+                    vb->setSigmaSpectrum(x);
-+
+                }
-+
+                if (mustWriteWt) {
-+
+                    vb->setWeight(vb->weight());
-+
+                }
-+
+                if (mustWriteSig) {
-+
+                    vb->setSigma(vb->sigma());
-+
+                }
-+
+                vb->setFlagCube(vb->flagCube());
-+
+                vb->setFlagRow(vb->flagRow());
-+
+                vb->writeChangesBack();
-138
+            }
-+
+            count += nrow;
-+
+            pm.update(count);
-141
+        }
--
+        // Add the column
--
+        String colWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
--
+        TableDesc tdWtSp;
--
+        tdWtSp.addColumn(ArrayColumnDesc<Float>(colWtSp, "weight spectrum", 2));
--
+        TiledShapeStMan wtSpStMan("TiledWgtSpectrum", dts);
--
+        _ms->addColumn(tdWtSp, wtSpStMan);
-142
+    }
--
+    else if (! _preview) {
--
+        // check to see if extant WEIGHT_SPECTRUM needs to be initialized
--
+        ArrayColumn<Float> col(*_ms, MS::columnName(MS::WEIGHT_SPECTRUM));
--
+        try {
--
+            col.get(0);
--
+            // its initialized, so even if we are using the full spw for
--
+            // binning, we still need to update WEIGHT_SPECTRUM
--
+            mustWriteWtSp = True;
--
+        }
--
+        catch (const AipsError& x) {
--
+            // its not initialized, so we aren't going to write to it unless
--
+            // chanbin has been specified to be less than the spw width
--
+            mustInitWtSp = mustWriteWtSp;
--
+        }
-+
+    if (_preview) {
-+
+        LogIO log(LogOrigin("StatWt", __func__));
-+
+        log << LogIO::NORMAL
-+
+            << "RAN IN PREVIEW MODE. NO WEIGHTS NOR FLAGS WERE CHANGED."
-+
+            << LogIO::POST;
-148
+    }
-+
+    factory->getTVI()->summarizeFlagging();
-+
+    Double mean, variance;
-+
+    factory->getTVI()->summarizeStats(mean, variance);
-+
+    Record ret;
-+
+    ret.define("mean", mean);
-+
+    ret.define("variance", variance);
-+
+    return ret;
-+
+}
-+
-+
+void StatWt::_constructVi(
-+
+    std::shared_ptr<vi::VisibilityIterator2>& vi,
-+
+    std::shared_ptr<vi::StatWtTVILayerFactory>& factory
-+
+) const {
-162
+    // default sort columns are from MSIter and are ARRAY_ID, FIELD_ID, DATA_DESC_ID, and TIME
-163
+    // I'm adding scan and state because, according to the statwt requirements, by default, scan
-164
+    // and state changes should mark boundaries in the weights computation
-165
+    std::vector<Int> scs;
-166
+    scs.push_back(MS::ARRAY_ID);
-167
+    if (! _combine.contains("scan")) {
-168
+        scs.push_back(MS::SCAN_NUMBER);
-169
+    }
-170
+    if (! _combine.contains("state")) {
-171
+        scs.push_back(MS::STATE_ID);
-178
+    Block<int> sort(scs.size());
-179
+    uInt i = 0;
-180
+    for (const auto& col: scs) {
-181
+        sort[i] = col;
-182
+        ++i;
-183
+    }
-184
+    vi::SortColumns sc(sort, False);
-185
+    vi::IteratingParameters ipar(_timeBinWidth, sc);
-186
+    vi::VisIterImpl2LayerFactory data(_ms, ipar, True);
-187
+    unique_ptr<Record> config(dynamic_cast<Record*>(_tviConfig.clone()));
--
+    vi::StatWtTVILayerFactory statWtLayerFactory(*config);
-+
+    factory.reset(new vi::StatWtTVILayerFactory(*config));
-189
+    Vector<vi::ViiLayerFactory*> facts(2);
-190
+    facts[0] = &data;
--
+    facts[1] = &statWtLayerFactory;
--
+    vi::VisibilityIterator2 vi(facts);
--
+    vi::VisBuffer2 *vb = vi.getVisBuffer();
--
+    Vector<Int> vr(1);
--
+    ProgressMeter pm(0, _ms->nrow(), "StatWt Progress");
--
+    uInt64 count = 0;
--
+    for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
--
+        for (vi.origin(); vi.more(); vi.next()) {
--
+            auto nrow = vb->nRows();
--
+            if (_preview) {
--
+                // just need to run the flags to accumulate
--
+                // flagging info
--
+                vb->flagCube();
--
+            }
--
+            else {
--
+                if (mustInitWtSp) {
--
+                    auto nchan = vb->nChannels();
--
+                    auto ncor = vb->nCorrelations();
--
+                    Cube<Float> newwtsp(ncor, nchan, nrow, 0);
--
+                    vb->initWeightSpectrum(newwtsp);
--
+                    vb->writeChangesBack();
--
+                }
--
+                if (mustWriteWtSp) {
--
+                    vb->setWeightSpectrum(vb->weightSpectrum());
--
+                }
--
+                vb->setWeight(vb->weight());
--
+                vb->setFlagCube(vb->flagCube());
--
+                vb->setFlagRow(vb->flagRow());
--
+                vb->writeChangesBack();
--
+            }
--
+            count += nrow;
--
+            pm.update(count);
--
+        }
--
+    }
--
+    if (_preview) {
--
+        LogIO log(LogOrigin("StatWt", __func__));
--
+        log << LogIO::NORMAL
--
+            << "RAN IN PREVIEW MODE. NO WEIGHTS NOR FLAGS WERE CHANGED."
--
+            << LogIO::POST;
--
+    }
--
+    statWtLayerFactory.getTVI()->summarizeFlagging();
--
+    Double mean, variance;
--
+    statWtLayerFactory.getTVI()->summarizeStats(mean, variance);
--
+    Record ret;
--
+    ret.define("mean", mean);
--
+    ret.define("variance", variance);
--
+    return ret;
-+
+    facts[1] = factory.get();
-+
+    vi.reset(new vi::VisibilityIterator2(facts));
-193
+}
-194
-195
+}
--

casa5

Commits

Jira Issues

Add shortcut