[casacore] r21572 committed - Merged in CASA changes

[casa...@googlecode.com]

Revision: 21572
Author: gervandiepen
Date: Tue Mar 3 12:22:11 2015 UTC
Log: Merged in CASA changes

https://code.google.com/p/casacore/source/detail?r=21572

Added:
/branches/nov14/lattices/LatticeMath/StatsTiledCollapser.h
/branches/nov14/lattices/LatticeMath/StatsTiledCollapser.tcc
/branches/nov14/ms/MeasurementSets/MSKeys.cc
/branches/nov14/ms/MeasurementSets/MSKeys.h
Modified:
/branches/nov14/CMakeLists.txt
/branches/nov14/casa/Arrays/Slicer.h
/branches/nov14/casa/BasicMath/Math.h
/branches/nov14/casa/Logging/LogOrigin.cc
/branches/nov14/casa/Logging/LogOrigin.h
/branches/nov14/casa/Utilities/CountedPtr.h
/branches/nov14/coordinates/Coordinates/DirectionCoordinate.cc
/branches/nov14/fits/FITS/blockio.cc
/branches/nov14/lattices/CMakeLists.txt
/branches/nov14/lattices/LatticeMath/LatticeStatistics.h
/branches/nov14/lattices/LatticeMath/LatticeStatistics.tcc
/branches/nov14/lattices/LatticeMath/LatticeStatsDataProvider.h
/branches/nov14/lattices/LatticeMath/LatticeStatsDataProvider.tcc
/branches/nov14/lattices/LatticeMath/LatticeStatsDataProviderBase.tcc
/branches/nov14/lattices/LatticeMath/MaskedLatticeStatsDataProvider.h
/branches/nov14/lattices/LatticeMath/MaskedLatticeStatsDataProvider.tcc
/branches/nov14/lattices/LatticeMath/test/tLatticeStatistics.cc
/branches/nov14/lattices/Lattices/SubLattice.h
/branches/nov14/mainpage.dox
/branches/nov14/measures/Measures/MDirection.cc
/branches/nov14/measures/Measures/MDirection.h
/branches/nov14/measures/Measures/MeasTableMul.h
/branches/nov14/ms/CMakeLists.txt
/branches/nov14/ms/MeasurementSets/MSMetaData.cc
/branches/nov14/ms/MeasurementSets/MSMetaData.h
/branches/nov14/ms/MeasurementSets/MSSummary.cc
/branches/nov14/ms/MeasurementSets/test/tMSMetaData.cc
/branches/nov14/scimath/Mathematics/ClassicalStatistics.h
/branches/nov14/scimath/Mathematics/ClassicalStatistics.tcc
/branches/nov14/scimath/Mathematics/FitToHalfStatistics.tcc
/branches/nov14/scimath/Mathematics/StatisticsAlgorithm.h
/branches/nov14/scimath/Mathematics/StatisticsAlgorithm.tcc

=======================================
--- /dev/null
+++ /branches/nov14/lattices/LatticeMath/StatsTiledCollapser.h Tue Mar 3
12:22:11 2015 UTC
@@ -0,0 +1,174 @@
+//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
+//# Associated Universities, Inc. Washington DC, USA.
+//#
+//# This library is free software; you can redistribute it and/or modify it
+//# under the terms of the GNU Library General Public License as published
by
+//# the Free Software Foundation; either version 2 of the License, or (at
your
+//# option) any later version.
+//#
+//# This library is distributed in the hope that it will be useful, but
WITHOUT
+//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
+//# License for more details.
+//#
+//# You should have received a copy of the GNU Library General Public
License
+//# along with this library; if not, write to the Free Software Foundation,
+//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
+//#
+//# Correspondence concerning AIPS++ should be addressed as follows:
+//# Internet email: aips2-...@nrao.edu.
+//# Postal address: AIPS++ Project Office
+//# National Radio Astronomy Observatory
+//# 520 Edgemont Road
+//# Charlottesville, VA 22903-2475 USA
+//#
+//# $Id: LatticeStatistics.h 20739 2009-09-29 01:15:15Z Malte.Marquarding $
+
+#ifndef LATTICES_STATSTILEDCOLLAPSER_H
+#define LATTICES_STATSTILEDCOLLAPSER_H
+
+
+//# Includes
+#include <casacore/casa/aips.h>
+
+namespace casacore {
+
+
+// <summary> Generate statistics, tile by tile, from a masked lattice
</summary>
+
+// NOTE this version was moved from LatticeStatistics (early Dec 2014
version)
+// and slightly modified mostly for style issues (no significant semantic
differences
+// from that version). For a large number of statistics sets that need to
be computed
+// simultaneously, this version is more efficient than using the new stats
framework,
+// because creating large numbers of eg ClassicalStatistics objects is
much less efficient
+// than the direct manipulation of pointers to primitive types that this
class does.
+//
+// <use visibility=export>
+//
+// <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
+// </reviewed>
+//
+// <prerequisite>
+// <li> <linkto class=LatticeApply>LatticeApply</linkto>
+// <li> <linkto class=TiledCollapser>TiledCollapser</linkto>
+// </prerequisite>
+//
+// <etymology>
+// This class is used by <src>LatticeStatistics</src> to generate
+// statistical sum from an input <src>MaskedLattice</src>.
+// The input lattice is iterated through in tile-sized chunks
+// and fed to an object of this class.
+// </etymology>
+//
+// <synopsis>
+// <src>StatsTiledCollapser</src> is derived from
<src>TiledCollapser</src> which
+// is a base class used to define methods. Objects of this base class are
+// used by <src>LatticeApply</src> functions. In this particular case,
+// we are interested in <src>LatticeApply::tiledApply</src>. This
function iterates
+// through a <src>MaskedLattice</src> and allows you to collapse one or
more
+// axes, computing some values from it, and placing those values into
+// an output <src>MaskedLattice</src>. It iterates through the input
+// lattice in optimal tile-sized chunks. <src>LatticeStatistics</src>
+// uses a <src>StatsTiledCollapser</src> object which it gives to
+// <src>LatticeApply::tiledApply</src> for digestion. After it has
+// done its work, <src>LatticeStatistics</src> then accesses the output
+// <src>Lattice</src> that it made.
+// </synopsis>
+//
+// <example>
+// <srcblock>
+//// Create collapser. Control information is passed in via the constructor
+//
+// StatsTiledCollapser<T> collapser(range_p, noInclude_p, noExclude_p,
+// fixedMinMax_p, blcParent_p);
+//
+//// This is the first output axis getting collapsed values. In
LatticeStatistics
+//// this is the last axis of the output lattice
+//
+// Int newOutAxis = outLattice.ndim()-1;
+//
+//// tiledApply does the work by passing the collapser data in chunks
+//// and by writing the results into the output lattice
+//
+// LatticeApply<T>::tiledApply(outLattice, inLattice,
+// collapser, collapseAxes,
+// newOutAxis);
+//
+// </srcblock>
+// In this example, a collapser is made and passed to LatticeApply.
+// Afterwards, the output Lattice is available for use.
+// The Lattices must all be the correct shapes on input to tiledApply
+// </example>
+//
+// <motivation>
+// The LatticeApply classes enable the ugly details of optimal
+// Lattice iteration to be hidden from the user.
+// </motivation>
+//
+// <todo asof="1998/05/10">
+// <li>
+// </todo>
+
+template <class T, class U=T>
+class StatsTiledCollapser : public TiledCollapser<T, U> {
+public:
+ // Constructor provides pixel selection range and whether that
+ // range is an inclusion or exclusion range. If
<src>fixedMinMax=True</src>
+ // and an inclusion range is given, the min and max is set to
+ // that inclusion range.
+ StatsTiledCollapser(
+ const Vector<T>& pixelRange, Bool noInclude,
+ Bool noExclude, Bool fixedMinMax
+ );
+
+ virtual ~StatsTiledCollapser() {}
+
+ // Initialize process, making some checks
+ virtual void init (uInt nOutPixelsPerCollapse);
+
+ // Initialiaze the accumulator
+ virtual void initAccumulator (uInt n1, uInt n3);
+
+ // Process the data in the current chunk.
+ virtual void process (
+ uInt accumIndex1, uInt accumIndex3,
+ const T* inData, const Bool* inMask,
+ uInt dataIncr, uInt maskIncr,
+ uInt nrval, const IPosition& startPos,
+ const IPosition& shape
+ );
+
+ // End the accumulation process and return the result arrays
+ virtual void endAccumulator(Array<U>& result,
+ Array<Bool>& resultMask,
+ const IPosition& shape);
+
+ // Can handle null mask
+ virtual Bool canHandleNullMask() const {return True;};
+
+ // Find the location of the minimum and maximum data values
+ // in the input lattice.
+ void minMaxPos(IPosition& minPos, IPosition& maxPos);
+
+private:
+ Vector<T> _range;
+ Bool _include, _exclude, _fixedMinMax, _isReal;
+ IPosition _minpos, _maxpos;
+
+ // Accumulators for sum, sum squared, number of points
+ // minimum, and maximum
+
+ CountedPtr<Block<U> > _sum, _sumSq, _npts,
+ _mean, _variance, _nvariance;
+ CountedPtr<Block<T> > _min, _max;
+ CountedPtr<Block<Bool> > _initMinMax;
+
+ uInt _n1, _n3;
+};
+
+}
+
+#ifndef CASACORE_NO_AUTO_TEMPLATES
+#include <casacore/lattices/LatticeMath/StatsTiledCollapser.tcc>
+#endif
+#endif
=======================================
--- /dev/null
+++ /branches/nov14/lattices/LatticeMath/StatsTiledCollapser.tcc Tue Mar 3
12:22:11 2015 UTC
@@ -0,0 +1,309 @@
+//# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2004
+//# Associated Universities, Inc. Washington DC, USA.
+//#
+//# This library is free software; you can redistribute it and/or modify it
+//# under the terms of the GNU Library General Public License as published
by
+//# the Free Software Foundation; either version 2 of the License, or (at
your
+//# option) any later version.
+//#
+//# This library is distributed in the hope that it will be useful, but
WITHOUT
+//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
+//# License for more details.
+//#
+//# You should have received a copy of the GNU Library General Public
License
+//# along with this library; if not, write to the Free Software Foundation,
+//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
+//#
+//# Correspondence concerning AIPS++ should be addressed as follows:
+//# Internet email: aips2-...@nrao.edu.
+//# Postal address: AIPS++ Project Office
+//# National Radio Astronomy Observatory
+//# 520 Edgemont Road
+//# Charlottesville, VA 22903-2475 USA
+//#
+//# $Id: LatticeStatistics.tcc 20652 2009-07-06 05:04:32Z
Malte.Marquarding $
+
+#include <casacore/lattices/LatticeMath/StatsTiledCollapser.h>
+
+namespace casacore {
+
+template <class T, class U>
+StatsTiledCollapser<T,U>::StatsTiledCollapser(
+ const Vector<T>& pixelRange,
+ Bool noInclude, Bool noExclude,
+ Bool fixedMinMax
+) : _range(pixelRange), _include(! noInclude),
+ _exclude(! noExclude), _fixedMinMax(fixedMinMax),
+ _isReal(isReal(whatType(&*(CountedPtr<T>(new T(0)))))),
+ _minpos(0), _maxpos(0) {}
+
+template <class T, class U>
+void StatsTiledCollapser<T,U>::init (uInt nOutPixelsPerCollapse) {
+ AlwaysAssert (nOutPixelsPerCollapse == LatticeStatsBase::NACCUM,
AipsError);
+}
+
+template <class T, class U>
+void StatsTiledCollapser<T,U>::initAccumulator (uInt n1, uInt n3) {
+ _sum = new Block<U>(n1*n3);
+ _sumSq = new Block<U>(n1*n3);
+ _npts = new Block<U>(n1*n3);
+ _mean = new Block<U>(n1*n3);
+ _variance = new Block<U>(n1*n3);
+ _nvariance = new Block<U>(n1*n3);
+
+ _min = new Block<T>(n1*n3);
+ _max = new Block<T>(n1*n3);
+ _initMinMax = new Block<Bool>(n1*n3);
+ _sum->set(0);
+ _sumSq->set(0);
+ _npts->set(0);
+ _mean->set(0);
+ _variance->set(0);
+ _nvariance->set(0);
+
+ _min->set(0);
+ _max->set(0);
+ _initMinMax->set(True);
+ _n1 = n1;
+ _n3 = n3;
+}
+
+template <class T, class U>
+void StatsTiledCollapser<T,U>::process (
+ uInt index1, uInt index3,
+ const T* pInData, const Bool* pInMask,
+ uInt dataIncr, uInt maskIncr,
+ uInt nrval, const IPosition& startPos,
+ const IPosition& shape
+) {
+ // Process the data in the current chunk. Everything in this
+ // chunk belongs in one output location in the storage
+ // lattices
+ uInt index = index1 + index3*_n1;
+ U& sum = (*_sum)[index];
+ U& sumSq = (*_sumSq)[index];
+ U& nPts = (*_npts)[index];
+ T& dataMin = (*_min)[index];
+ T& dataMax = (*_max)[index];
+ U& mean = (*_mean)[index];
+ U& variance = (*_variance)[index];
+ U& nvariance = (*_nvariance)[index];
+ Bool& minMaxInit = (*_initMinMax)[index];
+
+ // If these are != -1 after the accumulating, then
+ // the min and max were updated
+
+ Int minLoc = -1;
+ Int maxLoc = -1;
+ //
+ T useIt;
+ if (pInMask == 0) {
+ // All pixels are good
+ if (_include) {
+ T datum;
+ for (uInt i=0; i<nrval; i++) {
+ datum = *pInData;
+ useIt = LattStatsSpecialize::usePixelInc (_range(0), _range(1), datum);
+ LattStatsSpecialize::accumulate(
+ nPts, sum, mean, nvariance, variance, sumSq,
+ dataMin, dataMax, minLoc, maxLoc, minMaxInit,
+ False, *pInData, i, useIt
+ );
+ pInData += dataIncr;
+ }
+ if (_fixedMinMax) {
+ dataMin = _range(0);
+ dataMax = _range(1);
+ }
+ } else if (_exclude) {
+ T datum;
+ for (uInt i=0; i<nrval; i++) {
+ datum = *pInData;
+ useIt = LattStatsSpecialize::usePixelExc (
+ _range(0), _range(1), datum
+ );
+ LattStatsSpecialize::accumulate(
+ nPts, sum, mean, nvariance, variance, sumSq,
+ dataMin, dataMax, minLoc, maxLoc, minMaxInit,
+ False, *pInData, i, useIt
+ );
+ pInData += dataIncr;
+ }
+ }
+ else {
+ // All data accepted
+ LattStatsSpecialize::setUseItTrue(useIt);
+ for (uInt i=0; i<nrval; i++) {
+ LattStatsSpecialize::accumulate(
+ nPts, sum, mean, nvariance, variance, sumSq,
+ dataMin, dataMax, minLoc, maxLoc, minMaxInit,
+ False, *pInData, i, useIt
+ );
+ pInData += dataIncr;
+ }
+ }
+ }
+ else {
+ // Some pixels are masked
+ if (_include) {
+ T datum;
+ Bool mask;
+ for (uInt i=0; i<nrval; i++) {
+ datum = *pInData;
+ mask = *pInMask;
+ if (mask) {
+ useIt = LattStatsSpecialize::usePixelInc (
+ _range(0), _range(1), datum
+ );
+ LattStatsSpecialize::accumulate(
+ nPts, sum, mean, nvariance, variance, sumSq,
+ dataMin, dataMax, minLoc, maxLoc, minMaxInit,
+ False, *pInData, i, useIt
+ );
+ }
+ pInData += dataIncr;
+ pInMask += maskIncr;
+ }
+ if (_fixedMinMax) {
+ dataMin = _range(0);
+ dataMax = _range(1);
+ }
+ }
+ else if (_exclude) {
+ T datum;
+ Bool mask;
+ for (uInt i=0; i<nrval; i++) {
+ datum = *pInData;
+ mask = *pInMask;
+ if (mask) {
+ useIt = LattStatsSpecialize::usePixelExc (_range(0), _range(1),
datum);
+ LattStatsSpecialize::accumulate(
+ nPts, sum, mean, nvariance, variance, sumSq,
+ dataMin, dataMax, minLoc, maxLoc, minMaxInit,
+ False, *pInData, i, useIt
+ );
+ }
+ pInData += dataIncr;
+ pInMask += maskIncr;
+ }
+ }
+ else {
+ // All data accepted
+ LattStatsSpecialize::setUseItTrue(useIt);
+ for (uInt i=0; i<nrval; i++) {
+ if (*pInMask) {
+ LattStatsSpecialize::accumulate(
+ nPts, sum, mean, nvariance, variance, sumSq,
+ dataMin, dataMax, minLoc, maxLoc, minMaxInit,
+ False, *pInData, i, useIt
+ );
+ }
+ pInData += dataIncr;
+ pInMask += maskIncr;
+ }
+ }
+ }
+
+ // Update overall min and max location. These are never updated
+ // if fixedMinMax is true. These values are only meaningful for
+ // Float images. For Complex they are useless currently.
+
+ //DataType type = whatType(*T(0));
+ if (_isReal) {
+ if (minLoc != -1) {
+ _minpos = startPos + toIPositionInArray(minLoc, shape);
+ }
+ if (maxLoc != -1) {
+ _maxpos = startPos + toIPositionInArray(maxLoc, shape);
+ }
+ }
+}
+
+template <class T, class U>
+void StatsTiledCollapser<T,U>::endAccumulator(
+ Array<U>& result, Array<Bool>& resultMask,
+ const IPosition& shape
+) {
+
+ // Reshape arrays. The mask is always true. Any locations
+ // in the storage lattice for which there were no valid points
+ // will have the NPTS field set to zero. That is what
+ // we use to effectively mask it.
+
+ result.resize(shape);
+ resultMask.resize(shape);
+ resultMask.set(True);
+
+ Bool deleteRes;
+ U* res = result.getStorage (deleteRes);
+ U* resptr = res;
+
+ U* sumPtr = _sum->storage();
+ U* sumSqPtr = _sumSq->storage();
+ U* nPtsPtr = _npts->storage();
+ U* meanPtr = _mean->storage();
+ U* variancePtr = _variance->storage();
+ const T* minPtr = _min->storage();
+ const T* maxPtr = _max->storage();
+
+ uInt i,j;
+ U* resptr_root = resptr;
+ for (i=0; i<_n3; i++) {
+ resptr = resptr_root + (Int(LatticeStatsBase::NPTS) * _n1);
+ objcopy (resptr, nPtsPtr, _n1);
+ nPtsPtr += _n1;
+
+ resptr = resptr_root + (Int(LatticeStatsBase::SUM) * _n1);
+ objcopy (resptr, sumPtr, _n1);
+ sumPtr += _n1;
+
+ resptr = resptr_root + (Int(LatticeStatsBase::SUMSQ) * _n1);
+ objcopy (resptr, sumSqPtr, _n1);
+ sumSqPtr += _n1;
+
+ resptr = resptr_root + (Int(LatticeStatsBase::MEAN) * _n1);
+ objcopy (resptr, meanPtr, _n1);
+ meanPtr += _n1;
+
+ resptr = resptr_root + (Int(LatticeStatsBase::VARIANCE) * _n1);
+ objcopy (resptr, variancePtr, _n1);
+ variancePtr += _n1;
+
+ resptr = resptr_root + (Int(LatticeStatsBase::MIN) * _n1);
+ for (j=0; j<_n1; j++) {
+ convertScalar (*resptr++, *minPtr++);
+ }
+
+ resptr = resptr_root + (Int(LatticeStatsBase::MAX) * _n1);
+ for (j=0; j<_n1; j++) {
+ convertScalar (*resptr++, *maxPtr++);
+ }
+
+ resptr_root += _n1 * Int(LatticeStatsBase::NACCUM);
+ }
+
+ _sum = NULL;
+ _sumSq = NULL;
+ _npts = NULL;
+ _min = NULL;
+ _max = NULL;
+ _initMinMax = NULL;
+ _mean = NULL;
+ _variance = NULL;
+ _nvariance = NULL;
+
+ result.putStorage (res, deleteRes);
+}
+
+template <class T, class U>
+void StatsTiledCollapser<T,U>::minMaxPos(IPosition& minPos, IPosition&
maxPos)
+{
+ minPos.resize(_minpos.nelements());
+ minPos = _minpos;
+ maxPos.resize(_maxpos.nelements());
+ maxPos = _maxpos;
+}
+
+}
+
=======================================
--- /dev/null
+++ /branches/nov14/ms/MeasurementSets/MSKeys.cc Tue Mar 3 12:22:11 2015
UTC
@@ -0,0 +1,137 @@
+//# Copyright (C) 1998,1999,2000,2001
+//# Associated Universities, Inc. Washington DC, USA.
+//#
+//# This library is free software; you can redistribute it and/or modify it
+//# under the terms of the GNU Library General Public License as published
by
+//# the Free Software Foundation; either version 2 of the License, or (at
your
+//# option) any later version.
+//#
+//# This library is distributed in the hope that it will be useful, but
WITHOUT
+//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
+//# License for more details.
+//#
+//# You should have received a copy of the GNU Library General Public
License
+//# along with this library; if not, write to the Free Software Foundation,
+//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
+//#
+//# Correspondence concerning AIPS++ should be addressed as follows:
+//# Internet email: aips2-...@nrao.edu.
+//# Postal address: AIPS++ Project Office
+//# National Radio Astronomy Observatory
+//# 520 Edgemont Road
+//# Charlottesville, VA 22903-2475 USA
+//#
+
+#include <casacore/ms/MeasurementSets/MSKeys.h>
+
+#include <casacore/casa/BasicSL/String.h>
+
+namespace casacore {
+
+Bool operator<(const SubScanKey& lhs, const SubScanKey& rhs) {
+ if (lhs.obsID < rhs.obsID) {
+ return True;
+ }
+ else if (lhs.obsID == rhs.obsID) {
+ if (lhs.arrayID < rhs.arrayID) {
+ return True;
+ }
+ else if (lhs.arrayID == rhs.arrayID) {
+ if (lhs.scan < rhs.scan) {
+ return True;
+ }
+ else if (lhs.scan == rhs.scan) {
+ if (lhs.fieldID < rhs.fieldID) {
+ return True;
+ }
+ }
+ }
+ }
+ return False;
+}
+
+String toString(const SubScanKey& subScanKey) {
+ return toString(scanKey(subScanKey)) + " fieldID="
+ + String::toString(subScanKey.fieldID);
+}
+
+std::ostream& operator<<(std::ostream& os, const SubScanKey& subScanKey) {
+ os << toString(subScanKey) << endl;
+ return os;
+}
+
+
+String toString(const ScanKey& scanKey) {
+ return "observationID=" + String::toString(scanKey.obsID)
+ + " arrayID=" + String::toString(scanKey.arrayID)
+ + " scan number=" + String::toString(scanKey.scan);
+}
+
+
+Bool operator<(const ScanKey& lhs, const ScanKey& rhs) {
+ if (lhs.obsID < rhs.obsID) {
+ return True;
+ }
+ else if (lhs.obsID == rhs.obsID) {
+ if (lhs.arrayID < rhs.arrayID) {
+ return True;
+ }
+ else if (lhs.arrayID == rhs.arrayID) {
+ if (lhs.scan < rhs.scan) {
+ return True;
+ }
+ }
+ }
+ return False;
+}
+
+std::set<Int> scanNumbers(const std::set<ScanKey>& scanKeys) {
+ std::set<Int> scanNumbers;
+ std::set<ScanKey>::const_iterator iter = scanKeys.begin();
+ std::set<ScanKey>::const_iterator end = scanKeys.end();
+ while (iter != end) {
+ scanNumbers.insert(iter->scan);
+ ++iter;
+ }
+ return scanNumbers;
+}
+
+ostream& operator<<(ostream& os, const ScanKey& scanKey) {
+ os << toString(scanKey) << endl;
+ return os;
+}
+
+
+Bool operator<(const ArrayKey& lhs, const ArrayKey& rhs) {
+ if (lhs.obsID < rhs.obsID) {
+ return True;
+ }
+ else if (lhs.obsID == rhs.obsID) {
+ if (lhs.arrayID < rhs.arrayID) {
+ return True;
+ }
+ }
+ return False;
+}
+
+std::set<ScanKey> scanKeys(
+ const std::set<Int>& scans, const ArrayKey& arrayKey
+) {
+ std::set<ScanKey> scanKeys;
+ std::set<Int>::const_iterator iter = scans.begin();
+ std::set<Int>::const_iterator end = scans.end();
+ ScanKey scanKey;
+ scanKey.obsID = arrayKey.obsID;
+ scanKey.arrayID = arrayKey.arrayID;
+ while (iter != end) {
+ scanKey.scan = *iter;
+ scanKeys.insert(scanKey);
+ ++iter;
+ }
+ return scanKeys;
+}
+
+
+}
+
=======================================
--- /dev/null
+++ /branches/nov14/ms/MeasurementSets/MSKeys.h Tue Mar 3 12:22:11 2015 UTC
@@ -0,0 +1,99 @@
+//# Copyright (C) 1998,1999,2000,2001
+//# Associated Universities, Inc. Washington DC, USA.
+//#
+//# This library is free software; you can redistribute it and/or modify it
+//# under the terms of the GNU Library General Public License as published
by
+//# the Free Software Foundation; either version 2 of the License, or (at
your
+//# option) any later version.
+//#
+//# This library is distributed in the hope that it will be useful, but
WITHOUT
+//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
+//# License for more details.
+//#
+//# You should have received a copy of the GNU Library General Public
License
+//# along with this library; if not, write to the Free Software Foundation,
+//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
+//#
+//# Correspondence concerning AIPS++ should be addressed as follows:
+//# Internet email: aips2-...@nrao.edu.
+//# Postal address: AIPS++ Project Office
+//# National Radio Astronomy Observatory
+//# 520 Edgemont Road
+//# Charlottesville, VA 22903-2475 USA
+//#
+
+#ifndef MS_MSKEYS_H
+#define MS_MSKEYS_H
+
+#include <casacore/casa/aips.h>
+
+#include <set>
+#include <ostream>
+
+namespace casacore {
+
+class String;
+
+// A sub scan is a unique combination of observation ID, array ID, scan
number,
+// and field ID. Negative values are allowed to indicate all values of the
particular
+// ID are desired.
+struct SubScanKey {
+ Int obsID;
+ Int arrayID;
+ Int scan;
+ Int fieldID;
+};
+
+// define operator<() so it can be used as a key in std::map
+Bool operator<(const SubScanKey& lhs, const SubScanKey& rhs);
+
+String toString(const SubScanKey& subScanKey);
+
+std::ostream& operator<<(std::ostream& os, const SubScanKey& scanKey);
+
+// A scan is a unique combination of observation ID, array ID, and scan
number
+// Negative values are allowed to indicate all values of the particular
+// ID are desired.
+struct ScanKey {
+ Int obsID;
+ Int arrayID;
+ Int scan;
+};
+
+// create a ScanKey from a SubScanKey, just omits the SubScanKey's fieldID
+inline ScanKey scanKey(const SubScanKey& subScanKey) {
+ ScanKey key;
+ key.obsID = subScanKey.obsID;
+ key.arrayID = subScanKey.arrayID;
+ key.scan = subScanKey.scan;
+ return key;
+}
+
+String toString(const ScanKey& scanKey);
+
+// define operator<() so it can be used as a key in std::map
+Bool operator<(const ScanKey& lhs, const ScanKey& rhs);
+
+// extract all the unique scan numbers from the specified scans
+std::set<Int> scanNumbers(const std::set<ScanKey>& scanKeys);
+
+std::ostream& operator<<(std::ostream& os, const ScanKey& scanKey);
+
+// An ArrayKey is a unique combination of observation ID and array ID
+// Negative values are allowed to indicate all values of the particular
+// ID are desired.
+struct ArrayKey {
+ Int obsID;
+ Int arrayID;
+};
+
+// define operator<() so it can be used as a key in std::map
+Bool operator<(const ArrayKey& lhs, const ArrayKey& rhs);
+
+// construct scan keys given a set of scan numbers and an ArrayKey
+std::set<ScanKey> scanKeys(const std::set<Int>& scans, const ArrayKey&
arrayKey);
+
+}
+
+#endif
=======================================
--- /branches/nov14/CMakeLists.txt Wed Feb 11 14:00:31 2015 UTC
+++ /branches/nov14/CMakeLists.txt Tue Mar 3 12:22:11 2015 UTC
@@ -343,3 +343,32 @@
message (STATUS "Boost includes? ....... = ${Boost_INCLUDE_DIRS}")
endif (BUILD_PYTHON)

+
+# List of build variables and defaults.
+# BUILD_PYTHON NO
+# ENABLE_SHARED YES
+# ENABLE_RPATH YES
+# CXX11 NO
+# ENABLE_TABLELOCKING YES
+# USE_HDF5 NO
+# USE_FFTW NO
+# USE_THREADS NO
+# USE_OPENMP NO
+# USE_STACKTRACE NO
+# DATA_DIR ""
+#
+# List of possibly used external packages and where
+# CFITSIO fits
+# WCSLIB coordinates
+# DL casa (optional)
+# READLINE casa (optional)
+# HDF5 casa (optional)
+# BISON tables,images
+# FLEX tables,images
+# LAPACK scimath
+# BLAS scimath
+# FFTW scimath (optional)
+# BOOST python (Boost-Python only)
+# PYTHON python
+# NUMPY python
+
=======================================
--- /branches/nov14/casa/Arrays/Slicer.h Wed Dec 10 08:06:42 2014 UTC
+++ /branches/nov14/casa/Arrays/Slicer.h Tue Mar 3 12:22:11 2015 UTC
@@ -394,6 +394,24 @@
// Are all values fixed (i.e., no MimicSource given)?
Bool isFixed() const;

+ // Set the start and end positions. No explicit checking is done that
+ // the input parameters make sense, so you must be certain if you
+ // call these. These are useful if you have a loop with many iterations
+ // and you do not wish the overhead of creating a new Slicer object
+ // for each iteration if the only thing you are doing is adjusting
+ // the start and end positions. Other than for performance reasons,
+ // these methods should not be called and you should prefer the
+ // error checking provided by constructing a new Slicer object.
+ // Note that the length is not updated, so in principle care should
+ // be taken that the length does not change.
+ // <group>
+ void setStart (const IPosition& start)
+ { start_p = start; }
+ void setEnd (const IPosition& end)
+ { end_p = end; }
+ // </group>
+
+
private:
LengthOrLast asEnd_p;
IPosition start_p;
@@ -410,7 +428,7 @@

// Check the given start, end/length and stride.
// Fill in the length or end.
- // It also call <src>fillFixed</src> to fill the fixed flag.
+ // It also calls <src>fillFixed</src> to fill the fixed flag.
void fillEndLen();

// Fill in start, len and stride from a Slice.
=======================================
--- /branches/nov14/casa/BasicMath/Math.h Wed Dec 10 08:06:42 2014 UTC
+++ /branches/nov14/casa/BasicMath/Math.h Tue Mar 3 12:22:11 2015 UTC
@@ -297,7 +297,7 @@
}
inline Bool isNaN(Double val)
{
- return ( isnan(val) );
+ return ( std::isnan(val) );
}
// </group>

=======================================
--- /branches/nov14/casa/Logging/LogOrigin.cc Wed Dec 10 08:06:42 2014 UTC
+++ /branches/nov14/casa/Logging/LogOrigin.cc Tue Mar 3 12:22:11 2015 UTC
@@ -26,13 +26,13 @@
//# $Id$

#include <casacore/casa/Logging/LogOrigin.h>
-
+#include <casacore/casa/OS/EnvVar.h>
#include <casacore/casa/sstream.h>

namespace casacore { //# NAMESPACE CASACORE - BEGIN

LogOrigin::LogOrigin()
- : task_p(""), function_p(""), class_p(""), id_p(True), line_p(0),
file_p("")
+: task_p(""), function_p(""), class_p(""), id_p(True), line_p(0),
file_p(""), node_p(getNode())
{
// Nothing
}
@@ -40,7 +40,8 @@
LogOrigin::LogOrigin(const String &globalFunctionName, const
SourceLocation *where)
: task_p(""), function_p(globalFunctionName), class_p(""), id_p(True),
line_p(where ? where->lineNumber : 0),
- file_p(where ? where->fileName : "")
+ file_p(where ? where->fileName : ""),
+ node_p(getNode())
{
// Nothing
}
@@ -49,7 +50,8 @@
const SourceLocation *where)
: task_p(""), function_p(memberFuncName), class_p(className), id_p(True),
line_p(where ? where->lineNumber : 0),
- file_p(where ? where->fileName : "")
+ file_p(where ? where->fileName : ""),
+ node_p(getNode())
{
// Nothing
}
@@ -59,7 +61,8 @@
: task_p(""), function_p(memberFuncName), class_p(className),
id_p(id),
line_p(where ? where->lineNumber : 0),
- file_p(where ? where->fileName : "")
+ file_p(where ? where->fileName : ""),
+ node_p(getNode())
{
// Nothing
}
@@ -85,6 +88,7 @@
id_p = other.id_p;
line_p = other.line_p;
file_p = other.file_p;
+ node_p = other.node_p;
}

LogOrigin::~LogOrigin()
@@ -179,6 +183,8 @@
if(task_p.length())
nameTag = task_p + "::";
nameTag += className() + "::" + functionName();
+ if(node_p.length())
+ nameTag += "::" + node_p;
return nameTag;
}

@@ -230,6 +236,25 @@
return &permanent;
}

+// Get the OpenMPI rank from the current process.
+// This assumes that the MPI library implementation is openMPI.
+// If a program is started without mpirun, LogOrigin will print nothing
+// for node_p. If it is started with mpirun, each instance
+// running on an MPI server will print its rank at the LogOrigin
+String LogOrigin::getNode()
+{
+ // Note, MPI rank 0 means the MPIClient from mpi4py is running, not a
server,
+ // therefore it should print nothing
+ String rank = EnvironmentVariable::get("OMPI_COMM_WORLD_RANK");
+ if (! rank.empty()) {
+ if (rank == "0") {
+ rank = String();
+ } else {
+ rank = "MPIServer-" + rank;
+ }
+ }
+ return rank;
+}

} //# NAMESPACE CASACORE - END

=======================================
--- /branches/nov14/casa/Logging/LogOrigin.h Wed Dec 10 08:06:42 2014 UTC
+++ /branches/nov14/casa/Logging/LogOrigin.h Tue Mar 3 12:22:11 2015 UTC
@@ -94,6 +94,7 @@
class LogOrigin
{
public:
+
// The default constructor sets a null class name, function name,
object id,
// source file name, and sets the line number to zero.
LogOrigin();
@@ -155,8 +156,8 @@
// <src>where</src> will be defined with the <src>WHERE</src> macro.
LogOrigin &sourceLocation(const SourceLocation *where);

- // Returns <src>class\::function</src> for a member function, or
- // <src>\::function</src> for a global function.
+ // Returns <src>class::function</src> for a member function, or
+ // <src>::function</src> for a global function.
String fullName() const;

// Turn the entire origin into a String.
@@ -168,6 +169,7 @@

// Return true if the line number and file name are not set.
Bool isUnset() const;
+
private:
String task_p;
String function_p;
@@ -175,6 +177,10 @@
ObjectID id_p;
uInt line_p;
String file_p;
+ String node_p;
+
+ // Return a String with the MPI rank
+ String getNode();

// Provide common implementation for copy constructor and
// assignment operator.
=======================================
--- /branches/nov14/casa/Utilities/CountedPtr.h Wed Dec 31 15:18:51 2014 UTC
+++ /branches/nov14/casa/Utilities/CountedPtr.h Tue Mar 3 12:22:11 2015 UTC
@@ -79,6 +79,14 @@
template<class t>
class CountedPtr
{
+#ifdef AIPS_CXX11
+ typedef std::shared_ptr<t> PointerRep;
+ ///#elif HAVE_BOOST
+ /// typedef boost::shared_ptr<t> PointerRep;
+#else
+ typedef std::tr1::shared_ptr<t> PointerRep;
+#endif
+

protected:
// Helper class to make deletion of object optional.
@@ -114,12 +122,17 @@
: pointerRep_p (val, Deleter<t> (delit))
{}

+ // Create from a shared_ptr.
+ CountedPtr (const PointerRep& rep)
+ : pointerRep_p (rep)
+ {}
+
// This copy constructor allows <src>CountedPtr</src>s to be
// initialized from other <src>CountedPtr</src>s for which the pointer
TP*
// is convertible to T*.
template<typename TP>
CountedPtr(const CountedPtr<TP>& that)
- : pointerRep_p(that.pointerRep_p)
+ : pointerRep_p (that.pointerRep_p)
{}

// This destructor only deletes the really stored data when it was
@@ -248,19 +261,6 @@
private:
// Make all types of CountedPtr a friend for the templated operator=.
template<typename TP> friend class CountedPtr;
-
-#ifdef AIPS_CXX11
- typedef std::shared_ptr<t> PointerRep;
- ///#elif HAVE_BOOST
- /// typedef boost::shared_ptr<t> PointerRep;
-#else
- typedef std::tr1::shared_ptr<t> PointerRep;
-#endif
-
- // Create from a shared_ptr.
- CountedPtr (const PointerRep& rep)
- : pointerRep_p (rep)
- {}


PointerRep pointerRep_p;
=======================================
--- /branches/nov14/coordinates/Coordinates/DirectionCoordinate.cc Wed Dec
10 08:06:42 2014 UTC
+++ /branches/nov14/coordinates/Coordinates/DirectionCoordinate.cc Tue Mar
3 12:22:11 2015 UTC
@@ -1309,7 +1309,7 @@
}
for (uInt i=0; i<thisVal.nelements(); i++) {
if (!exclude[i]) {
- if (!casacore::near(thisVal[i],thatVal[i])) {
+ if (!casacore::near(thisVal[i],thatVal[i], tol)) {
oss << "The DirectionCoordinates have differing reference
values for axis "
<< i;
set_error(String(oss));
=======================================
--- /branches/nov14/fits/FITS/blockio.cc Mon Feb 9 07:42:55 2015 UTC
+++ /branches/nov14/fits/FITS/blockio.cc Tue Mar 3 12:22:11 2015 UTC
@@ -176,7 +176,12 @@
}

fits_clear_Fptr( fptr->Fptr, status); // clear Fptr address
- free((fptr->Fptr)->iobuffer); // free memory for I/O
buffers
+ // iobuffer was added with version 3.181... cfitsio 3.03-3.14 do not
have this...
+ // cfitsio 3.03 defines CFITSIO_VERSION to be 3.03 (which is greek to
CPP)...
+ // sometime after 3.181, a translator came along and added CFITSIO_MINOR
as a separate #define
+#ifdef CFITSIO_MINOR
+ free((fptr->Fptr)->iobuffer); // free memory for I/O
buffers
+#endif
free((fptr->Fptr)->headstart); // free memory for headstart
array
free((fptr->Fptr)->filename); // free memory for the filename
(fptr->Fptr)->filename = 0;
=======================================
--- /branches/nov14/lattices/CMakeLists.txt Wed Jan 28 15:31:24 2015 UTC
+++ /branches/nov14/lattices/CMakeLists.txt Tue Mar 3 12:22:11 2015 UTC
@@ -125,6 +125,8 @@
)

install (FILES
+LatticeMath/StatsTiledCollapser.h
+LatticeMath/StatsTiledCollapser.tcc
LatticeMath/CLIPNearest2D.h
LatticeMath/CLIPNearest2D.tcc
LatticeMath/CLInterpolator2D.h
=======================================
--- /branches/nov14/lattices/LatticeMath/LatticeStatistics.h Mon Feb 2
10:04:25 2015 UTC
+++ /branches/nov14/lattices/LatticeMath/LatticeStatistics.h Tue Mar 3
12:22:11 2015 UTC
@@ -37,6 +37,8 @@
#include <casacore/lattices/LatticeMath/TiledCollapser.h>
#include <casacore/lattices/LatticeMath/TiledCollapser.h>
#include <casacore/lattices/LEL/LatticeExprNode.h>
+#include <casacore/lattices/LatticeMath/LatticeStatsDataProvider.h>
+#include <casacore/lattices/LatticeMath/MaskedLatticeStatsDataProvider.h>
#include <casacore/scimath/Mathematics/NumericTraits.h>
#include <casacore/casa/System/PGPlotter.h>
#include <casacore/casa/Utilities/DataType.h>
@@ -203,6 +205,7 @@
typedef typename NumericTraits<T>::PrecisionType AccumType;

struct AlgConf {
+ StatisticsData::ALGORITHM algorithm;
// hinges-fences f factor
Double hf;
// fit to have center type
@@ -211,6 +214,10 @@
FitToHalfStatisticsData::USE_DATA ud;
// fit to half center value
AccumType cv;
+ // Chauvenet zscore
+ Double zs;
+ // Chauvenet max iterations
+ Int mi;
};

// Constructor takes the lattice and a <src>LogIO</src> object for logging.
@@ -371,24 +378,47 @@
// Did we construct with a logger ?
Bool hasLogger () const {return haveLogger_p;};

- // set the statistics algorithm to use. This will reset any previous
algorithm
- // configuration data, so if you want the algorithm to have a different
configuration
- // beyond its default, you will need to call the appropriate
configuration methods
- // after you have called this method.
- void setAlgorithm(
- StatisticsData::ALGORITHM algorithm
+ // configure object to use Classical Statistics
+ // The time, t_x, it takes to compute classical statistics using
algorithm x, can
+ // be modeled by
+ // t_x = n_sets*(a_x + b_x*n_el)
+ // where n_sets is the number of independent sets of data to compute
stats on,
+ // each containing n_el number of elements. a_x is the time it takes to
compute
+ // stats a a single set of data, and b_x is the time it takes to
accumulate
+ // a single point.
+ // The old algorithm was developed in the early history of the project,
I'm guessing
+ // by Neil Kileen, while the new algorithm was developed in 2015 by
Dave Mehringer
+ // as part of the stats framework project. The old algorithm is faster
in the regime
+ // of large n_sets and small n_el, while the new algorithm is faster in
the
+ // regime of small n_sets and large n_el.
+ // If one always wants to use one of these algorithms, that algorithm's
coefficients
+ // should be set to 0, while setting the other algorithm's coefficients
to positive
+ // values. Note that it's the relative, not the absolute, values of
these
+ // coeffecients that is important
+ // The version that takes no parameters uses the default values of the
coefficients;
+ // <group>
+ void configureClassical();
+
+ void configureClassical(Double aOld, Double bOld, Double aNew, Double
bNew);
+ // </group>
+
+ // configure to use fit to half algorithm.
+ void configureFitToHalf(
+ FitToHalfStatisticsData::CENTER
centerType=FitToHalfStatisticsData::CMEAN,
+ FitToHalfStatisticsData::USE_DATA
useData=FitToHalfStatisticsData::LE_CENTER,
+ AccumType centerValue=0
);

- // throws exception if the algorithm has not been set to fit to half.
- // <src>centerValue</src> is only used if <src>centerType</src>=CVALUE
- void configureFitToHalf(
- FitToHalfStatisticsData::CENTER
centerType=FitToHalfStatisticsData::CMEAN,
- FitToHalfStatisticsData::USE_DATA
useData=FitToHalfStatisticsData::LE_CENTER,
- AccumType centerValue=0
- );
+ // configure to use hinges-fences algorithm
+ void configureHingesFences(Double f);
+
+ // configure to use Chauvenet's criterion
+ void configureChauvenet(
+ Double zscore=-1, Int maxIterations=-1
+ );

- // throws exception if the algorithm has not been set to hinges-fences
- void configureHingesFences(Double f);
+ // get number of iterations associated with Chauvenet criterion
algorithm
+ std::map<String, uInt> getChauvenetNiter() const { return _chauvIters; }

protected:

@@ -472,7 +502,7 @@
private:

const MaskedLattice<T>* pInLattice_p;
- TempLattice<AccumType>* pStoreLattice_p;
+ CountedPtr<TempLattice<AccumType> > pStoreLattice_p;
Vector<Int> nxy_p, statsToPlot_p;
Vector<T> range_p;
PGPlotter plotter_p;
@@ -484,10 +514,20 @@
T minFull_p, maxFull_p;
Bool doneFullMinMax_p;

- vector<CountedPtr<StatisticsAlgorithm<AccumType, const T*, const Bool*>
> > _sa;
-
- StatisticsData::ALGORITHM _algorithm;
AlgConf _algConf;
+ std::map<String, uInt> _chauvIters;
+
+ Double _aOld, _bOld, _aNew, _bNew;
+
+ void _setDefaultCoeffs() {
+ // coefficients from timings run on PagedImages on
+ // etacarinae.cv.nrao.edu (dmehring's development
+ // machine)
+ _aOld = 4.7e-7;
+ _bOld = 2.3e-8;
+ _aNew = 1.6e-5;
+ _bNew = 1.5e-8;
+ }

// Summarize the statistics found over the entire lattice
virtual void summStats();
@@ -599,6 +639,15 @@
void _latticePosToStoragePos(
IPosition& storagePos, const IPosition& latticePos
);
+
+ CountedPtr<StatisticsAlgorithm<AccumType, const T*, const Bool*> >
_createStatsAlgorithm() const;
+
+ void _configureDataProviders(
+ LatticeStatsDataProvider<T>& lattDP,
+ MaskedLatticeStatsDataProvider<T>& maskedLattDP
+ ) const;
+
+ void _doStatsLoop(uInt nsets, CountedPtr<LattStatsProgress>
progressMeter);
};

// <summary> Generate statistics, tile by tile, from a masked lattice
</summary>
@@ -669,6 +718,7 @@
// <li>
// </todo>

+/*
template <class T, class U=T>
class StatsTiledCollapser : public TiledCollapser<T,U>
{
@@ -719,6 +769,7 @@
DataType _type;
vector<Bool> _first;
};
+*/

} //# NAMESPACE CASACORE - END

=======================================
--- /branches/nov14/lattices/LatticeMath/LatticeStatistics.tcc Mon Feb 16
07:05:15 2015 UTC
+++ /branches/nov14/lattices/LatticeMath/LatticeStatistics.tcc Tue Mar 3
12:22:11 2015 UTC
@@ -31,6 +31,7 @@
#include <casacore/lattices/LatticeMath/LatticeStatistics.h>
#include <casacore/lattices/LatticeMath/LattStatsSpecialize.h>
#include <casacore/lattices/LatticeMath/LattStatsProgress.h>
+#include <casacore/lattices/LatticeMath/StatsTiledCollapser.h>

#include <casacore/casa/Arrays/ArrayLogical.h>
#include <casacore/casa/Arrays/ArrayMath.h>
@@ -70,6 +71,7 @@

#include <casacore/casa/OS/Timer.h>

+#include <casacore/scimath/Mathematics/ChauvenetCriterionStatistics.h>
#include <casacore/scimath/Mathematics/ClassicalStatistics.h>
#include <casacore/scimath/Mathematics/FitToHalfStatistics.h>
#include <casacore/scimath/Mathematics/HingesFencesStatistics.h>
@@ -101,14 +103,14 @@
showProgress_p(showProgress),
forceDisk_p(forceDisk),
doneFullMinMax_p(False),
- _sa(), _algorithm(StatisticsData::CLASSICAL), _algConf() {
+ _algConf(), _chauvIters() {
nxy_p.resize(0);
statsToPlot_p.resize(0);
range_p.resize(0);
minPos_p.resize(0);
maxPos_p.resize(0);
blcParent_p.resize(0);
-
+ configureClassical();
if (setNewLattice(lattice)) {

// Cursor axes defaults to all
@@ -144,7 +146,7 @@
showProgress_p(showProgress),
forceDisk_p(forceDisk),
doneFullMinMax_p(False),
- _sa(), _algorithm(StatisticsData::CLASSICAL), _algConf()
+ _algConf(), _chauvIters()
{
nxy_p.resize(0);
statsToPlot_p.resize(0);
@@ -152,7 +154,7 @@
minPos_p.resize(0);
maxPos_p.resize(0);
blcParent_p.resize(0);
-//
+ configureClassical();
if (setNewLattice(lattice)) {

// Cursor axes defaults to all
@@ -168,7 +170,8 @@
template <class T>
LatticeStatistics<T>::LatticeStatistics(const LatticeStatistics<T> &other)
: pInLattice_p(0), pStoreLattice_p(0),
- _sa(), _algorithm(StatisticsData::CLASSICAL), _algConf()
+ _algConf(other._algConf), _chauvIters(other._chauvIters),
+ _aOld(other._aOld), _bOld(other._bOld), _aNew(other._aNew),
_bNew(other._bNew)
//
// Copy constructor. Storage lattice is not copied.
//
@@ -193,10 +196,11 @@
pInLattice_p = other.pInLattice_p->cloneML();
// Delete storage lattice

- if (pStoreLattice_p != 0) {
- delete pStoreLattice_p;
+ if (! pStoreLattice_p.null()) {
+ // delete pStoreLattice_p;
pStoreLattice_p = 0;
}
+
needStorageLattice_p = True;
// Do the rest

@@ -234,9 +238,12 @@
doneFullMinMax_p= other.doneFullMinMax_p;
minFull_p = other.minFull_p;
maxFull_p = other.maxFull_p;
- _sa.resize(0);
- _algorithm = other._algorithm;
_algConf = other._algConf;
+ _chauvIters = other._chauvIters;
+ _aNew = other._aNew;
+ _bNew = other._bNew;
+ _aOld = other._aOld;
+ _bOld = other._bOld;
}
return *this;
}
@@ -245,10 +252,9 @@
LatticeStatistics<T>::~LatticeStatistics() {
delete pInLattice_p;
pInLattice_p = 0;
- delete pStoreLattice_p;
- pStoreLattice_p = 0;
+ //delete pStoreLattice_p;
+ // pStoreLattice_p = 0;
}
-

template <class T>
Bool LatticeStatistics<T>::setAxes (const Vector<Int>& axes)
@@ -523,7 +529,9 @@
if (!goodParameterStatus_p) {
return False;
}
- if (needStorageLattice_p) generateStorageLattice();
+ if (needStorageLattice_p) {
+ generateStorageLattice();
+ }
if (type==LatticeStatsBase::NPTS) {
return retrieveStorageStatistic(stats, NPTS, dropDeg);
} else if (type==LatticeStatsBase::SUM) {
@@ -783,26 +791,30 @@
}

template <class T>
-void LatticeStatistics<T>::setAlgorithm(
- StatisticsData::ALGORITHM algorithm
+void LatticeStatistics<T>::configureClassical() {
+ _algConf.algorithm = StatisticsData::CLASSICAL;
+ _setDefaultCoeffs();
+}
+
+template <class T>
+void LatticeStatistics<T>::configureClassical(
+ Double aOld, Double bOld, Double aNew, Double bNew
) {
- if (_algorithm != algorithm) {
- _algorithm = algorithm;
- _algConf.hf = -1;
- _algConf.ct = FitToHalfStatisticsData::CMEAN;
- _algConf.ud = FitToHalfStatisticsData::LE_CENTER;
- _algConf.cv = 0;
- needStorageLattice_p = True;
- }
+ _algConf.algorithm = StatisticsData::CLASSICAL;
+ _aOld = aOld;
+ _bOld = bOld;
+ _aNew = aNew;
+ _bNew = bNew;
}
+

template <class T>
void LatticeStatistics<T>::configureHingesFences(Double f) {
- ThrowIf(
- _algorithm != StatisticsData::HINGESFENCES,
- "Logic Error: _algorithm is not set to HingesFences"
- );
- if (! near(f, _algConf.hf)) {
+ if (
+ _algConf.algorithm != StatisticsData::HINGESFENCES
+ || ! near(f, _algConf.hf)
+ ) {
+ _algConf.algorithm = StatisticsData::HINGESFENCES;
_algConf.hf = f;
needStorageLattice_p = True;
}
@@ -814,18 +826,16 @@
FitToHalfStatisticsData::USE_DATA useData,
AccumType centerValue
) {
- ThrowIf(
- _algorithm != StatisticsData::FITTOHALF,
- "Logic Error: _algorithm is not set to FITTOHALF"
- );
if (
- centerType != _algConf.ct
+ _algConf.algorithm != StatisticsData::FITTOHALF
+ || centerType != _algConf.ct
|| useData != _algConf.ud
|| (
centerType == FitToHalfStatisticsData::CVALUE
&& ! near(centerValue, _algConf.cv)
)
) {
+ _algConf.algorithm = StatisticsData::FITTOHALF;
_algConf.ct = centerType;
_algConf.ud = useData;
_algConf.cv = centerValue;
@@ -834,33 +844,53 @@
}

template <class T>
-Bool LatticeStatistics<T>::generateStorageLattice()
+void LatticeStatistics<T>::configureChauvenet(
+ Double zscore, Int maxIterations
+) {
+ if (
+ _algConf.algorithm != StatisticsData::CHAUVENETCRITERION
+ || ! near(zscore, _algConf.zs)
+ || maxIterations != _algConf.mi
+ ) {
+ _algConf.algorithm = StatisticsData::CHAUVENETCRITERION;
+ _algConf.zs = zscore;
+ _algConf.mi = maxIterations;
+ needStorageLattice_p = True;
+ }
+}

-// Iterate through the lattice and generate the storage lattice
-// The shape of the storage lattice is n1, n2, ..., NACCUM
-// where n1, n2 etc are the display axes
-{
-// Set the display axes vector (possibly already set in ::setAxes)
- displayAxes_p.resize(0);
- displayAxes_p = IPosition::otherAxes(pInLattice_p->ndim(),
- cursorAxes_p).asVector();
+template <class T>
+Bool LatticeStatistics<T>::generateStorageLattice() {

-// Work out dimensions of storage lattice (statistics accumulations
-// are along the last axis)
+ // Iterate through the lattice and generate the storage lattice
+ // The shape of the storage lattice is n1, n2, ..., NACCUM
+ // where n1, n2 etc are the display axes
+
+ // Set the display axes vector (possibly already set in ::setAxes)
+ displayAxes_p.resize(0);
+ displayAxes_p = IPosition::otherAxes(
+ pInLattice_p->ndim(), cursorAxes_p
+ ).asVector();
+
+ // Work out dimensions of storage lattice (statistics accumulations
+ // are along the last axis)
IPosition storeLatticeShape;
- LatticeStatsBase::setStorageImageShape(storeLatticeShape, True,
Int(LatticeStatsBase::NACCUM),
- displayAxes_p,
pInLattice_p->shape());
-// Set the storage lattice tile shape to the tile shape of the
-// axes of the parent lattice from which it is created.
-// For the statistics axis, set the tile shape to NACCUM (small).
+ LatticeStatsBase::setStorageImageShape(
+ storeLatticeShape, True, Int(LatticeStatsBase::NACCUM),
+ displayAxes_p, pInLattice_p->shape()
+ );
+
+ // Set the storage lattice tile shape to the tile shape of the
+ // axes of the parent lattice from which it is created.
+ // For the statistics axis, set the tile shape to NACCUM (small).

IPosition tileShape(storeLatticeShape.nelements(),1);
for (uInt i=0; i<tileShape.nelements()-1; i++) {
tileShape(i) = pInLattice_p->niceCursorShape()(displayAxes_p(i));
}
tileShape(tileShape.nelements()-1) =
storeLatticeShape(storeLatticeShape.nelements()-1);
-// Create storage lattice. If lattice is > 10% of available memory,
-// put it on disk.
+ // Create storage lattice. If lattice is > 10% of available memory,
+ // put it on disk.
uInt memory = HostInfo::memoryTotal()/1024;
Double useMemory = Double(memory)/10.0;
if (forceDisk_p) useMemory = 0.0;
@@ -868,183 +898,222 @@
os_p << LogIO::NORMAL1
<< "Creating new statistics storage lattice of shape " <<
storeLatticeShape << endl << LogIO::POST;
}
- delete pStoreLattice_p;
- pStoreLattice_p = new
TempLattice<AccumType>(TiledShape(storeLatticeShape,
- tileShape), useMemory);
-// Set up min/max location variables
+ pStoreLattice_p = new TempLattice<AccumType>(
+ TiledShape(storeLatticeShape,
+ tileShape), useMemory
+ );
+ // Set up min/max location variables

- //minPos_p.resize(pInLattice_p->shape().nelements());
- //maxPos_p.resize(pInLattice_p->shape().nelements());
- maxPos_p.resize(0);
- minPos_p.resize(0);
+ CountedPtr<LattStatsProgress> pProgressMeter = showProgress_p
+ ? new LattStatsProgress() : NULL;
+ Double timeOld = 0;
+ Double timeNew = 0;
+ uInt nsets = pStoreLattice_p->size()/storeLatticeShape.getLast(1)[0];

-
- CountedPtr<LattStatsProgress> pProgressMeter;
- if (showProgress_p) {
- pProgressMeter = new LattStatsProgress();
+ if (_algConf.algorithm == StatisticsData::CLASSICAL) {
+ uInt nel = pInLattice_p->size()/nsets;
+ timeOld = nsets*(_aOld + _bOld*nel);
+ timeNew = nsets*(_aNew + _bNew*nel);
}
- const uInt nCursorAxes = cursorAxes_p.nelements();
- const IPosition latticeShape(pInLattice_p->shape());
- IPosition cursorShape(pInLattice_p->ndim(),1);
- for (uInt i=0; i<nCursorAxes; i++) {
- cursorShape(cursorAxes_p(i)) = latticeShape(cursorAxes_p(i));
+ //Timer timer;
+ if (
+ _algConf.algorithm == StatisticsData::CLASSICAL
+ && timeOld < timeNew
+ ) {
+ //cout << "using old method" << endl;
+ // use older method for higher performance in the large loop count
+ // regime
+ //timer.mark();
+ minPos_p.resize(pInLattice_p->shape().nelements());
+ maxPos_p.resize(pInLattice_p->shape().nelements());
+ StatsTiledCollapser<T,AccumType> collapser(
+ range_p, noInclude_p, noExclude_p,
+ fixedMinMax_p
+ );
+ Int newOutAxis = pStoreLattice_p->ndim()-1;
+ SubLattice<AccumType> outLatt (*pStoreLattice_p, True);
+ LatticeApply<T,AccumType>::tiledApply(
+ outLatt, *pInLattice_p,
+ collapser, IPosition(cursorAxes_p),
+ newOutAxis,
+ pProgressMeter.null() ? NULL : &*pProgressMeter
+ );
+ collapser.minMaxPos(minPos_p, maxPos_p);
}
-
- IPosition axisPath = cursorAxes_p;
- axisPath.append(displayAxes_p);
- LatticeStepper stepper(latticeShape, cursorShape, axisPath);
- uInt stepperSteps = pProgressMeter.null()
- ? 0 : latticeShape.product()/cursorShape.product();
- DataRanges range;
- if (! noInclude_p || ! noExclude_p) {
- range.push_back(std::pair<T, T>(range_p[0], range_p[1]));
+ else {
+ //cout << "using new method" << endl;
+ //timer.mark();
+ _doStatsLoop(nsets, pProgressMeter);
}
+ pProgressMeter = NULL;
+ //cout << "time " << timer.real() << endl;

- T currentMax = 0;
- T currentMin = 0;
- T overallMax = 0;
- T overallMin = 0;
- Bool isReal = whatType(&currentMax);
- _sa.resize(storeLatticeShape.product()/storeLatticeShape.last());
- typename vector<CountedPtr<StatisticsAlgorithm<AccumType, const T*,
const Bool*> > >::iterator curSA = _sa.begin();
- for (stepper.reset(); ! stepper.atEnd(); stepper++, ++curSA) {
- IPosition curPos = stepper.position();
- IPosition posMax = locInStorageLattice(curPos, LatticeStatsBase::MAX);
- IPosition posMean = locInStorageLattice(curPos,
LatticeStatsBase::MEAN);
- IPosition posMin = locInStorageLattice(curPos, LatticeStatsBase::MIN);
- IPosition posNpts= locInStorageLattice(curPos,
LatticeStatsBase::NPTS);
- IPosition posSum = locInStorageLattice(curPos, LatticeStatsBase::SUM);
- IPosition posSumsq = locInStorageLattice(curPos,
LatticeStatsBase::SUMSQ);
- IPosition posVariance = locInStorageLattice(curPos,
LatticeStatsBase::VARIANCE);
-
- // Create SubLattice from chunk
- Slicer slicer(curPos, stepper.endPosition(), Slicer::endIsLast);
- SubLattice<T> subLat(*pInLattice_p, slicer);
- if (! pProgressMeter.null() && stepper.atStart()) {
- uInt nSublatticeSteps =
subLat.shape().product()/subLat.niceCursorShape().product();
- pProgressMeter->init(stepperSteps*nSublatticeSteps);
- }
-
- // we use T rather than AccumType for the first templated variable
because
- // we lose no accuracy and if AccumType=complex<double> and
T=complex<float>,
- // this code will not compile
-
- LatticeStatsDataProviderBase<T> *dataProvider = NULL;
- if (subLat.isMasked()) {
- dataProvider = new MaskedLatticeStatsDataProvider<T>(subLat);
- }
- else {
- dataProvider = new LatticeStatsDataProvider<T>(subLat);
- }
-
- /*
- if (! pProgressMeter.null()) {
- dataProvider->setProgressMeter(pProgressMeter);
- }
- */
- // FIXME having Bool variables that are true when a fundamental
property is negated
- // is incredibly confusing and certainly not best practice. Rename
and adjust
- // the meaning of noInclude_p and noExclude_p
- if (! noInclude_p || ! noExclude_p) {
- dataProvider->setRanges(range, ! noInclude_p);
- }
-
- // its annoying that valid implicit casting of CountedPtr's won't
compile, so
- // we have to do this the hard way. The ThrowIf statement can be
removed once
- // this has been thoroughly exercised.
- CountedPtr<StatsDataProvider<AccumType, const T*, const Bool*> > mydp
- = dynamic_cast<StatsDataProvider<AccumType, const T*, const Bool*>
*>(dataProvider);
- ThrowIf (mydp.null(), "Logic Error: dynamic cast failed");
- switch (_algorithm) {
- case StatisticsData::CLASSICAL:
- *curSA = new ClassicalStatistics<AccumType, const T*, const Bool*>();
- break;
- case StatisticsData::HINGESFENCES: {
- *curSA = new HingesFencesStatistics<AccumType, const T*, const
Bool*>(_algConf.hf);
- break;
- }
- case StatisticsData::FITTOHALF: {
- *curSA = new FitToHalfStatistics<AccumType, const T*, const Bool*>(
- _algConf.ct, _algConf.ud, _algConf.cv
- );
- break;
- }
- default:
- ThrowCc("Logic Error: Unhandled algorithm " +
String::toString(_algorithm));
- }
-
- (*curSA)->setDataProvider(mydp);
- if (! pProgressMeter.null()) {
- dataProvider->setProgressMeter(CountedPtr<LattStatsProgress>(NULL));
- }
- StatsData<AccumType> stats = (*curSA)->getStatistics();
- pStoreLattice_p->putAt(stats.mean, posMean);
- pStoreLattice_p->putAt(stats.npts, posNpts);
- pStoreLattice_p->putAt(stats.sum, posSum);
- pStoreLattice_p->putAt(stats.sumsq, posSumsq);
- pStoreLattice_p->putAt(stats.variance, posVariance);
- if (fixedMinMax_p && ! noInclude_p) {
- currentMax = range[0].second;
- }
- else if (! stats.max.null()) {
- currentMax = *stats.max;
- }
- pStoreLattice_p->putAt(currentMax, posMax);
- if (fixedMinMax_p && ! noInclude_p) {
- currentMin = range[0].first;
- }
- else if (! stats.min.null()) {
- currentMin = *stats.min;
- }
- pStoreLattice_p->putAt(currentMin, posMin);
- if (isReal) {
- // CAUTION The way this has worked in the past apparently for
- // lattices is that the max and min positions are representative
- // of the *entire* lattice, and were not stored on a sublattice
- // by sublattice basis. This is easy to fix now,
- // but for backward compatibility, I'm leaving this functionality as
- // it has been.
- if (! fixedMinMax_p || noInclude_p) {
- if (stepper.atStart()) {
- IPosition myMaxPos, myMinPos;
- dataProvider->minMaxPos(myMinPos, myMaxPos);
- if (myMinPos.size() > 0) {
- minPos_p = subLat.positionInParent(myMinPos);
- }
- if (myMaxPos.size() > 0) {
- maxPos_p = subLat.positionInParent(myMaxPos);
- }
- overallMin = currentMin;
- overallMax = currentMax;
- }
- else if (
- currentMax > overallMax || currentMin < overallMin
- ) {
- IPosition myMaxPos, myMinPos;
- dataProvider->minMaxPos(myMinPos, myMaxPos);
- if (currentMin < overallMin) {
- if (myMinPos.size() > 0) {
- minPos_p = subLat.positionInParent(myMinPos);
- }
- overallMin = currentMin;
- }
- if (currentMax > overallMax) {
- if (myMaxPos.size() > 0) {
- maxPos_p = subLat.positionInParent(myMaxPos);
- }
- overallMax = currentMax;
- }
- }
- }
- }
- }
- pProgressMeter = NULL;
// Do robust statistics separately as required.
generateRobust();
- needStorageLattice_p = False;
+ needStorageLattice_p = False;
doneSomeGoodPoints_p = False;
return True;
}
+
+template <class T>
+void LatticeStatistics<T>::_doStatsLoop(
+ uInt nsets, CountedPtr<LattStatsProgress> progressMeter
+) {
+ // use high performance method for low iteration count
+ maxPos_p.resize(0);
+ minPos_p.resize(0);
+ const uInt nCursorAxes = cursorAxes_p.nelements();
+ const IPosition latticeShape(pInLattice_p->shape());
+ IPosition cursorShape(pInLattice_p->ndim(),1);
+ for (uInt i=0; i<nCursorAxes; i++) {
+ cursorShape(cursorAxes_p(i)) = latticeShape(cursorAxes_p(i));
+ }
+ IPosition axisPath = cursorAxes_p;
+ axisPath.append(displayAxes_p);
+ LatticeStepper stepper(latticeShape, cursorShape, axisPath);
+ T currentMax = 0;
+ T currentMin = 0;
+ T overallMax = 0;
+ T overallMin = 0;
+ Bool isReal = whatType(&currentMax);
+
+ CountedPtr<StatisticsAlgorithm<AccumType, const T*, const Bool*> > sa =
_createStatsAlgorithm();
+ LatticeStatsDataProvider<T> lattDP;
+ MaskedLatticeStatsDataProvider<T> maskedLattDP;
+ LatticeStatsDataProviderBase<T> *dataProvider;
+ _configureDataProviders(lattDP, maskedLattDP);
+ Bool nsetsIsLarge = nsets > 50;
+ if (! progressMeter.null()) {
+ if (nsetsIsLarge) {
+ progressMeter->init(nsets);
+ }
+ else {
+ lattDP.setProgressMeter(progressMeter);
+ if (pInLattice_p->isMasked()) {
+ maskedLattDP.setProgressMeter(progressMeter);
+ }
+ }
+ }
+ _chauvIters.clear();
+ IPosition curPos, posMax, posMean, posMin, posNpts,
+ posSum, posSumsq, posVariance;
+ Slicer slicer(stepper.position(), stepper.endPosition(),
Slicer::endIsLast);
+ SubLattice<T> subLat(*pInLattice_p, slicer);
+ StatsData<AccumType> stats;
+ for (stepper.reset(); ! stepper.atEnd(); stepper++) {
+ curPos = stepper.position();
+ posMax = locInStorageLattice(curPos, LatticeStatsBase::MAX);
+ posMean = locInStorageLattice(curPos, LatticeStatsBase::MEAN);
+ posMin = locInStorageLattice(curPos, LatticeStatsBase::MIN);
+ posNpts = locInStorageLattice(curPos, LatticeStatsBase::NPTS);
+ posSum = locInStorageLattice(curPos, LatticeStatsBase::SUM);
+ posSumsq = locInStorageLattice(curPos, LatticeStatsBase::SUMSQ);
+ posVariance = locInStorageLattice(curPos, LatticeStatsBase::VARIANCE);
+ slicer.setStart(curPos);
+ slicer.setEnd(stepper.endPosition());
+ subLat.setRegion(slicer);
+ if (
+ stepper.atStart() && ! progressMeter.null()
+ && ! nsetsIsLarge
+ ) {
+ uInt nelem = subLat.size();
+ uInt nSublatticeSteps = nelem > 4096*4096
+ ? nelem/subLat.advisedMaxPixels()
+ : 1;
+ /*
+ cout << "nelem " << nelem << endl;
+ cout << "nice cursor shape " << subLat.niceCursorShape() << endl;
+ cout << "nsteps " << nsets*nSublatticeSteps << endl;
+ */
+ progressMeter->init(nsets*nSublatticeSteps);
+ }
+ if(subLat.isMasked()) {
+ maskedLattDP.setLattice(subLat);
+ dataProvider = &maskedLattDP;
+ }
+ else {
+ lattDP.setLattice(subLat);
+ dataProvider = &lattDP;
+ }
+ sa->setDataProvider(dataProvider);
+ stats = sa->getStatistics();
+ if (_algConf.algorithm == StatisticsData::CHAUVENETCRITERION) {
+ ChauvenetCriterionStatistics<AccumType, const T*, const Bool*> *ch
+ = dynamic_cast<ChauvenetCriterionStatistics<AccumType, const T*, const
Bool*> *>(
+ &*sa
+ );
+ ostringstream os;
+ os << curPos;
+ // using strings as keys rather than the IPosition objects directly
because for some reason,
+ // only one IPosition gets added to the map, and then no other ones get
added.
+ // I don't understand, things seem to work OK when I try this in
tIPosition, but not here.
+ _chauvIters[os.str()] = ch->getNiter();
+ }
+ pStoreLattice_p->putAt(stats.mean, posMean);
+ pStoreLattice_p->putAt(stats.npts, posNpts);
+ pStoreLattice_p->putAt(stats.sum, posSum);
+ pStoreLattice_p->putAt(stats.sumsq, posSumsq);
+ pStoreLattice_p->putAt(stats.variance, posVariance);
+ if (fixedMinMax_p && ! noInclude_p) {
+ currentMax = range_p[1];
+ }
+ else if (! stats.max.null()) {
+ currentMax = *stats.max;
+ }
+ pStoreLattice_p->putAt(currentMax, posMax);
+ if (fixedMinMax_p && ! noInclude_p) {
+ currentMin = range_p[0];
+ }
+ else if (! stats.min.null()) {
+ currentMin = *stats.min;
+ }
+ pStoreLattice_p->putAt(currentMin, posMin);
+ if (isReal) {
+ // CAUTION The way this has worked in the past apparently for
+ // lattices is that the max and min positions are representative
+ // of the *entire* lattice, and were not stored on a sublattice
+ // by sublattice basis. This is easy to fix now,
+ // but for backward compatibility, I'm leaving this functionality as
+ // it has been.
+ if (! fixedMinMax_p || noInclude_p) {
+ if (stepper.atStart()) {
+ IPosition myMaxPos, myMinPos;
+ dataProvider->minMaxPos(myMinPos, myMaxPos);
+ if (myMinPos.size() > 0) {
+ minPos_p = subLat.positionInParent(myMinPos);
+ }
+ if (myMaxPos.size() > 0) {
+ maxPos_p = subLat.positionInParent(myMaxPos);
+ }
+ overallMin = currentMin;
+ overallMax = currentMax;
+ }
+ else if (
+ currentMax > overallMax || currentMin < overallMin
+ ) {
+ IPosition myMaxPos, myMinPos;
+ dataProvider->minMaxPos(myMinPos, myMaxPos);
+ if (currentMin < overallMin) {
+ if (myMinPos.size() > 0) {
+ minPos_p = subLat.positionInParent(myMinPos);
+ }
+ overallMin = currentMin;
+ }
+ if (currentMax > overallMax) {
+ if (myMaxPos.size() > 0) {
+ maxPos_p = subLat.positionInParent(myMaxPos);
+ }
+ overallMax = currentMax;
+ }
+ }
+ }
+ }
+ if(! progressMeter.null() && nsetsIsLarge) {
+ (*progressMeter)++;
+ }
+ }
+}
+

template <class T>
void LatticeStatistics<T>::generateRobust () {
@@ -1060,20 +1129,38 @@
IPosition axisPath = cursorAxes_p;
axisPath.append(displayAxes_p);
LatticeStepper stepper(latticeShape, cursorShape, axisPath);
- std::set<Double> quartiles;
+ std::set<Double> fractions;
+ CountedPtr<StatisticsAlgorithm<AccumType, const T*, const Bool*> > sa;
+ LatticeStatsDataProvider<T> lattDP;
+ MaskedLatticeStatsDataProvider<T> maskedLattDP;
+
+ IPosition curPos, pos, pos2, pos3, posQ1, posQ3,
+ posNpts, posMax, posMin;
+ Slicer slicer;
+ SubLattice<T> subLat;
+ std::map<Double, AccumType> quantiles;
+ CountedPtr<uInt64> knownNpts;
+ CountedPtr<AccumType> knownMax, knownMin;
if (doRobust_p) {
- quartiles.insert(0.25);
- quartiles.insert(0.75);
+ fractions.insert(0.25);
+ fractions.insert(0.75);
+ sa = _createStatsAlgorithm();
+ _configureDataProviders(lattDP, maskedLattDP);
+ slicer = Slicer(stepper.position(), stepper.endPosition(),
Slicer::endIsLast);
+ subLat = SubLattice<T>(*pInLattice_p, slicer);
}
- std::map<Double, AccumType> quantileToValue;
- typename vector<CountedPtr<StatisticsAlgorithm<AccumType, const T*,
const Bool*> > >::iterator curSA = _sa.begin();
for (stepper.reset(); ! stepper.atEnd(); stepper++) {
- IPosition pos = locInStorageLattice(stepper.position(),
LatticeStatsBase::MEDIAN);
- IPosition pos2 = locInStorageLattice(stepper.position(),
LatticeStatsBase::MEDABSDEVMED);
- IPosition pos3 = locInStorageLattice(stepper.position(),
LatticeStatsBase::QUARTILE);
- IPosition posQ1 = locInStorageLattice(stepper.position(),
LatticeStatsBase::Q1);
- IPosition posQ3 = locInStorageLattice(stepper.position(),
LatticeStatsBase::Q3);
- if (! doRobust_p) {
+ curPos = stepper.position();
+ pos = locInStorageLattice(stepper.position(), LatticeStatsBase::MEDIAN);
+ pos2 = locInStorageLattice(stepper.position(),
LatticeStatsBase::MEDABSDEVMED);
+ pos3 = locInStorageLattice(stepper.position(),
LatticeStatsBase::QUARTILE);
+ posQ1 = locInStorageLattice(stepper.position(), LatticeStatsBase::Q1);
+ posQ3 = locInStorageLattice(stepper.position(), LatticeStatsBase::Q3);
+ if (doRobust_p) {
+ posNpts = locInStorageLattice(stepper.position(),
LatticeStatsBase::NPTS);
+ knownNpts = new uInt64((uInt64)abs(pStoreLattice_p->getAt(posNpts)));
+ }
+ if (! doRobust_p || *knownNpts == 0) {
// Stick zero in storage lattice (it's not initialized)
AccumType val(0);
pStoreLattice_p->putAt(val, pos);
@@ -1083,24 +1170,90 @@
pStoreLattice_p->putAt(val, posQ3);
continue;
}
- quantileToValue.clear();
+
+ posMax = locInStorageLattice(stepper.position(), LatticeStatsBase::MAX);
+ posMin = locInStorageLattice(stepper.position(), LatticeStatsBase::MIN);
+
+ quantiles.clear();
+
+ slicer.setStart(curPos);
+ slicer.setEnd(stepper.endPosition());
+ subLat.setRegion(slicer);
+ if (subLat.isMasked()) {
+ maskedLattDP.setLattice(subLat);
+ sa->setDataProvider(&maskedLattDP);
+ }
+ else {
+ lattDP.setLattice(subLat);
+ sa->setDataProvider(&lattDP);
+ }
// computing the median and the quartiles simultaneously minimizes
// the number of necessary data scans, as opposed to first calling
// getMedian() and getQuartiles() separately
- AccumType median = (*curSA)->getMedianAndQuantiles(quantileToValue,
quartiles);
- AccumType q1 = quantileToValue[0.25];
- AccumType q3 = quantileToValue[0.75];
- AccumType iqr = q3 - q1;
- AccumType medabsdevmed = (*curSA)->getMedianAbsDevMed();
- // Whack results into storage lattice
- pStoreLattice_p->putAt(median, pos);
- pStoreLattice_p->putAt(medabsdevmed, pos2);
- pStoreLattice_p->putAt(iqr, pos3);
- pStoreLattice_p->putAt(q1, posQ1);
- pStoreLattice_p->putAt(q3, posQ3);
- ++curSA;
+ knownMin = new AccumType(pStoreLattice_p->getAt(posMin));
+ knownMax = new AccumType(pStoreLattice_p->getAt(posMax));
+ pStoreLattice_p->putAt(
+ sa->getMedianAndQuantiles(
+ quantiles, fractions, knownNpts, knownMin, knownMax
+ ),
+ pos
+ );
+ pStoreLattice_p->putAt(sa->getMedianAbsDevMed(), pos2);
+ pStoreLattice_p->putAt(quantiles[0.75] - quantiles[0.25], pos3);
+ pStoreLattice_p->putAt(quantiles[0.25], posQ1);
+ pStoreLattice_p->putAt(quantiles[0.75], posQ3);
+ }
+}
+
+template <class T>
+CountedPtr<StatisticsAlgorithm<typename LatticeStatistics<T>::AccumType,
const T*, const Bool*> >
+LatticeStatistics<T>::_createStatsAlgorithm() const {
+ CountedPtr<StatisticsAlgorithm<AccumType, const T*, const Bool*> > sa;
+ switch (_algConf.algorithm) {
+ case StatisticsData::CLASSICAL:
+ sa = new ClassicalStatistics<AccumType, const T*, const Bool*>();
+ return sa;
+ case StatisticsData::HINGESFENCES: {
+ sa = new HingesFencesStatistics<AccumType, const T*, const
Bool*>(_algConf.hf);
+ return sa;
+ }
+ case StatisticsData::FITTOHALF: {
+ sa = new FitToHalfStatistics<AccumType, const T*, const Bool*>(
+ _algConf.ct, _algConf.ud, _algConf.cv
+ );
+ return sa;
+ }
+ case StatisticsData::CHAUVENETCRITERION: {
+ sa = new ChauvenetCriterionStatistics<AccumType, const T*, const Bool*>(
+ _algConf.zs, _algConf.mi
+ );
+ return sa;
+ }
+ default:
+ ThrowCc(
+ "Logic Error: Unhandled algorithm "
+ + String::toString(_algConf.algorithm)
+ );
+ }
+}
+
+template <class T>
+void LatticeStatistics<T>::_configureDataProviders(
+ LatticeStatsDataProvider<T>& lattDP,
+ MaskedLatticeStatsDataProvider<T>& maskedLattDP
+) const {
+ if (! noInclude_p || ! noExclude_p) {
+ DataRanges range;
+ if (! noInclude_p || ! noExclude_p) {
+ range.push_back(std::pair<T, T>(range_p[0], range_p[1]));
+ }
+ lattDP.setRanges(range, ! noInclude_p);
+ if (pInLattice_p->isMasked()) {
+ maskedLattDP.setRanges(range, ! noInclude_p);
+ }
}
}
+

template <class T>
void LatticeStatistics<T>::listMinMax(ostringstream& osMin,
@@ -1614,25 +1767,47 @@
Int layer = 0;
for ( pixelIterator.reset(); !pixelIterator.atEnd(); pixelIterator++ )
{
IPosition dPos = pixelIterator.position();
- if (displayAxes_p.nelements() == 2) {
- if (zAx == 1) {
- if (dPos[1] != zLayer) {
- continue;
- } else {
- layer = hLayer;
- }
- }
- if (hAx == 1) {
- if (dPos[1] != hLayer) {
- continue;
- } else {
- layer = zLayer;
- }
- }
- }
- if (displayAxes_p.nelements() == 1) {
- layer = zLayer;
- }
+ /*
+ if (displayAxes_p.nelements() == 2) {
+ if (zAx == 1) {
+ if (dPos[1] != zLayer)
+ continue;
+ else
+ layer = hLayer;
+ }
+ if (hAx == 1) {
+ if (dPos[1] != hLayer)
+ continue;
+ else
+ layer = zLayer;
+ }
+ }
+ if (displayAxes_p.nelements() == 1)
+ layer = zLayer;
+ */
+ if (displayAxes_p.nelements() == 2) {
+ if (zAx == 1) {
+ if (dPos[1] != zLayer) {
+ continue;
+ }
+ else {
+ layer = hLayer;
+ }
+ }
+ if (hAx == 1) {
+ if (dPos[1] != hLayer) {
+ continue;
+ }
+ else {
+ layer = zLayer;
+ }
+ }
+ }
+ if (displayAxes_p.nelements() == 1) {
+ layer = zLayer;
+ }
+
+
Matrix<AccumType> matrix(pixelIterator.matrixCursor());
for (uInt i=0; i<n1; i++) {
const AccumType& nPts = matrix(i,NPTS);
@@ -1936,12 +2111,12 @@
plotter_p.sch (1.2);
plotter_p.svp(0.1,0.9,0.1,0.9);
}
-
// Generate storage lattice if required

if (needStorageLattice_p) {
if (!generateStorageLattice()) return False;
}
+
// If we don't have any display axes just summarise the lattice statistics
if (displayAxes_p.nelements() == 0) {
summStats ();
@@ -1955,14 +2130,12 @@
// Allocate ordinate arrays for plotting and listing. Try to preserve
// the true Type of the data as long as we can. Eventually, for
// plotting we have to make it real valued
-
Matrix<AccumType> ord(n1,NSTATS);

// Iterate through storage lattice by planes (first and last axis of
storage lattice)
// Specify which axes are the matrix axes so that we can discard other
// degenerate axes with the matrixCursor function. n1 is only
// constrained to be n1 >= 1
-
IPosition cursorShape(pStoreLattice_p->ndim(),1);
cursorShape(0) = pStoreLattice_p->shape()(0);
cursorShape(pStoreLattice_p->ndim()-1) =
pStoreLattice_p->shape()(pStoreLattice_p->ndim()-1);
@@ -1976,25 +2149,11 @@

// Get beam area

- /*
- * dmehring 2012may23: Changing beam area from Double to Array<Double> to
support
- * per plane beams. However, I'm quite confused at what this method is
doing and so
- * am not sure how to integrate the beamArea array into it. Is plotting
of statistics
- * even supported any longer? It doesn't seem to be from the casapy
user interface;
- * in ImageAnalysis::statistics, the plotting device is explicitly
set "/NULL" independent
- * of user inputs. Until I'm sure how to handle the Array of beam areas
correctly and have
- * definitive proof plotting of statistics is actually still used, I'm
commenting out the
- * beam specific code here and in the nested loop below.
-
- Array<Double> beamArea;
- Bool hasBeam = _getBeamArea(beamArea);
- */
Bool hasBeam = False;
//
for (pixelIterator.reset(); !pixelIterator.atEnd(); pixelIterator++) {

// Convert accumulations to mean, sigma, and rms.
-
Matrix<AccumType> matrix(pixelIterator.matrixCursor()); //
Reference semantics
for (uInt i=0; i<n1; i++) {
const AccumType& nPts = matrix(i,NPTS);
@@ -2019,7 +2178,6 @@
ord(j,i) = matrix(j,i);
}
}
-

// Plot statistics

@@ -3088,7 +3246,6 @@
const IPosition shape = statsSliceShape();
Array<AccumType> stats(shape);
pStoreLattice_p->getSlice (stats, IPosition(1,0), shape,
IPosition(1,1));
-
IPosition pos(1);
pos(0) = NPTS;
AccumType nPts = stats(pos);
@@ -3099,7 +3256,6 @@

pos(0) = MEDABSDEVMED;
AccumType medAbsDevMed = stats(pos);
-//
pos(0) = QUARTILE;
AccumType quartile= stats(pos);

@@ -3116,17 +3272,13 @@

pos(0) = VARIANCE;
AccumType var = stats(pos);
-//
- // AccumType mean = LattStatsSpecialize::getMean(sum, nPts);
- // AccumType var = LattStatsSpecialize::getVariance(sum, sumSq, nPts);
AccumType rms = LattStatsSpecialize::getRms(sumSq, nPts);
AccumType sigma = LattStatsSpecialize::getSigma(var);
-//
+
pos(0) = MIN;
AccumType dMin = stats(pos);
pos(0) = MAX;
AccumType dMax = stats(pos);
-
// Do this check so that we only print the stats when we have values.
if (LattStatsSpecialize::hasSomePoints(nPts)) {
displayStats(
@@ -3282,14 +3434,14 @@
os.setf(ios::left, ios::adjustfield);
}

+/*
// StatsTiledCollapser

-
template <class T, class U>
StatsTiledCollapser<T,U>::StatsTiledCollapser(const Vector<T>& pixelRange,
Bool noInclude, Bool
noExclude,
Bool fixedMinMax)
-: /*range_p(pixelRange), */
+:
noInclude_p(noInclude),
noExclude_p(noExclude),
fixedMinMax_p(fixedMinMax),
@@ -3314,6 +3466,8 @@
template <class T, class U>
void StatsTiledCollapser<T,U>::initAccumulator (uInt n1, uInt n3)
{
+ PrecTimer timer;
+ timer.start();
n1_p = n1;
n3_p = n3;
_first = vector<Bool>(n1*n3, True);
@@ -3324,6 +3478,8 @@
iter->setCalculateAsAdded(True);
++iter;
}
+ timer.stop();
+ cout << timer.getReal() << endl;
}

template <class T, class U>
@@ -3335,17 +3491,14 @@
const IPosition& shape
***The diff for this file has been truncated for email.***
=======================================
--- /branches/nov14/lattices/LatticeMath/LatticeStatsDataProvider.h Mon Feb
16 07:05:15 2015 UTC
+++ /branches/nov14/lattices/LatticeMath/LatticeStatsDataProvider.h Tue
Mar 3 12:22:11 2015 UTC
@@ -38,9 +38,27 @@

template <class T> class LatticeStatsDataProvider
: public LatticeStatsDataProviderBase<T> {
+
public:

- LatticeStatsDataProvider(Lattice<T>& lattice);
+ // default constructor, must set lattice after construction but before
+ // using the object
+ LatticeStatsDataProvider();
+
+ // <src>iteratorLimitBytes</src> is related to the size of the lattice.
+ // If the lattice is greater than this size, then a lattice iterator will
+ // be used to step through the lattice. If less, then all the data in the
+ // values in the lattice are retrieved in a single chunk. The advantage of
+ // the iterator is that less memory is used. The disadvantage is there is
+ // a significant performace cost, so if the lattice is small, it is
better to
+ // get all its values in a single chunk and forgo the iterator. This is
particularly
+ // true when looping for a large number of iterations and creating a
+ // LatticeStatsDataProvider each loop (in that case, you probably will
want
+ // to create a single object before the loop and use setLattice() to
update
+ // its lattice).
+ LatticeStatsDataProvider(
+ const Lattice<T>& lattice, uInt iteratorLimitBytes=4096*4096
+ );

~LatticeStatsDataProvider();

@@ -73,6 +91,22 @@
// reset the provider to point to the first data set it manages.
void reset();

+ // set the lattice. Automatically resets the lattice iterator
+ // <src>iteratorLimitBytes</src> is related to the size of the lattice.
+ // If the lattice is greater than this size, then a lattice iterator will
+ // be used to step through the lattice. If less, then all the data in the
+ // values in the lattice are retrieved in a single chunk. The advantage of
+ // the iterator is that less memory is used. The disadvantage is there is
+ // a significant performace cost, so if the lattice is small, it is
better to
+ // get all its values in a single chunk and forgo the iterator. This is
particularly
+ // true when looping for a large number of iterations and creating a
+ // LatticeStatsDataProvider each loop (in that case, you probably will
want
+ // to create a single object before the loop and use setLattice() to
update
+ // its lattice).
+ void setLattice(
+ const Lattice<T>& lattice, uInt iteratorLimitBytes=4096*4096
+ );
+
// <group>
// see base class documentation.
void updateMaxPos(const std::pair<Int64, Int64>& maxpos);
@@ -81,10 +115,10 @@
// </group>

private:
- RO_LatticeIterator<T> _iter;
+ CountedPtr<RO_LatticeIterator<T> > _iter;
Array<T> _currentSlice;
const T* _currentPtr;
- Bool _delData;
+ Bool _delData, _atEnd;

void _freeStorage();

=======================================
--- /branches/nov14/lattices/LatticeMath/LatticeStatsDataProvider.tcc Mon
Feb 16 07:05:15 2015 UTC
+++ /branches/nov14/lattices/LatticeMath/LatticeStatsDataProvider.tcc Tue
Mar 3 12:22:11 2015 UTC
@@ -30,13 +30,21 @@
#include <casacore/lattices/LatticeMath/LatticeStatsDataProvider.h>

namespace casacore {
+
+template <class T>
+LatticeStatsDataProvider<T>::LatticeStatsDataProvider()
+ : LatticeStatsDataProviderBase<T>(),
+ _iter(), _currentSlice(),
+ _currentPtr(0), _delData(False), _atEnd(False) {}

template <class T>
LatticeStatsDataProvider<T>::LatticeStatsDataProvider(
- Lattice<T>& lattice
+ const Lattice<T>& lattice, uInt iteratorLimitBytes
) : LatticeStatsDataProviderBase<T>(),
- _iter(RO_LatticeIterator<T>(lattice)), _currentSlice(),
- _currentPtr(0), _delData(False) {}
+ _iter(), _currentSlice(),
+ _currentPtr(0), _delData(False), _atEnd(False) {
+ setLattice(lattice, iteratorLimitBytes);
+}

template <class T>
LatticeStatsDataProvider<T>::~LatticeStatsDataProvider() {}
@@ -44,14 +52,22 @@
template <class T>
void LatticeStatsDataProvider<T>::operator++() {
_freeStorage();
- ++_iter;
+ if (_iter.null()) {
+ _atEnd = True;
+ }
+ else {
+ ++(*_iter);
+ }
this->_updateProgress();
}

template <class T>
uInt LatticeStatsDataProvider<T>::estimatedSteps() const {
- IPosition lattShape = _iter.latticeShape();
- IPosition cursShape = _iter.cursor().shape();
+ if (_iter.null()) {
+ return 1;
+ }
+ IPosition lattShape = _iter->latticeShape();
+ IPosition cursShape = _iter->cursor().shape();
uInt ndim = lattShape.size();
uInt count = 1;
for (uInt i=0; i<ndim; i++) {
@@ -66,7 +82,10 @@

template <class T>
Bool LatticeStatsDataProvider<T>::atEnd() const {
- return _iter.atEnd();
+ if (_iter.null()) {
+ return _atEnd;
+ }
+ return _iter->atEnd();
}

template <class T>
@@ -77,12 +96,17 @@

template <class T>
uInt64 LatticeStatsDataProvider<T>::getCount() {
- return _iter.cursor().size();
+ if (_iter.null()) {
+ return _currentSlice.size();
+ }
+ return _iter->cursor().size();
}

template <class T>
const T* LatticeStatsDataProvider<T>::getData() {
- _currentSlice.assign(_iter.cursor());
+ if (! _iter.null()) {
+ _currentSlice.assign(_iter->cursor());
+ }
_currentPtr = _currentSlice.getStorage(_delData);
return _currentPtr;
}
@@ -100,25 +124,51 @@
template <class T>
void LatticeStatsDataProvider<T>::reset() {
LatticeStatsDataProviderBase<T>::reset();
- _iter.reset();
+ if (! _iter.null()) {
+ _iter->reset();
+ }
+}
+
+template <class T>
+void LatticeStatsDataProvider<T>::setLattice(
+ const Lattice<T>& lattice, uInt iteratorLimitBytes
+) {
+ finalize();
+ if (lattice.size() > iteratorLimitBytes/sizeof(T)) {
+ TileStepper stepper(
+ lattice.shape(), lattice.niceCursorShape(
+ lattice.advisedMaxPixels()
+ )
+ );
+ _iter = new RO_LatticeIterator<T>(lattice, stepper);
+ }
+ else {
+ _iter = NULL;
+ _currentSlice.assign(lattice.get());
+ _atEnd = False;
+ }
}

template <class T>
void LatticeStatsDataProvider<T>::updateMaxPos(
const std::pair<Int64, Int64>& maxpos
) {
- this->_updateMaxPos(
- _iter.position() + toIPositionInArray(maxpos.second,
_currentSlice.shape())
- );
+ IPosition p = toIPositionInArray(maxpos.second, _currentSlice.shape());
+ if (! _iter.null()) {
+ p += _iter->position();
+ }
+ this->_updateMaxPos(p);
}

template <class T>
void LatticeStatsDataProvider<T>::updateMinPos(
const std::pair<Int64, Int64>& minpos
) {
- this->_updateMinPos(
- _iter.position() + toIPositionInArray(minpos.second,
_currentSlice.shape())
- );
+ IPosition p = toIPositionInArray(minpos.second, _currentSlice.shape());
+ if (! _iter.null()) {
+ p += _iter->position();
+ }
+ this->_updateMinPos(p);
}

template <class T>
@@ -126,8 +176,6 @@
_currentSlice.freeStorage (_currentPtr, _delData);
_delData = False;
}
-
}
-

#endif
=======================================
--- /branches/nov14/lattices/LatticeMath/LatticeStatsDataProviderBase.tcc
Mon Feb 16 07:05:15 2015 UTC
+++ /branches/nov14/lattices/LatticeMath/LatticeStatsDataProviderBase.tcc
Tue Mar 3 12:22:11 2015 UTC
@@ -90,9 +90,12 @@
void LatticeStatsDataProviderBase<T>::reset() {
_minPos.resize(0);
_maxPos.resize(0);
- if (! _progressMeter.null()) {
- _progressMeter->init(_progressMeter->expectedNsteps());
+ /*
+ if (! _progressMeter.null()) {
+ cout << "reset progress" << endl;
+ _progressMeter->init(_progressMeter->expectedNsteps());
}
+ */
}

template <class T>
=======================================
--- /branches/nov14/lattices/LatticeMath/MaskedLatticeStatsDataProvider.h
Mon Feb 16 07:05:15 2015 UTC
+++ /branches/nov14/lattices/LatticeMath/MaskedLatticeStatsDataProvider.h
Tue Mar 3 12:22:11 2015 UTC
@@ -42,7 +42,25 @@
: public LatticeStatsDataProviderBase<T> {
public:

- MaskedLatticeStatsDataProvider(MaskedLattice<T>& lattice);
+ // default constructor. Must set lattice after construction but before
+ // using the object
+
+ MaskedLatticeStatsDataProvider();
+
+ // <src>iteratorLimitBytes</src> is related to the size of the lattice.
+ // If the lattice is greater than this size, then a lattice iterator will
+ // be used to step through the lattice. If less, then all the data in the
+ // values in the lattice are retrieved in a single chunk. The advantage of
+ // the iterator is that less memory is used. The disadvantage is there is
+ // a significant performace cost, so if the lattice is small, it is
better to
+ // get all its values in a single chunk and forgo the iterator. This is
particularly
+ // true when looping for a large number of iterations and creating a
+ // MaskedLatticeStatsDataProvider each loop (in that case, you probably
will want
+ // to create a single object before the loop and use setLattice() to
update
+ // its lattice).
+ MaskedLatticeStatsDataProvider(
+ MaskedLattice<T>& lattice, uInt iteratorLimitBytes=4096*4096
+ );

~MaskedLatticeStatsDataProvider();

@@ -74,6 +92,20 @@
// reset the provider to point to the first data set it manages.
void reset();

+ // set the lattice. Automatically resets the lattice iterator.
+ // <src>iteratorLimitBytes</src> is related to the size of the lattice.
+ // If the lattice is greater than this size, then a lattice iterator will
+ // be used to step through the lattice. If less, then all the data in the
+ // values in the lattice are retrieved in a single chunk. The advantage of
+ // the iterator is that less memory is used. The disadvantage is there is
+ // a significant performace cost, so if the lattice is small, it is
better to
+ // get all its values in a single chunk and forgo the iterator. This is
particularly
+ // true when looping for a large number of iterations and creating a
+ // MaskedLatticeStatsDataProvider each loop (in that case, you probably
will want
+ // to create a single object before the loop and use setLattice() to
update
+ // its lattice).
+ void setLattice(const MaskedLattice<T>& lattice, uInt
iteratorLimitBytes=4096*4096);
+
// <group>
// see base class documentation.
void updateMaxPos(const std::pair<Int64, Int64>& maxpos);
@@ -82,16 +114,16 @@
// </group>

private:
- RO_MaskedLatticeIterator<T> _iter;
+
+ CountedPtr<RO_MaskedLatticeIterator<T> > _iter;
Array<T> _currentSlice;
Array<Bool> _currentMaskSlice;
const T* _currentPtr;
const Bool* _currentMaskPtr;
- Bool _delData, _delMask;
+ Bool _delData, _delMask, _atEnd;

void _freeStorage();

- uInt _nsteps() const;
};

}
=======================================
--- /branches/nov14/lattices/LatticeMath/MaskedLatticeStatsDataProvider.tcc
Mon Feb 16 07:05:15 2015 UTC
+++ /branches/nov14/lattices/LatticeMath/MaskedLatticeStatsDataProvider.tcc
Tue Mar 3 12:22:11 2015 UTC
@@ -30,13 +30,21 @@
#include <casacore/lattices/LatticeMath/MaskedLatticeStatsDataProvider.h>

namespace casacore {
+
+template <class T>
+MaskedLatticeStatsDataProvider<T>::MaskedLatticeStatsDataProvider()
+ : LatticeStatsDataProviderBase<T>(),
+ _iter(), /* _ary(), _mask(), */ _currentSlice(), _currentMaskSlice(),
+ _currentPtr(0), _currentMaskPtr(0), _delData(False), _delMask(False),
_atEnd(False) {}

template <class T>
MaskedLatticeStatsDataProvider<T>::MaskedLatticeStatsDataProvider(
- MaskedLattice<T>& lattice
+ MaskedLattice<T>& lattice, uInt
) : LatticeStatsDataProviderBase<T>(),
- _iter(RO_MaskedLatticeIterator<T>(lattice)), _currentSlice(),
_currentMaskSlice(),
- _currentPtr(0), _currentMaskPtr(0), _delData(False), _delMask(False) {}
+ _iter(), _currentSlice(), _currentMaskSlice(),
+ _currentPtr(0), _currentMaskPtr(0), _delData(False), _delMask(False) {
+ setLattice(lattice);
+}

template <class T>
MaskedLatticeStatsDataProvider<T>::~MaskedLatticeStatsDataProvider() {}
@@ -44,14 +52,22 @@
template <class T>
void MaskedLatticeStatsDataProvider<T>::operator++() {
_freeStorage();
- ++_iter;
+ if (_iter.null()) {
+ _atEnd = True;
+ }
+ else {
+ ++(*_iter);
+ }
this->_updateProgress();
}

template <class T>
uInt MaskedLatticeStatsDataProvider<T>::estimatedSteps() const {
- IPosition lattShape = _iter.latticeShape();
- IPosition cursShape = _iter.cursor().shape();
+ if (_iter.null()) {
+ return 1;
+ }
+ IPosition lattShape = _iter->latticeShape();
+ IPosition cursShape = _iter->cursor().shape();
uInt ndim = lattShape.size();
uInt count = 1;
for (uInt i=0; i<ndim; i++) {
@@ -66,7 +82,13 @@

template <class T>
Bool MaskedLatticeStatsDataProvider<T>::atEnd() const {
- return _iter.atEnd();
+ if (_iter.null()) {
+ return _atEnd;
+
+ }
+ else {
+ return _iter->atEnd();
+ }
}

template <class T>
@@ -77,19 +99,28 @@

template <class T>
uInt64 MaskedLatticeStatsDataProvider<T>::getCount() {
- return _iter.cursor().size();
+ if (_iter.null()) {
+ return _currentSlice.size();
+ }
+ else {
+ return _iter->cursor().size();
+ }
}

template <class T>
const T* MaskedLatticeStatsDataProvider<T>::getData() {
- _currentSlice.assign(_iter.cursor());
+ if (! _iter.null()) {
+ _currentSlice.assign(_iter->cursor());
+ }
_currentPtr = _currentSlice.getStorage(_delData);
return _currentPtr;
}

template <class T>
const Bool* MaskedLatticeStatsDataProvider<T>::getMask() {
- _currentMaskSlice.assign(_iter.getMask());
+ if (! _iter.null()) {
+ _currentMaskSlice.assign(_iter->getMask());
+ }
_currentMaskPtr = _currentMaskSlice.getStorage(_delMask);
return _currentMaskPtr;
}
@@ -102,25 +133,53 @@
template <class T>
void MaskedLatticeStatsDataProvider<T>::reset() {
LatticeStatsDataProviderBase<T>::reset();
- _iter.reset();
+ if (! _iter.null()) {
+ _iter->reset();
+ }
+}
+
+template <class T>
+void MaskedLatticeStatsDataProvider<T>::setLattice(
+ const MaskedLattice<T>& lattice, uInt iteratorLimitBytes
+) {
+ finalize();
+ if (lattice.size() > iteratorLimitBytes/sizeof(T)) {
+ TileStepper stepper(
+ lattice.shape(), lattice.niceCursorShape(
+ lattice.advisedMaxPixels()
+ )
+ );
+ _iter = new RO_MaskedLatticeIterator<T>(lattice, stepper);
+ }
+ else {
+ _iter = NULL;
+ _currentSlice.assign(lattice.get());
+ _currentMaskSlice.assign(lattice.getMask());
+ _atEnd = False;
+ }
}
+

template <class T>
void MaskedLatticeStatsDataProvider<T>::updateMaxPos(
const std::pair<Int64, Int64>& maxpos
) {
- this->_updateMaxPos(
- _iter.position() + toIPositionInArray(maxpos.second,
_currentSlice.shape())
- );
+ IPosition p = toIPositionInArray(maxpos.second, _currentSlice.shape());
+ if (! _iter.null()) {
+ p += _iter->position();
+ }
+ this->_updateMaxPos(p);
}

template <class T>
void MaskedLatticeStatsDataProvider<T>::updateMinPos(
const std::pair<Int64, Int64>& minpos
) {
- this->_updateMinPos(
- _iter.position() + toIPositionInArray(minpos.second,
_currentSlice.shape())
- );
+ IPosition p = toIPositionInArray(minpos.second, _currentSlice.shape());
+ if (! _iter.null()) {
+ p += _iter->position();
+ }
+ this->_updateMinPos(p);
}

template <class T>
@@ -131,6 +190,7 @@
_delMask = False;
}

+/*
template <class T>
uInt MaskedLatticeStatsDataProvider<T>::_nsteps() const {
const IPosition trc = _iter.latticeShape() - 1;
@@ -143,6 +203,7 @@
}
return nsteps;
}
+*/

}

=======================================
--- /branches/nov14/lattices/LatticeMath/test/tLatticeStatistics.cc Mon
Feb 2 10:04:25 2015 UTC
+++ /branches/nov14/lattices/LatticeMath/test/tLatticeStatistics.cc Tue
Mar 3 12:22:11 2015 UTC
@@ -170,6 +170,7 @@
stats.getStatistic(npts, LatticeStatsBase::NPTS, False);
AlwaysAssert(npts.size() == 0, AipsError);
}
+#if 0
{
// using setAlgorithm()
ArrayLattice<Float> latt(data);
@@ -383,6 +384,7 @@
AlwaysAssert(maxPos.size() == 0, AipsError);

}
+#endif
}
catch (const AipsError& x) {
cerr << "aipserror: error " << x.getMesg() << endl;
=======================================
--- /branches/nov14/lattices/Lattices/SubLattice.h Wed Jan 28 10:01:12 2015
UTC
+++ /branches/nov14/lattices/Lattices/SubLattice.h Tue Mar 3 12:22:11 2015
UTC
@@ -287,24 +287,31 @@
}
}

+ // Set the region object using a slicer.
+ // Allows the region to be changed while keeping
+ // the same lattice, so that new SubLattice objects do not have to be
+ // created when one only wants to change the region of interest. Should
+ // only be called when performance is an issue; otherwise, just create
+ // a new SubLattice<T> object.
+ void setRegion (const Slicer& slicer);
+
protected:
- // Set the various pointer needed to construct the object.
- // One of the pointers should be zero.
- // It takes over the pointer and deletes the object in the destructor.
- void setPtr (Lattice<T>* latticePtr,
- MaskedLattice<T>* maskLatPtr,
- Bool writableIfPossible);
-
// Set the region object.
// It also fills in the parent pointer when the SubLattice is taken
// from a MaskedLattice.
// The default region is the entire lattice.
// <group>
void setRegion (const LatticeRegion& region);
- void setRegion (const Slicer& slicer);
void setRegion();
// </group>

+ // Set the various pointers needed to construct the object.
+ // One of the pointers should be zero.
+ // It takes over the pointer and deletes the object in the destructor.
+ void setPtr (Lattice<T>* latticePtr,
+ MaskedLattice<T>* maskLatPtr,
+ Bool writableIfPossible);
+
// Set the axes mapping from the specification.
void setAxesMap (const AxesSpecifier&);

=======================================
--- /branches/nov14/mainpage.dox Wed Feb 11 14:00:31 2015 UTC
+++ /branches/nov14/mainpage.dox Tue Mar 3 12:22:11 2015 UTC
@@ -225,7 +225,7 @@
/// </tr>
/// <tr>
/// <td><a href="group__python.html">python</a></td>
-/// <td>casa boost-python</td>
+/// <td>casa boost-python python numpy</td>
/// </tr>
/// </table>
/// <a
href="http://www.atnf.csiro.au/computing/software/miriad">mirlib</a>
=======================================
--- /branches/nov14/measures/Measures/MDirection.cc Wed Dec 10 08:06:42
2014 UTC
+++ /branches/nov14/measures/Measures/MDirection.cc Tue Mar 3 12:22:11
2015 UTC
@@ -31,6 +31,10 @@
#include <casacore/casa/Arrays/Vector.h>
#include <casacore/casa/Utilities/Register.h>
#include <casacore/measures/Measures/MDirection.h>
+#include <casacore/casa/Quanta/MVAngle.h>
+#include <casacore/casa/Quanta/MVTime.h>
+#include <casacore/casa/Quanta/QMath.h>
+

namespace casacore { //# NAMESPACE CASACORE - BEGIN

@@ -458,6 +462,34 @@
Measure *MDirection::clone() const {
return (new MDirection(*this));
}
+
+String MDirection::toString() const {
+ Quantity longitude = getValue().getLong("deg");
+ Quantity lat = getValue().getLat("deg");
+ Types type = castType(ref.getType());
+ String output;
+ if (
+ type == J2000 || type == JMEAN
+ || type == JTRUE || type == APP
+ || type == B1950 || type == B1950_VLA
+ || type == BMEAN || type == BTRUE
+ ) {
+ String ra = MVTime(longitude).string(MVTime::TIME, 12);
+ String dec = MVAngle(abs(lat)).string(MVAngle::ANGLE_CLEAN, 11);
+ dec.trim();
+ if (lat.getValue() < 0) {
+ dec = "-" + dec;
+ }
+ output = ra + " " + dec;
+ }
+ else {
+ String longStr = MVAngle(longitude).string(MVAngle::ANGLE, 11);
+ String latStr = MVAngle(lat).string(MVAngle::ANGLE, 11);
+ output = longStr + " " + latStr;
+ }
+ output += " " + showType(type);
+ return output;
+}

} //# NAMESPACE CASACORE - END

=======================================
--- /branches/nov14/measures/Measures/MDirection.h Wed Dec 10 08:06:42 2014
UTC
+++ /branches/nov14/measures/Measures/MDirection.h Tue Mar 3 12:22:11 2015
UTC
@@ -388,6 +388,10 @@
virtual Measure *clone() const;
// </group>

+ // Convert to a String in astronomer-friendly format based on
+ // reference frame
+ String toString() const;
+
private:

//# Data
=======================================
--- /branches/nov14/measures/Measures/MeasTableMul.h Wed Dec 10 08:06:42
2014 UTC
+++ /branches/nov14/measures/Measures/MeasTableMul.h Tue Mar 3 12:22:11
2015 UTC
@@ -68,7 +68,7 @@
// removed from the cache, while another thread is still using that
Matrix.
// This assumes that CountedPtr is compiled thread-safe.
//
- // The class provides two virtual function.
+ // The class provides two virtual functions.
// <ul>
// <li> <src>init</src> is called on the first access and makes it
possible
// for the derived class to precompute some variables. In
particular,
@@ -90,6 +90,7 @@
{
public:
MeasTableMul();
+ virtual ~MeasTableMul() {}
void clear();
CountedPtr<Matrix<Double> > getArray (Double time, Double epsilon);
virtual void init() = 0;
=======================================
--- /branches/nov14/ms/CMakeLists.txt Wed Jun 11 12:21:41 2014 UTC
+++ /branches/nov14/ms/CMakeLists.txt Tue Mar 3 12:22:11 2015 UTC
@@ -42,6 +42,7 @@
MeasurementSets/MSScanGram.cc
MeasurementSets/MSState.cc
MeasurementSets/MSReader.cc
+MeasurementSets/MSKeys.cc
MeasurementSets/MSTable2.cc
MeasurementSets/MSWeatherColumns.cc
MeasurementSets/MSArrayParse.cc
@@ -208,6 +209,7 @@
MeasurementSets/MSHistoryEnums.h
MeasurementSets/MSHistoryHandler.h
MeasurementSets/MSIter.h
+MeasurementSets/MSKeys.h
MeasurementSets/MSLister.h
MeasurementSets/MSMainColumns.h
MeasurementSets/MSMainEnums.h
=======================================
--- /branches/nov14/ms/MeasurementSets/MSMetaData.cc Mon Feb 16 07:03:44
2015 UTC
+++ /branches/nov14/ms/MeasurementSets/MSMetaData.cc Tue Mar 3 12:22:11
2015 UTC
@@ -30,20 +30,12 @@
#include <casacore/casa/Arrays/MaskArrMath.h>
#include <casacore/casa/OS/File.h>
#include <casacore/measures/Measures/MeasTable.h>
+#include <casacore/ms/MeasurementSets/MSKeys.h>
#include <casacore/ms/MeasurementSets/MSSpWindowColumns.h>
-#include <casacore/ms/MeasurementSets/MSPointingColumns.h>
#include <casacore/tables/Tables/ArrayColumn.h>
#include <casacore/tables/Tables/ScalarColumn.h>
#include <casacore/tables/TaQL/TableParse.h>
#include <casacore/tables/Tables/TableProxy.h>
-// DEBUG ONLY
-
-/*
-#include <iomanip>
-#include <casacore/casa/Arrays/ArrayIO.h>
-
-#include <casacore/casa/OS/PrecTimer.h>
-*/

#define _ORIGIN "MSMetaData::" + String(__FUNCTION__) + ": "

@@ -53,38 +45,52 @@
: _ms(ms), _cacheMB(0), _maxCacheMB(maxCacheSizeMB), _nStates(0),
_nACRows(0), _nXCRows(0), _nSpw(0), _nFields(0),
_nAntennas(0), _nObservations(0), _nScans(0), _nArrays(0),
- _nrows(0), _nPol(0), _nDataDescIDs(0), _uniqueIntents(),
_scanToSpwsMap(),
- _uniqueScanNumbers(),_uniqueFieldIDs(), _uniqueStateIDs(),
+ _nrows(0), _nPol(0), _nDataDescIDs(0),
+ _scanToSpwsMap(),
+ _scanToDDIDsMap(),
+ _dataDescIDToSpwMap(),
+ _dataDescIDToPolIDMap(), /*_scanToNRowsMap(),*/
+ _fieldToSpwMap(),
+ _scanToStatesMap(), _scanToFieldsMap(),
+ _fieldToStatesMap(), _stateToFieldsMap(),
+ _sourceToFieldsMap(),
+ // _arrayToScansMap(),
+ _antennaNameToIDMap(),
+ _intentToFieldIDMap(),
+ _intentToScansMap(),
+ _scanToTimeRangeMap(),
+ _scanSpwToIntervalMap(),
+ _uniqueIntents(),
+ /*_uniqueScanNumbers(), */ _uniqueFieldIDs(), _uniqueStateIDs(),
_avgSpw(), _tdmSpw(),
_fdmSpw(), _wvrSpw(), _sqldSpw(), _antenna1(), _antenna2(),
_scans(), _fieldIDs(), _stateIDs(), _dataDescIDs(),
_observationIDs(),
- _scanToNACRowsMap(), _scanToNXCRowsMap(),
+ _subScanToNACRowsMap(), _subScanToNXCRowsMap(),
_fieldToNACRowsMap(), _fieldToNXCRowsMap(),
- _dataDescIDToSpwMap(),
_scanToIntentsMap(), _stateToIntentsMap(),
_spwToIntentsMap(),
- _spwInfo(0), _fieldToSpwMap(),
- _spwToFieldIDsMap(0), _spwToScansMap(0),
- _scanToStatesMap(), _scanToFieldsMap(),
+ _spwInfo(0),
+ _spwToFieldIDsMap(), _obsToArraysMap(), _spwToScansMap(),
_fieldToScansMap(),
+
_fieldNames(0),
_antennaNames(0), _observatoryNames(0),
- _antennaNameToIDMap(), _times(),
- _scanToTimesMap(), _intentToFieldIDMap(),
+ _times(),
+ _scanToTimesMap(),
_fieldToTimesMap(), _observatoryPositions(0),
_antennaOffsets(0), _uniqueBaselines(0, 0),
_exposureTime(0), _nUnflaggedACRows(0),
_nUnflaggedXCRows(0), _unflaggedFieldNACRows(),
- _unflaggedFieldNXCRows(), _unflaggedScanNACRows(),
- _unflaggedScanNXCRows(),
+ _unflaggedFieldNXCRows(), _unflaggedSubScanNACRows(),
+ _unflaggedSubScanNXCRows(),
_taqlTableName(
File(ms->tableName()).exists() ? ms->tableName() : "$1"
),
_taqlTempTable(
File(ms->tableName()).exists() ? 0 : 1, ms
- ), _flagsColumn(), _scanToTimeRangeMap(),
- _scanSpwToIntervalMap(), _spwInfoStored(False) {}
+ ), _flagsColumn(),
+ _spwInfoStored(False) {}

MSMetaData::~MSMetaData() {}

@@ -129,9 +135,9 @@
vector<std::set<String> >::iterator sIter = stateToIntentsMap.begin();
for(
Vector<String>::const_iterator curIntentSet=intentSets.begin();
- curIntentSet!=end; curIntentSet++, sIter++
+ curIntentSet!=end; ++curIntentSet, ++sIter
) {
- Vector<String> intents = casacore::stringToVector(*curIntentSet, ',');
+ Vector<String> intents = stringToVector(*curIntentSet, ',');
*sIter = std::set <String>(intents.begin(), intents.end());
uniqueIntents.insert(intents.begin(), intents.end());
}
@@ -143,19 +149,19 @@
uInt count = 0;
for (
vector<std::set<String> >::const_iterator iter=stateToIntentsMap.begin();
- iter!=lastState; iter++, count++
+ iter!=lastState; ++iter, ++count
) {
std::set<String>::const_iterator lastIntent=iter->end();
for (
std::set<String>::const_iterator intent=iter->begin();
- intent!=lastIntent; intent++
+ intent!=lastIntent; ++intent
) {
mysize += intent->size();
}
}
for (
std::set<String>::const_iterator intent=uniqueIntents.begin();
- intent!=lastIntent; intent++
+ intent!=lastIntent; ++intent
) {
mysize += intent->size();
}
@@ -165,8 +171,12 @@
}
}

-std::set<Int> MSMetaData::getScanNumbers() const {
- // This method is responsible for setting _uniqueScanNumbers
+std::set<Int> MSMetaData::getScanNumbers(Int obsID, Int arrayID) const {
+ ArrayKey arrayKey;
+ arrayKey.obsID = obsID;
+ arrayKey.arrayID = arrayID;
+ return _getScanNumbers(arrayKey);
+ /*
if (_uniqueScanNumbers.size() > 0) {
return _uniqueScanNumbers;
}
@@ -176,11 +186,12 @@
_uniqueScanNumbers = myUniqueScans;
}
return myUniqueScans;
+ */
}

uInt MSMetaData::nScans() {
if (_nScans == 0) {
- _nScans = getScanNumbers().size();
+ _nScans = _getScanKeys().size();
}
return _nScans;
}
@@ -204,16 +215,15 @@
return _ms->nrow();
}
uInt MSMetaData::nRows(CorrelationType cType) {
-
if (cType == BOTH) {
return nRows();
}
uInt nACRows, nXCRows;
- CountedPtr<AOSFMapI> scanToNACRowsMap, scanToNXCRowsMap;
- CountedPtr<std::map<Int, uInt> > fieldToNACRowsMap, fieldToNXCRowsMap;
+ CountedPtr<std::map<SubScanKey, uInt> > subScanToNACRowsMap,
subScanToNXCRowsMap;
+ CountedPtr<vector<uInt> > fieldToNACRowsMap, fieldToNXCRowsMap;
_getRowStats(
- nACRows, nXCRows, scanToNACRowsMap,
- scanToNXCRowsMap, fieldToNACRowsMap,
+ nACRows, nXCRows, subScanToNACRowsMap,
+ subScanToNXCRowsMap, fieldToNACRowsMap,
fieldToNXCRowsMap
);
if (cType == AUTO) {
@@ -228,34 +238,40 @@
CorrelationType cType, Int arrayID, Int observationID,
Int scanNumber, Int fieldID
) const {
+ SubScanKey subScanKey;
+ subScanKey.obsID = observationID;
+ subScanKey.arrayID = arrayID;
+ subScanKey.scan = scanNumber;
+ subScanKey.fieldID = fieldID;
+ _checkSubScan(subScanKey);
uInt nACRows, nXCRows;
- CountedPtr<AOSFMapI> scanToNACRowsMap, scanToNXCRowsMap;
- CountedPtr<std::map<Int, uInt> > fieldToNACRowsMap, fieldToNXCRowsMap;
+ CountedPtr<std::map<SubScanKey, uInt> > subScanToNACRowsMap,
subScanToNXCRowsMap;
+ CountedPtr<vector<uInt> > fieldToNACRowsMap, fieldToNXCRowsMap;
_getRowStats(
- nACRows, nXCRows, scanToNACRowsMap,
- scanToNXCRowsMap, fieldToNACRowsMap,
+ nACRows, nXCRows, subScanToNACRowsMap,
+ subScanToNXCRowsMap, fieldToNACRowsMap,
fieldToNXCRowsMap
);
if (cType == AUTO) {
- return (*scanToNACRowsMap)[arrayID][observationID][scanNumber][fieldID];
+ return (*subScanToNACRowsMap)[subScanKey];
}
else if (cType == CROSS) {
- return (*scanToNXCRowsMap)[arrayID][observationID][scanNumber][fieldID];
+ return (*subScanToNXCRowsMap)[subScanKey];
}
else {
- return (*scanToNACRowsMap)[arrayID][observationID][scanNumber][fieldID]
- + (*scanToNXCRowsMap)[arrayID][observationID][scanNumber][fieldID];
+ return (*subScanToNACRowsMap)[subScanKey]
+ + (*subScanToNXCRowsMap)[subScanKey];
}
-
}

-uInt MSMetaData::nRows(CorrelationType cType, Int fieldID) const {
+uInt MSMetaData::nRows(CorrelationType cType, uInt fieldID) const {
+ _checkField(fieldID);
uInt nACRows, nXCRows;
- CountedPtr<AOSFMapI> scanToNACRowsMap, scanToNXCRowsMap;
- CountedPtr<std::map<Int, uInt> > fieldToNACRowsMap, fieldToNXCRowsMap;
+ CountedPtr<std::map<SubScanKey, uInt> > subScanToNACRowsMap,
subScanToNXCRowsMap;
+ CountedPtr<vector<uInt> > fieldToNACRowsMap, fieldToNXCRowsMap;
_getRowStats(
- nACRows, nXCRows, scanToNACRowsMap,
- scanToNXCRowsMap, fieldToNACRowsMap,
+ nACRows, nXCRows, subScanToNACRowsMap,
+ subScanToNXCRowsMap, fieldToNACRowsMap,
fieldToNXCRowsMap
);
if (cType == AUTO) {
@@ -271,11 +287,11 @@

Double MSMetaData::nUnflaggedRows() const {
Double nACRows, nXCRows;
- CountedPtr<AOSFMapD> scanToNACRowsMap, scanToNXCRowsMap;
- CountedPtr<std::map<Int, Double> > fieldToNACRowsMap, fieldToNXCRowsMap;
+ CountedPtr<std::map<SubScanKey, Double> > subScanToNACRowsMap,
subScanToNXCRowsMap;
+ CountedPtr<vector<Double> > fieldToNACRowsMap, fieldToNXCRowsMap;
_getUnflaggedRowStats(
- nACRows, nXCRows, scanToNACRowsMap,
- scanToNXCRowsMap, fieldToNACRowsMap,
+ nACRows, nXCRows, subScanToNACRowsMap,
+ subScanToNXCRowsMap, fieldToNACRowsMap,
fieldToNXCRowsMap
);
return nACRows + nXCRows;
@@ -285,11 +301,11 @@
return nUnflaggedRows();
}
Double nACRows, nXCRows;
- CountedPtr<AOSFMapD> scanToNACRowsMap, scanToNXCRowsMap;
- CountedPtr<std::map<Int, Double> > fieldToNACRowsMap, fieldToNXCRowsMap;
+ CountedPtr<std::map<SubScanKey, Double> > subScanToNACRowsMap,
subScanToNXCRowsMap;
+ CountedPtr<vector<Double> > fieldToNACRowsMap, fieldToNXCRowsMap;
_getUnflaggedRowStats(
- nACRows, nXCRows, scanToNACRowsMap,
- scanToNXCRowsMap, fieldToNACRowsMap,
+ nACRows, nXCRows, subScanToNACRowsMap,
+ subScanToNXCRowsMap, fieldToNACRowsMap,
fieldToNXCRowsMap
);
if (cType == AUTO) {
@@ -301,36 +317,42 @@
}

Double MSMetaData::nUnflaggedRows(
- CorrelationType cType, Int arrayID, Int observationID,
- Int scanNumber, Int fieldID
+ CorrelationType cType, Int arrayID, uInt observationID,
+ Int scanNumber, uInt fieldID
) const {
+ SubScanKey subScanKey;
+ subScanKey.obsID = observationID;
+ subScanKey.arrayID = arrayID;
+ subScanKey.scan = scanNumber;
+ subScanKey.fieldID = fieldID;
+ _checkSubScan(subScanKey);
Double nACRows, nXCRows;
- CountedPtr<AOSFMapD> scanToNACRowsMap, scanToNXCRowsMap;
- CountedPtr<std::map<Int, Double> > fieldToNACRowsMap, fieldToNXCRowsMap;
+ CountedPtr<std::map<SubScanKey, Double> > subScanToNACRowsMap,
subScanToNXCRowsMap;
+ CountedPtr<vector<Double> > fieldToNACRowsMap, fieldToNXCRowsMap;
_getUnflaggedRowStats(
- nACRows, nXCRows, scanToNACRowsMap,
- scanToNXCRowsMap, fieldToNACRowsMap,
+ nACRows, nXCRows, subScanToNACRowsMap,
+ subScanToNXCRowsMap, fieldToNACRowsMap,
fieldToNXCRowsMap
);
if (cType == AUTO) {
- return (*scanToNACRowsMap)[arrayID][observationID][scanNumber][fieldID];
+ return (*subScanToNACRowsMap)[subScanKey];
}
else if (cType == CROSS) {
- return (*scanToNXCRowsMap)[arrayID][observationID][scanNumber][fieldID];
+ return (*subScanToNXCRowsMap)[subScanKey];
}
else {
- return (*scanToNACRowsMap)[arrayID][observationID][scanNumber][fieldID]
- + (*scanToNXCRowsMap)[arrayID][observationID][scanNumber][fieldID];
+ return (*subScanToNACRowsMap)[subScanKey]
+ + (*subScanToNXCRowsMap)[subScanKey];
}
}

Double MSMetaData::nUnflaggedRows(CorrelationType cType, Int fieldID)
const {
Double nACRows, nXCRows;
- CountedPtr<AOSFMapD> scanToNACRowsMap, scanToNXCRowsMap;
- CountedPtr<std::map<Int, Double> > fieldToNACRowsMap, fieldToNXCRowsMap;
+ CountedPtr<std::map<SubScanKey, Double> > subScanToNACRowsMap,
subScanToNXCRowsMap;
+ CountedPtr<vector<Double> > fieldToNACRowsMap, fieldToNXCRowsMap;
_getUnflaggedRowStats(
- nACRows, nXCRows, scanToNACRowsMap,
- scanToNXCRowsMap, fieldToNACRowsMap,
+ nACRows, nXCRows, subScanToNACRowsMap,
+ subScanToNXCRowsMap, fieldToNACRowsMap,
fieldToNXCRowsMap
);
if (cType == AUTO) {
@@ -344,6 +366,7 @@
}
}

+/*
void MSMetaData::_getRowStats(
uInt& nACRows, uInt& nXCRows,
CountedPtr<AOSFMapI>& scanToNACRowsMap,
@@ -372,10 +395,10 @@
myScanToNXCRowsMap, myFieldToNACRowsMap,
myFieldToNXCRowsMap
);
- scanToNACRowsMap.reset(myScanToNACRowsMap);
- scanToNXCRowsMap.reset(myScanToNXCRowsMap);
- fieldToNACRowsMap.reset(myFieldToNACRowsMap);
- fieldToNXCRowsMap.reset(myFieldToNXCRowsMap);
+ scanToNACRowsMap = myScanToNACRowsMap;
+ scanToNXCRowsMap = myScanToNXCRowsMap;
+ fieldToNACRowsMap = myFieldToNACRowsMap;
+ fieldToNXCRowsMap = myFieldToNXCRowsMap;

Float newSize = _cacheMB + sizeof(Int)*(
2 + 2*scanToNACRowsMap->size()
@@ -394,6 +417,9 @@
_fieldToNXCRowsMap = fieldToNXCRowsMap;
}
}
+*/
+
+/*

void MSMetaData::_getRowStats(
uInt& nACRows, uInt& nXCRows,
@@ -428,19 +454,19 @@
scanToNXCRowsMap->clear();
for (
std::set<Int>::const_iterator arrNum=uniqueArrIDs.begin();
- arrNum!=lastUniqueArrID; arrNum++
+ arrNum!=lastUniqueArrID; ++arrNum
) {
for (
std::set<Int>::const_iterator obsNum=uniqueObsIDs.begin();
- obsNum!=lastUniqueObsID; obsNum++
+ obsNum!=lastUniqueObsID; ++obsNum
) {
for (
std::set<Int>::const_iterator scanNum=scanNumbers.begin();
- scanNum!=lastScan; scanNum++
+ scanNum!=lastScan; ++scanNum
) {
for (
std::set<Int>::const_iterator fieldNum=uniqueFieldIDs.begin();
- fieldNum!=lastUniqueFieldID; fieldNum++
+ fieldNum!=lastUniqueFieldID; ++fieldNum
) {
(*scanToNACRowsMap)[*arrNum][*obsNum][*scanNum][*fieldNum] = 0;
(*scanToNXCRowsMap)[*arrNum][*obsNum][*scanNum][*fieldNum] = 0;
@@ -450,7 +476,7 @@
}
for (
std::set<Int>::const_iterator fieldNum=uniqueFieldIDs.begin();
- fieldNum!=lastUniqueFieldID; fieldNum++
+ fieldNum!=lastUniqueFieldID; ++fieldNum
) {
(*fieldToNACRowsMap)[*fieldNum] = 0;
(*fieldToNXCRowsMap)[*fieldNum] = 0;
@@ -472,23 +498,128 @@

}
if (*a1Iter == *a2Iter) {
- nACRows++;
- (*scanToNACRowsMap)[*arIter][*oIter][*sIter][*fIter]++;
- (*fieldToNACRowsMap)[*fIter]++;
+ ++nACRows;
+ ++(*scanToNACRowsMap)[*arIter][*oIter][*sIter][*fIter];
+ ++(*fieldToNACRowsMap)[*fIter];
}
else {
- nXCRows++;
- (*scanToNXCRowsMap)[*arIter][*oIter][*sIter][*fIter]++;
- (*fieldToNXCRowsMap)[*fIter]++;
+ ++nXCRows;
+ ++(*scanToNXCRowsMap)[*arIter][*oIter][*sIter][*fIter];
+ ++(*fieldToNXCRowsMap)[*fIter];
}
- a1Iter++;
- a2Iter++;
- sIter++;
- fIter++;
- arIter++;
- oIter++;
+ ++a1Iter;
+ ++a2Iter;
+ ++sIter;
+ ++fIter;
+ ++arIter;
+ ++oIter;
}
}
+*/
+
+void MSMetaData::_getRowStats(
+ uInt& nACRows, uInt& nXCRows,
+ std::map<SubScanKey, uInt>*& subScanToNACRowsMap,
+ std::map<SubScanKey, uInt>*& subScanToNXCRowsMap,
+ vector<uInt>*& fieldToNACRowsMap,
+ vector<uInt>*& fieldToNXCRowsMap
+) const {
+ nACRows = 0;
+ nXCRows = 0;
+ subScanToNACRowsMap = new std::map<SubScanKey, uInt>();
+ subScanToNXCRowsMap = new std::map<SubScanKey, uInt>();
+ uInt myNFields = nFields();
+ fieldToNACRowsMap = new vector<uInt>(myNFields, 0);
+ fieldToNXCRowsMap = new vector<uInt>(myNFields, 0);
+ std::set<SubScanKey> subScanKeys = _getSubScanKeys();
+ std::set<SubScanKey>::const_iterator subIter = subScanKeys.begin();
+ std::set<SubScanKey>::const_iterator subEnd = subScanKeys.end();
+ while (subIter != subEnd) {
+ (*subScanToNACRowsMap)[*subIter] = 0;
+ (*subScanToNXCRowsMap)[*subIter] = 0;
+ ++subIter;
+ }
+ CountedPtr<Vector<Int> > ant1, ant2;
+ _getAntennas(ant1, ant2);
+ CountedPtr<Vector<Int> > scans = _getScans();
+ CountedPtr<Vector<Int> > fieldIDs = _getFieldIDs();
+ CountedPtr<Vector<Int> > obsIDs = _getObservationIDs();
+ CountedPtr<Vector<Int> > arrIDs = _getArrayIDs();
+ Vector<Int>::const_iterator aEnd = ant1->end();
+ Vector<Int>::const_iterator a1Iter = ant1->begin();
+ Vector<Int>::const_iterator a2Iter = ant2->begin();
+ Vector<Int>::const_iterator sIter = scans->begin();
+ Vector<Int>::const_iterator fIter = fieldIDs->begin();
+ Vector<Int>::const_iterator oIter = obsIDs->begin();
+ Vector<Int>::const_iterator arIter = arrIDs->begin();
+ SubScanKey subScanKey;
+ while (a1Iter!=aEnd) {
+ subScanKey.obsID = *oIter;
+ subScanKey.arrayID = *arIter;
+ subScanKey.scan = *sIter;
+ subScanKey.fieldID = *fIter;
+ if (*a1Iter == *a2Iter) {
+ ++nACRows;
+ ++(*subScanToNACRowsMap)[subScanKey];
+ ++(*fieldToNACRowsMap)[*fIter];
+ }
+ else {
+ ++nXCRows;
+ ++(*subScanToNXCRowsMap)[subScanKey];
+ ++(*fieldToNXCRowsMap)[*fIter];
+ }
+ ++a1Iter;
+ ++a2Iter;
+ ++sIter;
+ ++fIter;
+ ++arIter;
+ ++oIter;
+ }
+}
+
+void MSMetaData::_getRowStats(
+ uInt& nACRows, uInt& nXCRows,
+ CountedPtr<std::map<SubScanKey, uInt> >& scanToNACRowsMap,
+ CountedPtr<std::map<SubScanKey, uInt> >& scanToNXCRowsMap,
+ CountedPtr<vector<uInt> >& fieldToNACRowsMap,
+ CountedPtr<vector<uInt> >& fieldToNXCRowsMap
+) const {
+ // this method is responsible for setting _nACRows, _nXCRows,
_subScanToNACRowsMap,
+ // _subScanToNXCRowsMap, _fieldToNACRowsMap, _fieldToNXCRowsMap
+ if (_nACRows > 0 || _nXCRows > 0) {
+ nACRows = _nACRows;
+ nXCRows = _nXCRows;
+ scanToNACRowsMap = _subScanToNACRowsMap;
+ scanToNXCRowsMap = _subScanToNXCRowsMap;
+ fieldToNACRowsMap = _fieldToNACRowsMap;
+ fieldToNXCRowsMap = _fieldToNXCRowsMap;
+ return;
+ }
+
+ std::map<SubScanKey, uInt> *myScanToNACRowsMap, *myScanToNXCRowsMap;
+ vector<uInt> *myFieldToNACRowsMap, *myFieldToNXCRowsMap;
+ _getRowStats(
+ nACRows, nXCRows, myScanToNACRowsMap,
+ myScanToNXCRowsMap, myFieldToNACRowsMap,
+ myFieldToNXCRowsMap
+ );
+ scanToNACRowsMap = myScanToNACRowsMap;
+ scanToNXCRowsMap = myScanToNXCRowsMap;
+ fieldToNACRowsMap = myFieldToNACRowsMap;
+ fieldToNXCRowsMap = myFieldToNXCRowsMap;
+ uInt size = 2*(
+ sizeof(uInt) + _sizeof(*scanToNACRowsMap)
+ + _sizeof(*fieldToNACRowsMap)
+ );
+ if (_cacheUpdated(size)) {
+ _nACRows = nACRows;
+ _nXCRows = nXCRows;
+ _subScanToNACRowsMap = scanToNACRowsMap;
+ _subScanToNXCRowsMap = scanToNXCRowsMap;
+ _fieldToNACRowsMap = fieldToNACRowsMap;
+ _fieldToNXCRowsMap = fieldToNXCRowsMap;
+ }
+}

void MSMetaData::_getAntennas(
CountedPtr<Vector<Int> >& ant1,
@@ -508,7 +639,6 @@
ROScalarColumn<Int> ant2Col(*_ms, ant2ColName);
Vector<Int> a2 = ant2Col.getColumn();

- //MSMetaData::_getAntennas(a1, a2, *_ms);
ant1.reset(new Vector<Int>(a1));
ant2.reset(new Vector<Int>(a2));

@@ -537,10 +667,10 @@
static const String obsColName =
MeasurementSet::columnName(MSMainEnums::OBSERVATION_ID);
CountedPtr<Vector<Int> > obsIDs(
new Vector<Int>(ROScalarColumn<Int>(*_ms, obsColName).getColumn())
- );
+ );
if (_cacheUpdated(sizeof(Int)*obsIDs->size())) {
_observationIDs = obsIDs;
- }
+ }
return obsIDs;
}

@@ -557,6 +687,23 @@
}
return arrIDs;
}
+
+/*
+CountedPtr<Vector<Int> > MSMetaData::_getDataDescIDs() const {
+ if (! _dataDescIDs.null() && ! _dataDescIDs->empty()) {
+ return _dataDescIDs;
+ }
+ static const String ddColName =
MeasurementSet::columnName(MSMainEnums::DATA_DESC_ID);
+ ROScalarColumn<Int> ddCol(*_ms, ddColName);
+ CountedPtr<Vector<Int> > dataDescIDs(
+ new Vector<Int>(ddCol.getColumn())
+ );
+ if (_cacheUpdated(sizeof(Int)*dataDescIDs->size())) {
+ _dataDescIDs = dataDescIDs;
+ }
+ return dataDescIDs;
+}
+*/

CountedPtr<Vector<Int> > MSMetaData::_getFieldIDs() const {
if (_fieldIDs && ! _fieldIDs->empty()) {
@@ -571,6 +718,78 @@
}
return fields;
}
+
+/*
+CountedPtr<Vector<Int> > MSMetaData::_getObservationIDs() const {
+ if (! _observationIDs.null() && _observationIDs->size() > 0) {
+ return _observationIDs;
+ }
+ static const String obsColName =
MeasurementSet::columnName(MSMainEnums::OBSERVATION_ID);
+ CountedPtr<Vector<Int> > obsIDs(
+ new Vector<Int>(ROScalarColumn<Int>(*_ms, obsColName).getColumn())
+ );
+ if (_cacheUpdated(sizeof(Int)*obsIDs->size())) {
+ _observationIDs = obsIDs;
+ }
+ return obsIDs;
+}
+*/
+
+/*
+std::map<Int, std::set<Int> > MSMetaData::_getScanToAntennasMap() const {
+ // responsible for setting _scanToAntennasMap
+ if (! _scanToAntennasMap.empty()) {
+ return _scanToAntennasMap;
+ }
+ Vector<Int> scans = *_getScans();
+ CountedPtr<Vector<Int> > ant1, ant2;
+ _getAntennas(ant1, ant2);
+ Vector<Int>::const_iterator sIter = scans.begin();
+ Vector<Int>::const_iterator sEnd = scans.end();
+ Vector<Int>::const_iterator a1Iter = ant1->begin();
+ Vector<Int>::const_iterator a2Iter = ant2->begin();
+ std::map<Int, std::set<Int> > mymap;
+ while (sIter != sEnd) {
+ mymap[*sIter].insert(*a1Iter);
+ mymap[*sIter].insert(*a2Iter);
+ ++sIter;
+ ++a1Iter;
+ ++a2Iter;
+ }
+ if (_cacheUpdated(_sizeof(mymap))) {
+ _scanToAntennasMap = mymap;
+ }
+ return mymap;
+}
+*/
+
+/*
+std::map<Int, uInt> MSMetaData::_getScansToNRowsMap() const {
+ // this method is responsible for setting _scanToNRowsMap
+ if (! _scanToNRowsMap.empty()) {
+ return _scanToNRowsMap;
+ }
+ std::map<Int, uInt> mymap;
+ std::set<Int> scanNumbers = getScanNumbers();
+ std::set<Int>::const_iterator nIter = scanNumbers.begin();
+ std::set<Int>::const_iterator nEnd = scanNumbers.end();
+ while (nIter != nEnd) {
+ mymap[*nIter] = 0;
+ ++nIter;
+ }
+ Vector<Int> scans = *_getScans();
+ Vector<Int>::const_iterator sIter = scans.begin();
+ Vector<Int>::const_iterator sEnd = scans.end();
+ while (sIter != sEnd) {
+ ++mymap[*sIter];
+ ++sIter;
+ }
+ if (_cacheUpdated(2*sizeof(Int)*(mymap.size()))) {
+ _scanToNRowsMap = mymap;
+ }
+ return mymap;
+}
+*/

CountedPtr<Vector<Int> > MSMetaData::_getStateIDs() const {
if (_stateIDs && _stateIDs->size() > 0) {
@@ -609,75 +828,125 @@
return dataDescIDs;
}

-std::set<Int> MSMetaData::getScansForState(const Int stateID) {
+std::set<Int> MSMetaData::getScansForState(
+ Int stateID, Int obsID, Int arrayID
+) {
if (! _hasStateID(stateID)) {
return std::set<Int>();
}
- std::set<Int> uniqueScans;
- std::map<Int, std::set<Int> > myScanToStatesMap = _getScanToStatesMap();
- CountedPtr<Vector<Int> > scans = _getScans();
- uniqueScans.insert(scans->begin(), scans->end());
+ std::map<ScanKey, std::set<Int> > myScanToStatesMap =
_getScanToStatesMap();
+ ArrayKey arrayKey;
+ arrayKey.obsID = obsID;
+ arrayKey.arrayID = arrayID;
+ std::set<ScanKey> scanKeys = _getScanKeys(arrayKey);
+ //CountedPtr<Vector<Int> > scans = _getScans();
+ std::set<ScanKey>::const_iterator iter = scanKeys.begin();
+ std::set<ScanKey>::const_iterator end = scanKeys.end();
+ std::set<Int> stateIDs, scansForState;
+ while (iter != end) {
+ stateIDs = myScanToStatesMap[*iter];
+ if (stateIDs.find(stateID) != stateIDs.end()) {
+ scansForState.insert(iter->scan);
+ }
+ ++iter;
+ }
+
+
+
+ // uniqueScans.insert(scans->begin(), scans->end());
+ /*
std::set<Int>::const_iterator lastScan = uniqueScans.end();
std::set<Int> scansForState;
for (
std::set<Int>::const_iterator scanNum=uniqueScans.begin();
- scanNum!=lastScan; scanNum++
+ scanNum!=lastScan; ++scanNum
) {
std::set<Int> statesSet = myScanToStatesMap.find(*scanNum)->second;
if (statesSet.find(stateID) != statesSet.end()) {
scansForState.insert(*scanNum);
}
}
+ */
return scansForState;
}

-std::map<Int, std::set<Int> > MSMetaData::_getScanToStatesMap() const {
+std::map<ScanKey, std::set<Int> > MSMetaData::_getScanToStatesMap() const {
if (! _scanToStatesMap.empty()) {
return _scanToStatesMap;
}
- std::map<Int, std::set<Int> > myScanToStatesMap;
- uInt mySize = 0;
-
+ std::map<ScanKey, std::set<Int> > myScanToStatesMap;
if (nStates() == 0) {
std::set<Int> empty;
- std::set<Int> uniqueScans = getScanNumbers();
- std::set<Int>::const_iterator end = uniqueScans.end();
+ std::set<ScanKey> scanKeys = this->_getScanKeys();
+ //std::set<SubScanKey> subScanKeys;
+ //_getSubScanKeys(subScanKeys, scanKeys);
+ std::set<ScanKey>::const_iterator end = scanKeys.end();
for (
- std::set<Int>::const_iterator scanNum=uniqueScans.begin();
- scanNum!=end; scanNum++
+ std::set<ScanKey>::const_iterator scanKey=scanKeys.begin();
+ scanKey!=end; ++scanKey
) {
- myScanToStatesMap[*scanNum] = empty;
+ myScanToStatesMap[*scanKey] = empty;
}
}
else {
+ map<SubScanKey, SubScanProperties> subScanProps =
_getSubScanProperties();
+ //map<ScanKey, ScanProperties> scanProps;
+ //map<ArrayKey, ArrayProperties> arrayProps;
+
+ //_getSubScanProperties(subScanProps, scanProps, arrayProps);
+ map<SubScanKey, SubScanProperties>::const_iterator iter =
subScanProps.begin();
+ map<SubScanKey, SubScanProperties>::const_iterator end =
subScanProps.end();
+ ScanKey key;
+ while (iter != end) {
+ SubScanKey subKey = iter->first;
+ key.obsID = subKey.obsID;
+ key.arrayID = subKey.arrayID;
+ key.scan = subKey.scan;
+ SubScanProperties subProps = iter->second;
+ myScanToStatesMap[scanKey(subKey)].insert(subProps.stateIDs.begin(),
subProps.stateIDs.end());
+ ++iter;
+ }
+
+ /*
+ const Vector<Int> obs = *_getObservationIDs();
+ const Vector<Int> array = *_getArrayIDs();
const Vector<Int> scans = *(_getScans());
const Vector<Int> states = *(_getStateIDs());
Vector<Int>::const_iterator curScan = scans.begin();
Vector<Int>::const_iterator lastScan = scans.end();
+ Vector<Int>::const_iterator curArray = array.begin();
+ Vector<Int>::const_iterator curObs = obs.begin();
Vector<Int>::const_iterator curStateID = states.begin();
- //std::map<Int, std::set<Int> > myScanToStatesMap;
while (curScan != lastScan) {
- myScanToStatesMap[*curScan].insert(*curStateID);
- curScan++;
- curStateID++;
+ ScanKey key;
+ key.obsID = *curObs;
+ key.arrayID = *curArray;
+ key.scan = *curScan;
+ myScanToStatesMap[key].insert(*curStateID);
+ ++curScan;
+ ++curStateID;
+ ++curArray;
+ ++curObs;
}
+ */
}
- std::map<Int, std::set<Int> >::const_iterator end =
myScanToStatesMap.end();
+ std::map<ScanKey, std::set<Int> >::const_iterator end =
myScanToStatesMap.end();
+ uInt mySize = sizeof(ScanKey)*myScanToStatesMap.size();
for (
- std::map<Int, std::set<Int> >::const_iterator
iter=myScanToStatesMap.begin();
- iter!=end; iter++
+ std::map<ScanKey, std::set<Int> >::const_iterator
iter=myScanToStatesMap.begin();
+ iter!=end; ++iter
) {
- mySize += iter->second.size() + 1;
+ mySize += sizeof(Int)*iter->second.size();
}
- if (_cacheUpdated(mySize*sizeof(Int)/1e6)) {
+ if (_cacheUpdated(mySize)) {
_scanToStatesMap = myScanToStatesMap;
}
return myScanToStatesMap;
}

void MSMetaData::_getScansAndIntentsMaps(
- std::map<Int, std::set<String> >& scanToIntentsMap,
- std::map<String, std::set<Int> >& intentToScansMap
+ std::map<ScanKey, std::set<String> >& scanToIntentsMap,
+ std::map<String, std::set<ScanKey> >& intentToScansMap
) const {
// This method is responsible for setting _scanToIntentsMap and
_intentToScansMap
if (! _scanToIntentsMap.empty() && ! _intentToScansMap.empty()) {
@@ -690,27 +959,27 @@
_getStateToIntentsMap(
stateToIntentsMap, uniqueIntents
);
- std::map<Int, std::set<Int> > scanToStatesMap = _getScanToStatesMap();
- std::map<Int, std::set<Int> >::const_iterator end = scanToStatesMap.end();
+ std::map<ScanKey, std::set<Int> > scanToStatesMap = _getScanToStatesMap();
+ std::map<ScanKey, std::set<Int> >::const_iterator end =
scanToStatesMap.end();
std::set<Int> states;
std::set<String> intents;
for (
- std::map<Int, std::set<Int> >::const_iterator
iter=scanToStatesMap.begin();
- iter!=end; iter++
+ std::map<ScanKey, std::set<Int> >::const_iterator
iter=scanToStatesMap.begin();
+ iter!=end; ++iter
) {
- uInt scan = iter->first;
+ ScanKey scan = iter->first;
states = iter->second;
std::set<Int>::const_iterator endState = states.end();
for (
std::set<Int>::const_iterator myState=states.begin();
- myState!=endState; myState++
+ myState!=endState; ++myState
) {
intents = stateToIntentsMap[*myState];
scanToIntentsMap[scan].insert(intents.begin(), intents.end());
std::set<String>::const_iterator endIntent = intents.end();
for (
std::set<String>::const_iterator myIntent=intents.begin();
- myIntent!=endIntent; myIntent++
+ myIntent!=endIntent; ++myIntent
) {
intentToScansMap[*myIntent].insert(scan);
}
@@ -722,35 +991,91 @@
}
}

+template <class T>
+uInt MSMetaData::_sizeof(const std::map<T, std::set<String> >& m) {
+ uInt size = sizeof(T) * m.size();
+ typename std::map<T, std::set<String> >::const_iterator iter = m.begin();
+ typename std::map<T, std::set<String> >::const_iterator end = m.end();
+ while (iter != end) {
+ std::set<String>::const_iterator end2 = iter->second.end();
+ for (
+ std::set<String>::const_iterator iter2=iter->second.begin();
+ iter2!=end2; ++iter2
+ ) {
+ size += iter2->size();
+ }
+ ++iter;
+ }
+ return size;
+}
+
+/*
uInt MSMetaData::_sizeof(const std::map<Int, std::set<String> >& m) {
uInt size = sizeof(Int) * m.size();
std::map<Int, std::set<String> >::const_iterator end = m.end();
for (
std::map<Int, std::set<String> >::const_iterator iter=m.begin();
- iter!=end; iter++
+ iter!=end; ++iter
) {
std::set<String>::const_iterator end2 = iter->second.end();
for (
std::set<String>::const_iterator iter2=iter->second.begin();
- iter2!=end2; iter2++
+ iter2!=end2; ++iter2
) {
size += iter2->size();
}
}
return size;
}
+
+uInt MSMetaData::_sizeof(const std::map<ScanKey, std::set<String> >& m) {
+ uInt size = sizeof(ScanKey)*m.size();
+ std::map<ScanKey, std::set<String> >::const_iterator iter = m.begin();
+ std::map<ScanKey, std::set<String> >::const_iterator end = m.end();
+ while (iter != end) {
+ std::set<String> strings = iter->second;
+ std::set<String>::const_iterator sIter = strings.begin();
+ std::set<String>::const_iterator sEnd = strings.end();
+ while (sIter != sEnd) {
+ size += sIter->size();
+ ++sIter;
+ }
+ ++iter;
+ }
+ return size;
+}
+*/
+
+template <class T, class U>
+uInt MSMetaData::_sizeof(const std::map<T, std::set<U> >& m) {
+ uInt size = sizeof(T)*m.size();
+ typename std::map<T, std::set<U> >::const_iterator iter = m.begin();
+ typename std::map<T, std::set<U> >::const_iterator end = m.end();
+ uInt nElements = 0;
+ while (iter != end) {
+ nElements += iter->second.size();
+ ++iter;
+ }
+ size += sizeof(U)*nElements;
+ return size;
+}
+
+template <class T, class U>
+uInt MSMetaData::_sizeof(const std::map<T, U>& m) {
+ return m.size()*(sizeof(T) + sizeof(U));
+}

uInt MSMetaData::_sizeof(const vector<std::set<String> >& m) {
uInt size = sizeof(Int) * m.size();
vector<std::set<String> >::const_iterator end = m.end();
for (
vector<std::set<String> >::const_iterator iter=m.begin();
- iter!=end; iter++
+ iter!=end; ++iter
) {
std::set<String>::const_iterator end2 = iter->end();
for (
std::set<String>::const_iterator iter2=iter->begin();
- iter2!=end2; iter2++
+ iter2!=end2; ++iter2
) {
size += iter2->size();
}
@@ -763,12 +1088,17 @@
uInt size = 0;
for (
vector<String>::const_iterator iter=m.begin();
- iter!=end; iter++
+ iter!=end; ++iter
) {
size += iter->length();
}
return size;
}
+
+template <class T>
+uInt MSMetaData::_sizeof(const vector<T>& v) {
+ return v.size()*sizeof(T);
+}

uInt MSMetaData::_sizeof(const Quantum<Vector<Double> >& m) {
return (sizeof(Double)+10)*m.getValue().size();
@@ -780,7 +1110,7 @@
typename std::map<String, std::set<T> >::const_iterator end = m.end();
for (
typename std::map<String, std::set<T> >::const_iterator iter=m.begin();
- iter!=end; iter++
+ iter!=end; ++iter
) {
size += iter->first.size();
setssize += iter->second.size();
@@ -796,20 +1126,20 @@
uInt qsize = 20;
for (
vector<std::map<Int, Quantity> >::const_iterator iter = m.begin();
- iter!=end; iter++
+ iter!=end; ++iter
) {
size += iter->size()*(2*intsize + qsize);
}
return size;
}
-
+/*
uInt MSMetaData::_sizeof(const std::map<Double, std::set<Int> >& m) {
uInt setssize = 0;
uInt size = sizeof(Double) * m.size();
std::map<Double, std::set<Int> >::const_iterator end = m.end();
for (
std::map<Double, std::set<Int> >::const_iterator iter=m.begin();
- iter!=end; iter++
+ iter!=end; ++iter
) {
setssize += iter->second.size();
}
@@ -823,18 +1153,19 @@
std::map<Int, std::set<Double> >::const_iterator end = m.end();
for (
std::map<Int, std::set<Double> >::const_iterator iter=m.begin();
- iter!=end; iter++
+ iter!=end; ++iter
) {
setssize += iter->second.size();
}
size += sizeof(Double) * setssize;
return size;
}
+*/

-std::set<String> MSMetaData::getIntentsForScan(const Int scan) const {
- _checkScan(scan, getScanNumbers());
***The diff for this file has been truncated for email.***
=======================================
--- /branches/nov14/ms/MeasurementSets/MSMetaData.h Mon Feb 16 07:03:44
2015 UTC
+++ /branches/nov14/ms/MeasurementSets/MSMetaData.h Tue Mar 3 12:22:11
2015 UTC
@@ -23,7 +23,6 @@
//# 520 Edgemont Road
//# Charlottesville, VA 22903-2475 USA
//#
-//# $Id$

#ifndef MS_MSMETADATA_H
#define MS_MSMETADATA_H
@@ -32,18 +31,21 @@
#include <casacore/casa/Quanta/QVector.h>
#include <casacore/measures/Measures/MPosition.h>
#include <casacore/ms/MeasurementSets/MeasurementSet.h>
+#include <casacore/ms/MeasurementSets/MSPointingColumns.h>
#include <casacore/casa/Utilities/CountedPtr.h>
#include <map>

namespace casacore {

template <class T> class ArrayColumn;
-
+class ArrayKey;
+class ScanKey;
+class SubScanKey;

// <summary>
// Class to interrogate an MS for metadata. Interrogation happens on
demand
-// and resulting metadata are stored if the cache has not exceeded the
specified
-// limit.
+// and resulting metadata are stored for use by subsequent queries if the
+// cache has not exceeded the specified limit.
// </summary>

class MSMetaData {
@@ -69,31 +71,59 @@
// or else subsequent method calls on this object will result in a
segmentation
// fault; the pointer is not copied.
// <src>maxCacheSizeMB</src> is the maximum cache size in megabytes. <=0
means
- // do not use a cache.
+ // do not use a cache, in which case, each method call will have to
(re)query
+ // the MS. It is highly recommended to use a cache of reasonable size for
the
+ // specified MS if multiple methods are going to be called.
MSMetaData(const MeasurementSet *const &ms, const Float maxCacheSizeMB);

virtual ~MSMetaData();

- // number of unique states (number of rows from the STATE table)
- uInt nStates() const;
+ // get the antenna diameters
+ Quantum<Vector<Double> > getAntennaDiameters();

- // get unique scan numbers
- std::set<Int> getScanNumbers() const;
+ // get the antenna ID for the antenna with the specified name.
+ vector<uInt> getAntennaIDs( const vector<String>& antennaNames);

- std::set<Int> getScansForState(const Int stateID);
+ // get the name of the antenna for the specified antenna ID
+ vector<String> getAntennaNames(
+ std::map<String, uInt>& namesToIDsMap,
+ const vector<uInt>& antennaIDs=vector<uInt>(0)
+ );

- std::set<String> getIntentsForScan(const Int scan) const;
+ // get the antenna stations for the specified antenna IDs
+ vector<String> getAntennaStations(const vector<uInt>& antennaIDs);
+
+ // get the antenna stations for the specified antenna names
+ vector<String> getAntennaStations(const vector<String>& antennaNames);
+
+ std::map<String, std::set<Int> > getIntentToFieldsMap();
+
+ std::map<String, std::set<ScanKey> > getIntentToScansMap();
+
+ std::map<String, std::set<uInt> > getIntentToSpwsMap();
+
+ std::set<String> getIntentsForScan(const ScanKey& scan) const;

// get all intents, in no particular (nor guaranteed) order.
std::set<String> getIntents() const;

+ // get a set of intents corresponding to a specified field
+ std::set<String> getIntentsForField(Int fieldID);
+
+ // get a set of intents corresponding to the specified spectral window
+ std::set<String> getIntentsForSpw(const uInt spw);
+ // get unique scan numbers
+
+ std::set<Int> getScanNumbers(Int obsID, Int arrayID) const;
+
+ std::set<Int> getScansForState(
+ Int stateID, Int obsID, Int arrayID
+ );
+
// get a set of spectral windows for which the specified <src>intent</src>
// applies.
virtual std::set<uInt> getSpwsForIntent(const String& intent);

- // get number of spectral windows
- uInt nSpw(Bool includewvr) const;
-
// get the number of visibilities
uInt nRows() const;

@@ -104,22 +134,16 @@
Int scanNumber, Int fieldID
) const;

- uInt nRows(CorrelationType cType, Int fieldID) const;
+ uInt nRows(CorrelationType cType, uInt fieldID) const;

- // get a set of intents corresponding to the specified spectral window
- std::set<String> getIntentsForSpw(const uInt spw);
+ // get number of spectral windows
+ uInt nSpw(Bool includewvr) const;

- // get a set of intents corresponding to a specified field
- std::set<String> getIntentsForField(Int fieldID);
-
- std::map<String, std::set<Int> > getIntentToFieldsMap();
-
- std::map<String, std::set<Int> > getIntentToScansMap();
-
- std::map<String, std::set<uInt> > getIntentToSpwsMap();
+ // number of unique states (number of rows from the STATE table)
+ uInt nStates() const;

// get the number of fields.
- uInt nFields();
+ uInt nFields() const;

// get a set of spectral windows corresponding to the specified fieldID
std::set<uInt> getSpwsForField(const Int fieldID);
@@ -146,31 +170,15 @@
uInt nDataDescriptions() const;

// get the set of spectral windows for the specified scan number.
- std::set<uInt> getSpwsForScan(const Int scan) const;
+ std::set<uInt> getSpwsForScan(const ScanKey& scan) const;

// get the set of scan numbers for the specified spectral window.
- std::set<Int> getScansForSpw(const uInt spw);
+ std::set<Int> getScansForSpw(uInt spw, Int obsID, Int arrayID) const;

// get the number of antennas in the ANTENNA table
uInt nAntennas() const;

- // get the name of the antenna for the specified antenna ID
- vector<String> getAntennaNames(
- std::map<String, uInt>& namesToIDsMap,
- const vector<uInt>& antennaIDs=vector<uInt>(0)
- );

- // get the antenna ID for the antenna with the specified name.
- vector<uInt> getAntennaIDs( const vector<String>& antennaNames);
-
- // get the antenna stations for the specified antenna IDs
- vector<String> getAntennaStations(const vector<uInt>& antennaIDs);
-
- // get the antenna stations for the specified antenna names
- vector<String> getAntennaStations(const vector<String>& antennaNames);
-
- // get the antenna diameters
- Quantum<Vector<Double> > getAntennaDiameters();

// ALMA-specific. get set of spectral windows used for TDM. These are
windows that have
// 64, 128, or 256 channels
@@ -190,43 +198,69 @@
std::set<uInt> getSQLDSpw();

// Get the scans which fail into the specified time range (center-tol to
center+tol)
- std::set<Int> getScansForTimes(const Double center, const Double tol);
+ std::set<Int> getScansForTimes(
+ Double center, Double tol,
+ Int obsID, Int arrayID
+ ) const;

// Get the times for the specified scans
- std::set<Double> getTimesForScans(std::set<Int> scans) const;
+ std::set<Double> getTimesForScans(std::set<ScanKey> scans) const;

// get the times for the specified scan.
// The return values come from the TIME column.
- std::set<Double> getTimesForScan(const Int scan);
+ std::set<Double> getTimesForScan(const ScanKey& scan) const;

// get the time range for the specified scan. The pair will contain
// the start and stop time of the scan, determined from
min(TIME(x)-0.5*INTERVAL(x)) and
// max(TIME(x)-0.5*INTERVAL(x))
- std::pair<Double, Double> getTimeRangeForScan(Int scan) const;
+ std::pair<Double, Double> getTimeRangeForScan(const ScanKey& scanKey)
const;

// get the times for the specified scan
// std::set<Double> getTimesForScan(const uInt scan) const;

- // get the stateIDs associated with the specified scan number.
- std::set<Int> getStatesForScan(const Int scan);
+ // get the stateIDs associated with the specified scan.
+ std::set<Int> getStatesForScan(uInt obsID, uInt arrayID, Int scan) const;

- // get the scans associated with the specified intent
- std::set<Int> getScansForIntent(const String& intent);
+ // get the position of the specified antenna relative to the observatory
position.
+ // the three vector returned represents the longitudinal, latitudinal,
and elevation
+ // offsets (elements 0, 1, and 2 respectively). The longitude and
latitude offsets are
+ // measured along the surface of a sphere centered at the earth's center
and whose surface
+ // intersects the position of the observatory.
+ Quantum<Vector<Double> > getAntennaOffset(uInt which);

- // get the scan numbers associated with the specified field ID.
- std::set<Int> getScansForFieldID(const Int fieldID);
+ Quantum<Vector<Double> > getAntennaOffset(const String& name);

- // get the scan numbers associated with the specified field. Subclasses
should not implement or override.
- std::set<Int> getScansForField(const String& field);
+ vector<Quantum<Vector<Double> > > getAntennaOffsets() const;
+ // get the positions of the specified antennas. If <src>which</src> is
empty, return
+ // all antenna positions.
+
+ vector<MPosition> getAntennaPositions(
+ const vector<uInt>& which=std::vector<uInt>(0)
+ ) const;
+
+ // <src>names</src> cannot be empty.
+ vector<MPosition> getAntennaPositions(const vector<String>& names);
+
+ std::map<uInt, Double> getAverageIntervalsForScan(const ScanKey& scan)
const;
+
+ vector<uInt> getBBCNos() const;
+
+ std::map<uInt, std::set<uInt> > getBBCNosToSpwMap(SQLDSwitch sqldSwitch);
+
+ vector<vector<Double> > getEdgeChans();
+
+

// get the field IDs for the specified field name. Case insensitive.
- std::set<Int> getFieldIDsForField(const String& field);
+ std::set<Int> getFieldIDsForField(const String& field) const;

// get field IDs associated with the specified scan number.
- std::set<Int> getFieldsForScan(const Int scan);
+ std::set<Int> getFieldsForScan(const ScanKey& scan) const;

// get the field IDs associated with the specified scans
- std::set<Int> getFieldsForScans(const std::set<Int>& scans);
+ std::set<Int> getFieldsForScans(
+ const std::set<Int>& scans, Int obsID, Int arrayID
+ ) const;

// get the field IDs associated with the specified intent.
std::set<Int> getFieldsForIntent(const String& intent);
@@ -245,12 +279,6 @@
// Get the fields which fail into the specified time range (center-tol to
center+tol)
std::set<Int> getFieldsForTimes(Double center, Double tol);

- // get the times for which the specified field was observed
- std::set<Double> getTimesForField(Int fieldID);
-
- // get the time stamps associated with the specified intent
- std::set<Double> getTimesForIntent(const String& intent) const;
-
// get telescope names in the order they are listed in the OBSERVATION
table. These are
// the telescopes (observatories), not the antenna names.
vector<String> getObservatoryNames();
@@ -258,26 +286,30 @@
// get the position of the specified telescope (observatory).
MPosition getObservatoryPosition(uInt which) const;

- // get the positions of the specified antennas. If <src>which</src> is
empty, return
- // all antenna positions.
- vector<MPosition> getAntennaPositions(
- const vector<uInt>& which=std::vector<uInt>(0)
+ // get the scans associated with the specified intent
+ std::set<Int> getScansForIntent(
+ const String& intent, Int obsID, Int arrayID
) const;

- // <src>names</src> cannot be empty.
- vector<MPosition> getAntennaPositions(const vector<String>& names);
+ // get the scan numbers associated with the specified field ID.
+ std::set<Int> getScansForFieldID(Int fieldID, Int obsID, Int arrayID)
const;
+
+ // get the scan numbers associated with the specified field. Subclasses
should not implement or override.
+ std::set<Int> getScansForField(const String& field, Int obsID, Int
arrayID) const;
+
+ // The first value of the pair is spw, the second is polarization ID.
+ std::map<std::pair<uInt, uInt>, Int> getSpwIDPolIDToDataDescIDMap();

- // get the position of the specified antenna relative to the observatory
position.
- // the three vector returned represents the longitudinal, latitudinal,
and elevation
- // offsets (elements 0, 1, and 2 respectively). The longitude and
latitude offsets are
- // measured along the surface of a sphere centered at the earth's center
and whose surface
- // intersects the position of the observatory.
- Quantum<Vector<Double> > getAntennaOffset(uInt which);
+ vector<String> getSpwNames() const;

- Quantum<Vector<Double> > getAntennaOffset(const String& name);
+ // get a data structure, consumable by users, representing a summary of
the dataset
+ Record getSummary() const;

- vector<Quantum<Vector<Double> > > getAntennaOffsets() const;
+ // get the times for which the specified field was observed
+ std::set<Double> getTimesForField(Int fieldID);

+ // get the time stamps associated with the specified intent
+ std::set<Double> getTimesForIntent(const String& intent) const;
Bool hasBBCNo() const;

//std::map<Double, Double> getExposuresForTimes() const;
@@ -302,8 +334,8 @@
Double nUnflaggedRows(CorrelationType cType) const;

Double nUnflaggedRows(
- CorrelationType cType, Int arrayID, Int observationID,
- Int scanNumber, Int fieldID
+ CorrelationType cType, Int arrayID, uInt observationID,
+ Int scanNumber, uInt fieldID
) const;

Double nUnflaggedRows(CorrelationType cType, Int fieldID) const;
@@ -324,18 +356,6 @@

vector<uInt> nChans() const;

- vector<vector<Double> > getEdgeChans();
-
- vector<uInt> getBBCNos() const;
-
- std::map<uInt, std::set<uInt> > getBBCNosToSpwMap(SQLDSwitch sqldSwitch);
-
- vector<String> getSpwNames() const;
-
- std::map<uInt, Double> getAverageIntervalsForScan(Int scan) const;
-
- // The first value of the pair is spw, the second is polarization ID.
- std::map<std::pair<uInt, uInt>, Int> getSpwIDPolIDToDataDescIDMap();

uInt nPol();

@@ -344,15 +364,16 @@
vector<std::map<Int, Quantity> > getFirstExposureTimeMap();

// get polarization IDs for the specified scan and spwid
- std::set<uInt> getPolarizationIDs(Int scan, uInt spwid);
+ std::set<uInt> getPolarizationIDs(uInt obsID, Int arrayID, Int scan, uInt
spwid) const;

// get unique field IDs that exist in the main table.
- std::set<Int> getUniqueFiedIDs();
+ std::set<Int> getUniqueFiedIDs() const;

// get the pointing directions associated with antenna1 and antenna2 for
// the specified row of the main MS table
std::pair<MDirection, MDirection> getPointingDirection(
- Int& ant1, Int& ant2, Double& time, uInt row
+ Int& ant1, Int& ant2, Double& time, uInt row,
+ Bool interpolate=false, Int initialguess=0
) const;

// get the time range for the entire dataset. min(TIME(x) -
0.5*INTERVAL(x)) to
@@ -361,10 +382,6 @@

private:

- // (array_id, observation_id, scan_number, field_id) -> stuff mappings
- typedef std::map<Int, std::map<Int, std::map<Int, std::map<Int, uInt> > >
> AOSFMapI;
- typedef std::map<Int, std::map<Int, std::map<Int, std::map<Int, Double> >
> > AOSFMapD;
-
struct SpwProperties {
Double bandwidth;
QVD chanfreqs;
@@ -382,6 +399,39 @@
uInt bbcno;
String name;
};
+
+ /*
+ struct ArrayProperties {
+ uInt nrows;
+ std::set<Int> scans;
+ };
+
+ struct ScanProperties {
+ Double beginTime;
+ std::set<Int> ddIDs;
+ Double endTime;
+ std::set<Int> fieldIDs;
+ uInt nrows;
+ std::set<Int> stateIDs;
+ };
+ */
+
+ struct SubScanProperties {
+ std::set<Int> antennas;
+ Double beginTime;
+ std::set<uInt> ddIDs;
+ Double endTime;
+ uInt nrows;
+ std::set<Int> stateIDs;
+ };
+
+ // (array_id, observation_id, scan_number, field_id) -> stuff mappings
+ //typedef std::map<Int, std::map<Int, std::map<Int, std::map<Int, uInt> >
> > AOSFMapI;
+ //typedef std::map<Int, std::map<Int, std::map<Int, std::map<Int, Double>
> > > AOSFMapD;
+ //typedef std::map<Int, std::map<Int, std::map<Int, std::map<Int,
SubScanProperties> > > > AOSFMapSSP;
+
+ //typedef std::map<std::pair<Int, Int>, std::pair<Int, Int> > BigKey;
+

// The general pattern is that a mutable gets set only once, on demand,
when its
// setter is called for the first time. If this pattern is broken,
defective behavior
@@ -392,29 +442,43 @@
const Float _maxCacheMB;
mutable uInt _nStates, _nACRows, _nXCRows, _nSpw, _nFields, _nAntennas,
_nObservations, _nScans, _nArrays, _nrows, _nPol, _nDataDescIDs;
+ mutable std::map<ScanKey, std::set<uInt> > _scanToSpwsMap,
_scanToDDIDsMap;
+ mutable std::map<Int, uInt> _dataDescIDToSpwMap, _dataDescIDToPolIDMap
/*, _scanToNRowsMap */;
+ std::map<Int, std::set<uInt> > _fieldToSpwMap;
+ mutable std::map<ScanKey, std::set<Int> > _scanToStatesMap,
_scanToFieldsMap;
+ mutable std::map<Int, std::set<Int> > _fieldToStatesMap,
_stateToFieldsMap, _sourceToFieldsMap /*,
+ _arrayToScansMap, _scanToAntennasMap */;
+ std::map<std::pair<uInt, uInt>, Int> _spwPolIDToDataDescIDMap;
+ std::map<String, uInt> _antennaNameToIDMap;
+ //mutable std::map<ArrayKey, ArrayProperties> _arrayProperties;
+ //mutable std::map<ScanKey, ScanProperties> _scanProperties;
+ mutable std::map<SubScanKey, SubScanProperties> _subScanProperties;
+
+ mutable std::map<String, std::set<Int> > _intentToFieldIDMap;
+ mutable std::map<String, std::set<ScanKey> > _intentToScansMap;
+ mutable std::map<ScanKey, std::pair<Double, Double> > _scanToTimeRangeMap;
+ mutable std::map<ScanKey, std::map<uInt, Double> > _scanSpwToIntervalMap;
+ mutable std::map<std::pair<ScanKey, uInt>, std::set<uInt> >
_scanSpwToPolIDMap;
mutable std::set<String> _uniqueIntents;
- mutable std::map<Int, std::set<uInt> > _scanToSpwsMap, _scanToDDIDsMap;
- mutable std::set<Int> _uniqueScanNumbers, _uniqueFieldIDs,
_uniqueStateIDs;
+ mutable std::set<Int> /*_uniqueScanNumbers, */ _uniqueFieldIDs,
_uniqueStateIDs;
mutable std::set<uInt> _avgSpw, _tdmSpw, _fdmSpw, _wvrSpw, _sqldSpw;
mutable CountedPtr<Vector<Int> > _antenna1, _antenna2, _scans, _fieldIDs,
_stateIDs, _dataDescIDs, _observationIDs, _arrayIDs;
- mutable CountedPtr<AOSFMapI> _scanToNACRowsMap, _scanToNXCRowsMap;
- mutable CountedPtr<std::map<Int, uInt> > _fieldToNACRowsMap,
_fieldToNXCRowsMap;
- mutable std::map<Int, uInt> _dataDescIDToSpwMap, _dataDescIDToPolIDMap;
- std::map<std::pair<uInt, uInt>, Int> _spwPolIDToDataDescIDMap;
- mutable std::map<Int, std::set<String> > _scanToIntentsMap;
+ //mutable CountedPtr<AOSFMapI> _scanToNACRowsMap, _scanToNXCRowsMap;
+ mutable CountedPtr<std::map<SubScanKey, uInt> > _subScanToNACRowsMap,
_subScanToNXCRowsMap;
+ //mutable CountedPtr<std::map<Int, uInt> > _fieldToNACRowsMap,
_fieldToNXCRowsMap;
+
+ mutable CountedPtr<vector<uInt> > _fieldToNACRowsMap, _fieldToNXCRowsMap;
+ mutable std::map<ScanKey, std::set<String> > _scanToIntentsMap;
mutable vector<std::set<String> > _stateToIntentsMap, _spwToIntentsMap,
_fieldToIntentsMap;
mutable vector<SpwProperties> _spwInfo;
- std::map<Int, std::set<uInt> > _fieldToSpwMap;
- mutable vector<std::set<Int> > _spwToFieldIDsMap, _spwToScansMap,
_ddidToScansMap;
- mutable std::map<Int, std::set<Int> > _scanToStatesMap, _scanToFieldsMap,
_fieldToScansMap,
- _fieldToStatesMap, _stateToFieldsMap, _sourceToFieldsMap;
+ mutable vector<std::set<Int> > _spwToFieldIDsMap, _obsToArraysMap;
+ mutable vector<std::set<ScanKey> > _spwToScansMap, _ddidToScansMap,
_fieldToScansMap;
+
mutable vector<String> _fieldNames, _antennaNames, _observatoryNames,
_stationNames;
- std::map<String, uInt> _antennaNameToIDMap;
mutable CountedPtr<Vector<Double> > _times;
CountedPtr<Quantum<Vector<Double> > > _exposures;
- mutable CountedPtr<std::map<Int, std::set<Double> > > _scanToTimesMap;
- mutable std::map<String, std::set<Int> > _intentToFieldIDMap,
_intentToScansMap;
+ mutable CountedPtr<std::map<ScanKey, std::set<Double> > > _scanToTimesMap;
std::map<String, std::set<uInt> > _intentToSpwsMap;
mutable std::map<String, std::set<Double> > _intentToTimesMap;

@@ -427,16 +491,21 @@
Matrix<Bool> _uniqueBaselines;
Quantity _exposureTime;
mutable Double _nUnflaggedACRows, _nUnflaggedXCRows;
- mutable CountedPtr<std::map<Int, Double> > _unflaggedFieldNACRows,
_unflaggedFieldNXCRows;
- mutable CountedPtr<AOSFMapD> _unflaggedScanNACRows, _unflaggedScanNXCRows;
+ mutable CountedPtr<vector<Double> > _unflaggedFieldNACRows,
_unflaggedFieldNXCRows;
+ mutable CountedPtr<std::map<SubScanKey, Double> >
_unflaggedSubScanNACRows, _unflaggedSubScanNXCRows;
+ // mutable CountedPtr<AOSFMapD> _unflaggedScanNACRows,
_unflaggedScanNXCRows;
const String _taqlTableName;
const vector<const Table*> _taqlTempTable;
mutable CountedPtr<ArrayColumn<Bool> > _flagsColumn;
- mutable std::map<Int, std::pair<Double, Double> > _scanToTimeRangeMap;
- mutable std::map<Int, std::map<uInt, Double> > _scanSpwToIntervalMap;
+
mutable Bool _spwInfoStored;
vector<std::map<Int, Quantity> > _firstExposureTimeMap;
- std::map<std::pair<Int, uInt>, std::set<uInt> > _scanSpwToPolIDMap;
+ mutable vector<Int> _numCorrs;
+
+ mutable std::set<ArrayKey> _arrayKeys;
+ mutable std::set<ScanKey> _scanKeys;
+ mutable std::set<SubScanKey> _subscans;
+ mutable std::map<ScanKey, std::set<SubScanKey> > _scanToSubScans;

// disallow copy constructor and = operator
MSMetaData(const MSMetaData&);
@@ -456,11 +525,19 @@
// set metadata from OBSERVATION table
void _setObservation(const MeasurementSet& ms);

- static void _checkScan(const Int scan, const std::set<Int> allScans);
+ void _checkField(uInt fieldID) const;
+
+ //static void _checkScan(const Int scan, const std::set<Int> allScans);
+
+ void _checkScan(const ScanKey& key) const;
+
+ void _checkScans(const std::set<ScanKey>& scanKeys) const;
+
+ void _checkSubScan(const SubScanKey& key) const;

Bool _hasIntent(const String& intent) const;

- Bool _hasFieldID(Int fieldID);
+ Bool _hasFieldID(Int fieldID) const;

Bool _hasStateID(Int stateID);

@@ -491,78 +568,31 @@

CountedPtr<ArrayColumn<Bool> > _getFlags() const;

- std::map<Int, std::set<Int> > _getScanToStatesMap() const;
+ //std::map<Int, std::set<Int> > _getScanToStatesMap() const;

Bool _cacheUpdated(const Float incrementInBytes) const;

- void _getStateToIntentsMap(
- vector<std::set<String> >& statesToIntentsMap,
- std::set<String>& uniqueIntents
- ) const;
+ static void _checkTolerance(const Double tol);

- vector<SpwProperties> _getSpwInfo(
- std::set<uInt>& avgSpw, std::set<uInt>& tdmSpw,
- std::set<uInt>& fdmSpw, std::set<uInt>& wvrSpw,
- std::set<uInt>& sqldSpw
- ) const;
+ //std::map<Int, std::set<Int> > _getArrayIDToScansMap() const;

- static uInt _sizeof(const std::map<Int, std::set<uInt> >& map);
+ std::map<Int, uInt> _getDataDescIDToSpwMap() const;

- static uInt _sizeof(const std::map<Int, std::set<Int> >& map);
-
- static uInt _sizeof(const vector<std::set<Int> >& v);
-
- template <class T> static uInt _sizeof(const std::map<String, std::set<T>
>& map);
-
- static uInt _sizeof(const vector<std::map<Int, Quantity> >& map);
-
- static uInt _sizeof(const std::map<std::pair<Int, uInt>, std::set<uInt>
>& map);
-
- void _getFieldsAndSpwMaps(
- std::map<Int, std::set<uInt> >& fieldToSpwMap,
- vector<std::set<Int> >& spwToFieldMap
- );
-
- void _getScansAndDDIDMaps(
- std::map<Int, std::set<uInt> >& scanToDDIDMap,
- vector<std::set<Int> >& ddIDToScanMap
- ) const;
-
- void _getScansAndSpwMaps(
- std::map<Int, std::set<uInt> >& scanToSpwMap,
- vector<std::set<Int> >& spwToScanMap
- ) const;
+ std::map<Int, uInt> _getDataDescIDToPolIDMap() const;

void _getFieldsAndIntentsMaps(
vector<std::set<String> >& fieldToIntentsMap,
std::map<String, std::set<Int> >& intentToFieldsMap
);

- static uInt _sizeof(const std::map<Int, std::set<String> >& m);
-
- static uInt _sizeof(const vector<std::set<String> >& m);
-
- static uInt _sizeof(const vector<String>& m);
-
- static uInt _sizeof(const Quantum<Vector<Double> >& m);
-
- static uInt _sizeof(const std::map<Int, std::set<Double> >& m);
-
- static uInt _sizeof(const std::map<Double, std::set<Int> >& m);
-
- void _getScansAndIntentsMaps(
- std::map<Int, std::set<String> >& scanToIntentsMap,
- std::map<String, std::set<Int> >& intentToScansMap
+ void _getFieldsAndScansMaps(
+ vector<std::set<ScanKey> >& fieldToScansMap,
+ std::map<ScanKey, std::set<Int> >& scanToFieldsMap
) const;

- void _getSpwsAndIntentsMaps(
- vector<std::set<String> >& spwToIntentsMap,
- std::map<String, std::set<uInt> >& intentToSpwsMap
- );
-
- void _getFieldsAndScansMaps(
- std::map<Int, std::set<Int> >& fieldToScansMap,
- std::map<Int, std::set<Int> >& scanToFieldsMap
+ void _getFieldsAndSpwMaps(
+ std::map<Int, std::set<uInt> >& fieldToSpwMap,
+ vector<std::set<Int> >& spwToFieldMap
);

void _getFieldsAndStatesMaps(
@@ -574,19 +604,33 @@
CountedPtr<std::map<Int, std::set<Double> > >& fieldToTimesMap,
CountedPtr<std::map<Double, std::set<Int> > >& timesToFieldMap
);
+ vector<String> _getFieldNames() const;

- std::map<Int, uInt> _getDataDescIDToSpwMap() const;
+ std::map<String, std::set<Double> > _getIntentsToTimesMap() const;

- std::map<Int, uInt> _getDataDescIDToPolIDMap();
+ vector<String> _getAntennaNames(
+ std::map<String, uInt>& namesToIDsMap
+ );

- vector<String> _getFieldNames() const;
+ vector<MPosition> _getAntennaPositions() const;

- vector<String> _getStationNames();
+ std::map<Double, Double> _getTimeToTotalBWMap(
+ const Vector<Double>& times, const Vector<Int>& ddIDs
+ );

- std::map<String, std::set<Double> > _getIntentsToTimesMap() const;
+ MDirection _getInterpolatedDirection(
+ const ROMSPointingColumns& pCols, const Int& index,
+ const Double& time
+ ) const;

- CountedPtr<std::map<Int, std::set<Double> > > _getScanToTimesMap() const;
+ // number of correlations from the polarzation table.
+ vector<Int> _getNumCorrs() const;

+ vector<std::set<Int> > _getObservationIDToArrayIDsMap() const;
+
+ vector<MPosition> _getObservatoryPositions();
+
+ /*
void _getRowStats(
uInt& nACRows, uInt& nXCRows,
CountedPtr<AOSFMapI>& scanToNACRowsMap,
@@ -594,48 +638,144 @@
CountedPtr<std::map<Int, uInt> >& fieldToNACRowsMap,
CountedPtr<std::map<Int, uInt> >& fieldToNXCRowsMap
) const;
+ */
+ /*
+ void _getRowStats(
+ uInt& nACRows, uInt& nXCRows,
+ AOSFMapI*& scanToNACRowsMap,
+ AOSFMapI*& scanToNXCRowsMap,
+ std::map<Int, uInt>*& fieldToNACRowsMap,
+ std::map<Int, uInt>*& fieldToNXCRowsMap
+ ) const ;
+ */

- void _getUnflaggedRowStats(
- Double& nACRows, Double& nXCRows,
- CountedPtr<AOSFMapD>& scanToNACRowsMap,
- CountedPtr<AOSFMapD>& scanToNXCRowsMap,
- CountedPtr<std::map<Int, Double> >& fieldToNACRowsMap,
- CountedPtr<std::map<Int, Double> >& fieldToNXCRowsMap
+ void _getRowStats(
+ uInt& nACRows, uInt& nXCRows,
+ std::map<SubScanKey, uInt>*& subScanToNACRowsMap,
+ std::map<SubScanKey, uInt>*& subScanToNXCRowsMap,
+ vector<uInt>*& fieldToNACRowsMap,
+ vector<uInt>*& fieldToNXCRowsMap
+ ) const;
+
+ void _getRowStats(
+ uInt& nACRows, uInt& nXCRows,
+ CountedPtr<std::map<SubScanKey, uInt> >& scanToNACRowsMap,
+ CountedPtr<std::map<SubScanKey, uInt> >& scanToNXCRowsMap,
+ CountedPtr<vector<uInt> >& fieldToNACRowsMap,
+ CountedPtr<vector<uInt> >& fieldToNXCRowsMap
+ ) const;
+
+ /*
+ static std::set<Int> _getScans(
+ Int obsID, Int arrayID, std::set<ScanKey>& scanKeys
+ );
+ */
+
+ // get all ScanKeys in the dataset
+ std::set<ScanKey> _getScanKeys() const;
+
+ // get all ScanKeys in the dataset that have the specified arrayKey
+ std::set<ScanKey> _getScanKeys(const ArrayKey& arrayKey) const;
+
+ // get the scan keys in the specified set that have the associated
arrayKey
+ std::set<ScanKey> _getScanKeys(
+ const std::set<ScanKey>& scanKeys, const ArrayKey& arrayKey
+ ) const;
+
+
+ // get all valid scan numbers associated with the specified arrayKey
+ std::set<Int> _getScanNumbers(const ArrayKey& arrayKey) const;
+
+ // get the scan numbers associated with the scanKeys. The members of
+ // scanKeys do not have to have the same ArrayKey, and no warning is
+ // given if they are not.
+ std::set<Int> _getScanNumbers(const std::set<ScanKey>& scanKeys) const;
+
+
+ void _getScansAndDDIDMaps(
+ std::map<ScanKey, std::set<uInt> >& scanToDDIDMap,
+ vector<std::set<ScanKey> >& ddIDToScanMap
+ ) const;
+
+ void _getScansAndIntentsMaps(
+ std::map<ScanKey, std::set<String> >& scanToIntentsMap,
+ std::map<String, std::set<ScanKey> >& intentToScansMap
+ ) const;
+
+ void _getScansAndSpwMaps(
+ std::map<ScanKey, std::set<uInt> >& scanToSpwMap,
+ vector<std::set<ScanKey> >& spwToScanMap
) const;

- void _getTimesAndInvervals(
- std::map<Int, std::pair<Double, Double> >& scanToTimeRangeMap,
- std::map<Int, std::map<uInt, Double> >& scanSpwToIntervalMap
+ // std::map<Int, std::set<Int> > _getScanToAntennasMap() const;
+
+ // std::map<Int, uInt> _getScansToNRowsMap() const;
+
+ std::map<ScanKey, std::set<Int> > _getScanToStatesMap() const;
+
+ std::map<ScanKey, std::set<SubScanKey> > _getScanToSubScansMap() const;
+
+
+ CountedPtr<std::map<ScanKey, std::set<Double> > > _getScanToTimesMap()
const;
+
+
+ vector<SpwProperties> _getSpwInfo(
+ std::set<uInt>& avgSpw, std::set<uInt>& tdmSpw,
+ std::set<uInt>& fdmSpw, std::set<uInt>& wvrSpw,
+ std::set<uInt>& sqldSpw
) const;
+ void _getSpwsAndIntentsMaps(
+ vector<std::set<String> >& spwToIntentsMap,
+ std::map<String, std::set<uInt> >& intentToSpwsMap
+ );

vector<SpwProperties> _getSpwInfo2(
std::set<uInt>& avgSpw, std::set<uInt>& tdmSpw, std::set<uInt>& fdmSpw,
std::set<uInt>& wvrSpw, std::set<uInt>& sqldSpw
) const;

- static void _checkTolerance(const Double tol);
+ void _getStateToIntentsMap(
+ vector<std::set<String> >& statesToIntentsMap,
+ std::set<String>& uniqueIntents
+ ) const;
+
+ vector<String> _getStationNames();

- vector<MPosition> _getObservatoryPositions();
+ std::map<SubScanKey, SubScanProperties> _getSubScanProperties() const;

- vector<String> _getAntennaNames(
- std::map<String, uInt>& namesToIDsMap
- );
+ std::set<SubScanKey> _getSubScanKeys() const;

- vector<MPosition> _getAntennaPositions() const;
+ // get subscans related to the given scan
+ std::set<SubScanKey> _getSubScanKeys(const ScanKey& scanKey) const;

- static std::map<Int, uInt> _toUIntMap(const Vector<Int>& v);
+ //In scanSpwToIntervalMap, the key is a scan, spw pair.
+ void _getTimesAndInvervals(
+ std::map<ScanKey, std::pair<Double, Double> >& scanToTimeRangeMap,
+ std::map<ScanKey, std::map<uInt, Double> >& scanSpwToIntervalMap
+ ) const;

- std::map<Double, Double> _getTimeToTotalBWMap(
- const Vector<Double>& times, const Vector<Int>& ddIDs
- );
+ void _getUnflaggedRowStats(
+ Double& nACRows, Double& nXCRows,
+ CountedPtr<std::map<SubScanKey, Double> >& subScanToNACRowsMap,
+ CountedPtr<std::map<SubScanKey, Double> >& subScanToNXCRowsMap,
+ CountedPtr<vector<Double> >& fieldToNACRowsMap,
+ CountedPtr<vector<Double> >& fieldToNXCRowsMap
+ ) const;

- void _getRowStats(
- uInt& nACRows, uInt& nXCRows,
- AOSFMapI*& scanToNACRowsMap,
- AOSFMapI*& scanToNXCRowsMap,
- std::map<Int, uInt>*& fieldToNACRowsMap,
- std::map<Int, uInt>*& fieldToNXCRowsMap
- )const ;
+ void _getUnflaggedRowStats(
+ Double& nACRows, Double& nXCRows,
+ vector<Double>*& fieldNACRows, vector<Double>*& fieldNXCRows,
+ std::map<SubScanKey, Double>*& scanNACRows,
+ std::map<SubScanKey, Double>*& scanNXCRows
+ ) const;
+ /*
+ void _getUnflaggedRowStats(
+ Double& nACRows, Double& nXCRows,
+ CountedPtr<AOSFMapD>& scanToNACRowsMap,
+ CountedPtr<AOSFMapD>& scanToNXCRowsMap,
+ CountedPtr<std::map<Int, Double> >& fieldToNACRowsMap,
+ CountedPtr<std::map<Int, Double> >& fieldToNXCRowsMap
+ ) const;

void _getUnflaggedRowStats(
Double& nACRows, Double& nXCRows,
@@ -643,8 +783,56 @@
AOSFMapD*& scanNACRows,
AOSFMapD*& scanNXCRows
) const;
+ */
+ template <class T>
+ static uInt _sizeof(const std::map<T, std::set<String> >& m);
+
+ /*
+ static uInt _sizeof(const std::map<Int, std::set<String> >& m);
+
+ static uInt _sizeof(const std::map<ScanKey, std::set<String> >& m);
+ */
+
+ template <class T, class U>
+ static uInt _sizeof(const std::map<T, std::set<U> >& m);
+
+ template <class T, class U>
+ static uInt _sizeof(const std::map<T, U>& m);
+
+ static uInt _sizeof(const vector<std::set<String> >& m);
+
+ static uInt _sizeof(const vector<String>& m);
+
+ template <class T>
+ static uInt _sizeof(const vector<T>& v);
+
+ static uInt _sizeof(const Quantum<Vector<Double> >& m);
+
+ /*
+ static uInt _sizeof(const std::map<Int, std::set<Double> >& m);
+
+ static uInt _sizeof(const std::map<Double, std::set<Int> >& m);
+
+ static uInt _sizeof(const std::map<Int, std::set<uInt> >& map);
+
+ static uInt _sizeof(const std::map<Int, std::set<Int> >& map);
+ */
+
+ template <class T>
+ static uInt _sizeof(const vector<std::set<T> >& v);
+
+ template <class T> static uInt _sizeof(const std::map<String, std::set<T>
>& map);
+
+ static uInt _sizeof(const vector<std::map<Int, Quantity> >& map);
+
+ static uInt _sizeof(const std::map<std::pair<Int, uInt>, std::set<uInt>
>& map);
+
+ static std::map<Int, uInt> _toUIntMap(const Vector<Int>& v);

};
+
+
+
}

#endif
=======================================
--- /branches/nov14/ms/MeasurementSets/MSSummary.cc Tue Feb 10 08:33:29
2015 UTC
+++ /branches/nov14/ms/MeasurementSets/MSSummary.cc Tue Mar 3 12:22:11
2015 UTC
@@ -25,6 +25,7 @@
//#
//# $Id$
//#
+
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays.h>
#include <casacore/casa/Arrays/ArrayMath.h>
@@ -50,6 +51,7 @@
#include <casacore/ms/MeasurementSets.h>
#include <casacore/ms/MeasurementSets/MeasurementSet.h>
#include <casacore/ms/MeasurementSets/MSColumns.h>
+#include <casacore/ms/MeasurementSets/MSKeys.h>
#include <casacore/ms/MeasurementSets/MSMetaData.h>
#include <casacore/ms/MeasurementSets/MSSummary.h>
#include <casacore/ms/MeasurementSets/MSRange.h>
@@ -237,8 +239,6 @@
os << "The MAIN table is empty: there are no data!!!" << endl <<
LogIO::POST;
return;
}
-
- // Make objects
std::pair<Double, Double> timerange = _msmd->getTimeRange();
Double startTime = timerange.first;
Double stopTime = timerange.second;
@@ -351,7 +351,6 @@
}
os << "SpwIds Average Interval(s) ScanIntent"
<< endl;
}
-
/* CAS-2751. Sort the table by scan then field (not timestamp)
* so that scans are not listed for every different DDID
*/
@@ -379,13 +378,12 @@
Int thisnrow(0);
Double meanIntTim(0.0);

-
os.output().precision(3);
Int subsetscan=0;
// Iterate over scans/fields:
std::set<uInt> spw;
+ ScanKey lastScanKey;
while (!stiter.pastEnd()) {
-
// ms table at this scan
Table t(stiter.table());
Int nrow=t.nrow();
@@ -398,6 +396,10 @@
ROTableVector<Int> ddicol(t,"DATA_DESC_ID");
ROTableVector<Int> stidcol(t,"STATE_ID");

+ lastScanKey.obsID = obsid;
+ lastScanKey.arrayID = arrid;
+ lastScanKey.scan = lastscan;
+
// this timestamp
//Double thistime(timecol(0));

@@ -408,7 +410,6 @@
fldids.resize(1,False);
fldids(0)=fldcol(0);
nfld=1;
-
ddids.resize(1,False);
ddids(0)=ddicol(0);
nddi=1;
@@ -416,7 +417,6 @@
stids.resize(1, False);
stids(0) = stidcol(0);
nst=1;
-
nVisPerField_(fldids(0))+=nrow;

// fill field and ddi lists for this scan
@@ -442,11 +442,12 @@
}

// If not first timestamp, check if scan changed, etc.
+ ScanKey thisScanKey = lastScanKey;
+ thisScanKey.scan = thisscan;
if (!firsttime) {
// Has state changed?
Bool samefld;
samefld=fldids.conform(lastfldids) && !anyNE(fldids,lastfldids);
-
Bool sameddi;
sameddi=ddids.conform(lastddids) && !anyNE(ddids,lastddids);

@@ -473,13 +474,10 @@
else {
meanIntTim=0.0;
}
-
// this MJD
day=floor(MVTime(btime/C::day).day());
-
- spw = _msmd->getSpwsForScan(lastscan);
+ spw = _msmd->getSpwsForScan(lastScanKey);
String name=fieldnames(lastfldids(0));
-
if (verbose) {
// Print out last scan's times, fields, ddis
os.output().setf(ios::right, ios::adjustfield);
@@ -502,7 +500,6 @@
os.output().setf(ios::left, ios::adjustfield);
if (name.length()>20) name.replace(19,1,'*');
os.output().width(widthField); os << name.at(0,20);
-
//os.output().width(widthnrow); os << thisnrow;
os.output().width(widthnrow);
os.output().setf(ios::right, ios::adjustfield);
@@ -516,7 +513,6 @@

os.output().width(widthNUnflaggedRow);
os << xx.str();
}
-
//os.output().width(widthInttim); os << meanIntTim;
os.output().width(widthLead); os << " ";

@@ -528,7 +524,7 @@
*/
os.output().width(widthLead); os << " ";

- std::map<uInt, Double> intToScanMap =
_msmd->getAverageIntervalsForScan(lastscan);
+ std::map<uInt, Double> intToScanMap =
_msmd->getAverageIntervalsForScan(lastScanKey);
os << "[";
for (
std::set<uInt>::const_iterator iter=spw.begin();
@@ -548,7 +544,7 @@
//os << obsMode;

//os << endl;
- std::set<String> intents = _msmd->getIntentsForScan(lastscan);
+ std::set<String> intents = _msmd->getIntentsForScan(lastScanKey);
if (intents.size() > 0) {
os << intents;
}
@@ -583,9 +579,7 @@
subsetscan=0;
}
}
-
- // new btime:
- btime = _msmd->getTimeRangeForScan(thisscan).first;
+ btime = _msmd->getTimeRangeForScan(thisScanKey).first;
// next last day is this day
lastday=day;

@@ -593,12 +587,12 @@
meanIntTim=0.;
++recLength;
}
-
// etime=thistime;
//etime=timecol(nrow-1); //CAS-2751
- etime = _msmd->getTimeRangeForScan(thisscan).second;
- } else {
- std::pair<Double, Double> timeRange =
_msmd->getTimeRangeForScan(thisscan);
+ etime = _msmd->getTimeRangeForScan(thisScanKey).second;
+ }
+ else {
+ std::pair<Double, Double> timeRange =
_msmd->getTimeRangeForScan(thisScanKey);
// initialize btime and etime
//btime=thistime;
btime = timeRange.first;
@@ -623,21 +617,19 @@
stiter.next();
} // end of time iteration

- if (thisnrow>0)
+ meanIntTim = thisnrow > 0 ? meanIntTim/thisnrow : 0;
+ /*
+ if
meanIntTim/=thisnrow;
else
meanIntTim=0.0;
-
+ */
// this MJD
day=floor(MVTime(btime/C::day).day());
-
- // Spws
- //spwids.resize(lastddids.nelements());
- /*
- for (uInt iddi=0; iddi<spwids.nelements();++iddi)
- spwids(iddi)=specwindids(lastddids(iddi));
- */
- spw = _msmd->getSpwsForScan(lastscan);
+ lastScanKey.obsID = obsid;
+ lastScanKey.arrayID = arrid;
+ lastScanKey.scan = lastscan;
+ spw = _msmd->getSpwsForScan(lastScanKey);
String name=fieldnames(lastfldids(0));
if (verbose) {
// Print out final scan's times, fields, ddis
@@ -674,14 +666,9 @@
// os.output().width(widthInttim); os << meanIntTim;
os.output().width(widthLead); os << " ";
showContainer (os.output(), spw, ", ", "[", "]");
- /*
- if (spw.size() <= 9) {
- os.output().width(28 - (3*spw.size())); os << " ";
- }
- */
os.output().width(widthLead); os << " ";

- std::map<uInt, Double> intToScanMap =
_msmd->getAverageIntervalsForScan(lastscan);
+ std::map<uInt, Double> intToScanMap =
_msmd->getAverageIntervalsForScan(lastScanKey);
os << "[";
for (
std::set<uInt>::const_iterator iter=spw.begin();
@@ -700,7 +687,7 @@
//}
//os << obsMode;
//os << endl;
- set<String> intents = _msmd->getIntentsForScan(lastscan);
+ set<String> intents = _msmd->getIntentsForScan(lastScanKey);
if (intents.size() > 0) {
os << intents;
}
@@ -951,6 +938,7 @@
subScanRecord.define("nRow", thisnrow);
subScanRecord.define("IntegrationTime", meanIntTim);
subScanRecord.define("SpwIds", spwids);
+ subScanRecord.define("DDIds", lastddids);
scanRecord.defineRecord(String::toString(subsetscan), subScanRecord);
if(!outRec.isDefined(scanrecid)){
outRec.defineRecord(scanrecid, scanRecord);
@@ -1022,6 +1010,7 @@
subScanRecord.define("nRow", thisnrow);
subScanRecord.define("IntegrationTime", meanIntTim);
subScanRecord.define("SpwIds", spwids);
+ subScanRecord.define("DDIds", lastddids);
scanRecord.defineRecord(String::toString(subsetscan), subScanRecord);
if(!outRec.isDefined(scanrecid)){
outRec.defineRecord(scanrecid, scanRecord);
@@ -1672,6 +1661,8 @@
ROMSSpWindowColumns msSWC(pMS->spectralWindow());
// Create a MS-data_desc-columns object
ROMSDataDescColumns msDDC(pMS->dataDescription());
+ // Create a MS-polin-columns object
+ ROMSPolarizationColumns msPOLC(pMS->polarization());

if (msDDC.nrow()<=0 or msSWC.nrow()<=0) {
//The DATA_DESCRIPTION or SPECTRAL_WINDOW table is empty
@@ -1687,6 +1678,7 @@
for (uInt i=0; i<ddId.nelements(); i++) uddId(i)=ddId(i);
// now get the corresponding spectral windows and pol setups
Vector<Int> spwIds = msDDC.spectralWindowId().getColumnCells(uddId);
+ Vector<Int> polIds = msDDC.polarizationId().getColumnCells(uddId);
//const Int option=Sort::HeapSort | Sort::NoDuplicates;
//const Sort::Order order=Sort::Ascending;
//Int nSpw=GenSort<Int>::sort (spwIds, order, option);
@@ -1695,6 +1687,7 @@
// For each row of the DataDesc subtable, write the info
for (uInt i=0; i<ddId.nelements(); i++) {
Int dd=ddId(i);
+ Int pol=polIds(i);
Int spw = msDDC.spectralWindowId()(dd);

Record ddRec;
@@ -1705,6 +1698,8 @@
ddRec.define("ChanWidth", msSWC.chanWidth()(spw)(IPosition(1,0)));
ddRec.define("TotalWidth", msSWC.totalBandwidth()(spw));
ddRec.define("RefFreq", msSWC.refFrequency()(spw));
+ ddRec.define("PolId", polIds(i));
+ ddRec.define("NumCorr", msPOLC.numCorr()(pol));

outRec.defineRecord(String::toString(dd), ddRec);
}
=======================================
--- /branches/nov14/ms/MeasurementSets/test/tMSMetaData.cc Wed Dec 24
11:46:02 2014 UTC
+++ /branches/nov14/ms/MeasurementSets/test/tMSMetaData.cc Tue Mar 3
12:22:11 2015 UTC
@@ -34,6 +34,7 @@
#include <casacore/casa/OS/Directory.h>
#include <casacore/casa/OS/EnvVar.h>
#include <casacore/casa/Quanta/QLogical.h>
+#include <casacore/ms/MeasurementSets/MSKeys.h>
#include <casacore/ms/MeasurementSets/MeasurementSet.h>
#include <casacore/measures/Measures/MDirection.h>

@@ -71,13 +72,16 @@
}

void testIt(MSMetaData& md) {
+ ArrayKey arrayKey;
+ arrayKey.obsID = 0;
+ arrayKey.arrayID = 0;
cout << "*** test nStates()" << endl;
AlwaysAssert(md.nStates() == 43, AipsError);
cout << "*** cache size " << md.getCache() << endl;

cout << "*** test getScansForState()" << endl;
for (uInt stateID=0; stateID<md.nStates(); stateID++) {
- std::set<Int> scans = md.getScansForState(stateID);
+ std::set<Int> scans = md.getScansForState(stateID, 0, 0);
std::set<Int> expec;
if (stateID < 5) {
uInt myints[]= {1, 5, 8};
@@ -120,7 +124,7 @@
cout << "*** cache size " << md.getCache() << endl;

cout << "*** test getScanNumbers()" << endl;
- std::set<Int> scans = md.getScanNumbers();
+ std::set<Int> scans = md.getScanNumbers(0, 0);
uInt myints[] = {
1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32
@@ -133,12 +137,15 @@
}
std::set<String> uniqueIntents;
cout << "*** test getIntentsForScan()" << endl;
-
+ ScanKey scanKey;
+ scanKey.obsID = 0;
+ scanKey.arrayID = 0;
for (
std::set<Int>::const_iterator scanNum = scans.begin();
scanNum!=scans.end(); scanNum++
) {
- std::set<String> intents = md.getIntentsForScan(*scanNum);
+ scanKey.scan = *scanNum;
+ std::set<String> intents = md.getIntentsForScan(scanKey);
std::set<String> exp;

if (*scanNum == 1 || *scanNum == 5 || *scanNum == 8) {
@@ -387,9 +394,13 @@
}
{
cout << "*** test nScans()" << endl;
+ cout << "nscans " << md.nScans() << endl;
AlwaysAssert(md.nScans() == 32, AipsError);
- std::set<Int> scanNumbers = md.getScanNumbers();
+ std::set<Int> scanNumbers = md.getScanNumbers(0, 0);
cout << "*** test getSpwsForScan() and getPolarizationIDs()" << endl;
+ ScanKey scanKey;
+ scanKey.obsID = 0;
+ scanKey.arrayID = 0;
for (
std::set<Int>::const_iterator scan=scanNumbers.begin();
scan!=scanNumbers.end(); scan++
@@ -423,13 +434,14 @@
};
exp.insert(myints, myints+9);
}
- AlwaysAssert(md.getSpwsForScan(*scan) == exp, AipsError);
+ scanKey.scan = *scan;
+ AlwaysAssert(md.getSpwsForScan(scanKey) == exp, AipsError);
for (
std::set<uInt>::const_iterator spw=exp.begin();
spw!=exp.end(); spw++
) {
std::set<uInt> exppols;
- std::set<uInt> pols = md.getPolarizationIDs(*scan, *spw);
+ std::set<uInt> pols = md.getPolarizationIDs(0, 0, *scan, *spw);
if (*spw == 0) {
exppols.insert(1);
}
@@ -472,8 +484,7 @@
else {
// empty set
}
-
- AlwaysAssert(md.getScansForSpw(i) == exp, AipsError);
+ AlwaysAssert(md.getScansForSpw(i, 0, 0) == exp, AipsError);
}
{
cout << "*** test nAntennas()" << endl;
@@ -544,18 +555,17 @@
std::set<Int> exp;
exp.insert(27);
AlwaysAssert(
- md.getScansForTimes(4.84282937e+09,20) == exp,
- AipsError
+ md.getScansForTimes(4.84282937e+09, 20, 0, 0) == exp,
+ AipsError
);
exp.insert(24);
exp.insert(25);
exp.insert(26);
exp.insert(28);
AlwaysAssert(
- md.getScansForTimes(4.84282937e+09,200) == exp,
- AipsError
+ md.getScansForTimes(4.84282937e+09, 200, 0, 0) == exp,
+ AipsError
);
-
}
{
cout << "*** test getTimesForScans()" << endl;
@@ -582,7 +592,7 @@
myscans.insert(3);
myscans.insert(6);
AlwaysAssert(
- allNearAbs(md.getTimesForScans(myscans), expec, 0.1),
+ allNearAbs(md.getTimesForScans(scanKeys(myscans, arrayKey)), expec,
0.1),
AipsError
);
}
@@ -603,14 +613,17 @@
std::set<Int> myscans;
myscans.insert(3);
AlwaysAssert(
- allNearAbs(md.getTimesForScans(myscans), expec, 0.1),
+ allNearAbs(
+ md.getTimesForScans(scanKeys(myscans, arrayKey)),
+ expec, 0.1
+ ),
AipsError
);
}
{
cout << "*** test getStatesForScan()" << endl;
std::set<Int> expec;
- std::set<Int> scanNumbers = md.getScanNumbers();
+ std::set<Int> scanNumbers = md.getScanNumbers(0, 0);
for (
std::set<Int>::const_iterator curScan=scanNumbers.begin();
curScan!=scanNumbers.end(); curScan++
@@ -665,7 +678,7 @@
Int mine[] = {33, 34, 35, 36};
expec.insert(mine, mine+4);
}
- std::set<Int> got = md.getStatesForScan(*curScan);
+ std::set<Int> got = md.getStatesForScan(0, 0, *curScan);
AlwaysAssert(got == expec, AipsError);
}
cout << "*** cache size " << md.getCache() << endl;
@@ -731,9 +744,11 @@
Int mine[] = {12, 16, 20, 23, 27, 30};
expec.insert(mine, mine+6);
}
- std::set<Int> got = md.getScansForIntent(*intent);
- AlwaysAssert(md.getScansForIntent(*intent) == expec, AipsError);
- AlwaysAssert(md.getIntentToScansMap()[*intent] == expec, AipsError);
+ AlwaysAssert(md.getScansForIntent(*intent, 0, 0) == expec, AipsError);
+ AlwaysAssert(
+ casa::scanNumbers(md.getIntentToScansMap()[*intent]) == expec,
+ AipsError
+ );
}
}
{
@@ -780,7 +795,7 @@
default:
throw AipsError("bad fieldID");
}
- AlwaysAssert(md.getScansForFieldID(i) == expec, AipsError);
+ AlwaysAssert(md.getScansForFieldID(i, 0, 0) == expec, AipsError);
}
}
{
@@ -850,14 +865,13 @@
default:
throw AipsError("bad fieldID");
}
- AlwaysAssert(md.getScansForField(name) == expec, AipsError);
+ AlwaysAssert(md.getScansForField(name, 0, 0) == expec, AipsError);
}
cout << "*** cache size " << md.getCache() << endl;
-
}
{
cout << "*** test getFieldsForScan() and getFieldsForScans()" << endl;
- std::set<Int> scans = md.getScanNumbers();
+ std::set<Int> scans = md.getScanNumbers(0, 0);
std::set<Int> expec2;
std::set<Int> curScanSet;
for (
@@ -901,12 +915,16 @@
expec.insert(5);
expec2.insert(5);
}
+ ScanKey scanKey;
+ scanKey.obsID = 0;
+ scanKey.arrayID = 0;
+ scanKey.scan = *curScan;
AlwaysAssert(
- md.getFieldsForScan(*curScan) == expec,
+ md.getFieldsForScan(scanKey) == expec,
AipsError
);
AlwaysAssert(
- md.getFieldsForScans(curScanSet) == expec2,
+ md.getFieldsForScans(curScanSet, 0, 0) == expec2,
AipsError
);
}
@@ -1424,6 +1442,10 @@
AlwaysAssert(nTimes == exp, AipsError);
intent++;
}
+ {
+ cout << "*** test getSummary()" << endl;
+ //cout << md.getSummary() << endl;
+ }
}
{
cout << "*** cache size " << md.getCache() << endl;
=======================================
--- /branches/nov14/scimath/Mathematics/ClassicalStatistics.h Mon Feb 16
07:05:15 2015 UTC
+++ /branches/nov14/scimath/Mathematics/ClassicalStatistics.h Tue Mar 3
12:22:11 2015 UTC
@@ -177,7 +177,7 @@
virtual void setCalculateAsAdded(Bool c);

// An exception will be thrown if setCalculateAsAdded(True) has been
called.
- void setDataProvider(CountedPtr<StatsDataProvider<AccumType,
InputIterator, MaskIterator> > dataProvider);
+ void setDataProvider(StatsDataProvider<AccumType, InputIterator,
MaskIterator> *dataProvider);

void setStatsToCalculate(std::set<StatisticsData::STATS>& stats);

=======================================
--- /branches/nov14/scimath/Mathematics/ClassicalStatistics.tcc Mon Feb 16
07:05:15 2015 UTC
+++ /branches/nov14/scimath/Mathematics/ClassicalStatistics.tcc Tue Mar 3
12:22:11 2015 UTC
@@ -166,7 +166,11 @@
std::set<uInt64> medianIndices;
quantiles.clear();
CountedPtr<uInt64> mynpts = knownNpts.null() ? new uInt64(getNPts()) :
knownNpts;
- if (_getStatsData().median.null()) {
+ ThrowIf(
+ *mynpts == 0,
+ "No valid data found"
+ );
+ if (_getStatsData().median.null()) {
medianIndices = _medianIndices(mynpts);
}
std::map<Double, uInt64> quantileToIndex =
StatisticsData::indicesFromFractions(
@@ -257,7 +261,8 @@
"Value of all quantiles must be between 0 and 1 (noninclusive)"
);
uInt64 mynpts = knownNpts.null() ? getNPts() : *knownNpts;
- std::map<Double, uInt64> quantileToIndexMap =
StatisticsData::indicesFromFractions(
+ ThrowIf(mynpts == 0, "No valid data found");
+ std::map<Double, uInt64> quantileToIndexMap =
StatisticsData::indicesFromFractions(
mynpts, fractions
);
// This seemingly convoluted way of doing things with maps is necessary
because
@@ -296,7 +301,7 @@
Bool c
) {
ThrowIf (
- ! this->_getDataProvider().null() && c,
+ this->_getDataProvider() && c,
"Logic Error: It is nonsensical to call " + String(__func__) + " method "
"with a True value if one is using a data provider"
);
@@ -310,7 +315,7 @@

template <class AccumType, class InputIterator, class MaskIterator>
void ClassicalStatistics<AccumType, InputIterator,
MaskIterator>::setDataProvider(
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
) {
ThrowIf(
_calculateAsAdded,
@@ -429,7 +434,7 @@
_initIterators();
_getStatsData().masked = False;
_getStatsData().weighted = False;
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
while (True) {
_initLoopVars();
@@ -512,7 +517,7 @@
_updateMaxMin(mymin, mymax, minpos, maxpos, _myStride, _idataset);
}
++_idataset;
- if (! dataProvider.null()) {
+ if (dataProvider) {
++(*dataProvider);
if (dataProvider->atEnd()) {
dataProvider->finalize();
@@ -798,7 +803,7 @@
++iDesc;
}
_initIterators();
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
while (True) {
_initLoopVars();
@@ -867,7 +872,7 @@
binDesc, maxLimit
);
}
- if (! dataProvider.null()) {
+ if (dataProvider) {
++(*dataProvider);
if (dataProvider->atEnd()) {
dataProvider->finalize();
@@ -893,7 +898,7 @@
) {
//cout << __func__ << endl;
_initIterators();
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
while (True) {
_initLoopVars();
@@ -954,7 +959,7 @@
ary, _myData, _myCount, _myStride
);
}
- if (! dataProvider.null()) {
+ if (dataProvider) {
++(*dataProvider);
if (dataProvider->atEnd()) {
dataProvider->finalize();
@@ -1003,7 +1008,7 @@
++iLimits;
}
_initIterators();
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
uInt currentCount = 0;
while (True) {
@@ -1072,7 +1077,7 @@
includeLimits, maxCount
);
}
- if (! dataProvider.null()) {
+ if (dataProvider) {
++(*dataProvider);
if (dataProvider->atEnd()) {
dataProvider->finalize();
@@ -1325,7 +1330,7 @@
) {
//cout << __func__ << endl;
_initIterators();
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
CountedPtr<AccumType> mymax;
CountedPtr<AccumType> mymin;
@@ -1387,7 +1392,7 @@
// with it. No filtering of the data is necessary.
_minMax(mymin, mymax, _myData, _myCount, _myStride);
}
- if (! dataProvider.null()) {
+ if (dataProvider) {
++(*dataProvider);
if (dataProvider->atEnd()) {
dataProvider->finalize();
@@ -1416,7 +1421,7 @@
Int64 ClassicalStatistics<AccumType, InputIterator,
MaskIterator>::_doNpts() {
//cout << __func__ << endl;
_initIterators();
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
uInt64 npts = 0;
while (True) {
@@ -1479,7 +1484,7 @@
npts, _myData, _myCount, _myStride
);
}
- if (! dataProvider.null()) {
+ if (dataProvider) {
++(*dataProvider);
if (dataProvider->atEnd()) {
dataProvider->finalize();
@@ -1832,7 +1837,7 @@
const std::set<uInt64>& indices, Bool persistSortedArray
) {
std::map<uInt64, AccumType> indexToValue;
- if (
+ if (
_valuesFromSortedArray(
indexToValue, knownNpts, indices, maxArraySize, persistSortedArray
)
@@ -1878,10 +1883,10 @@
template <class AccumType, class InputIterator, class MaskIterator>
void ClassicalStatistics<AccumType, InputIterator,
MaskIterator>::_initIterators() {
ThrowIf(
- this->_getData().size() == 0 && this->_getDataProvider().null(),
+ this->_getData().size() == 0 && ! this->_getDataProvider(),
"No data sets have been added"
);
- if (! this->_getDataProvider().null()) {
+ if (this->_getDataProvider()) {
this->_getDataProvider()->reset();
}
else {
@@ -1907,9 +1912,9 @@

template <class AccumType, class InputIterator, class MaskIterator>
void ClassicalStatistics<AccumType, InputIterator,
MaskIterator>::_initLoopVars() {
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
- if (! dataProvider.null()) {
+ if (dataProvider) {
_myData = dataProvider->getData();
_myCount = dataProvider->getCount();
_myStride = dataProvider->getStride();
@@ -1957,7 +1962,7 @@
) {
//cout << __func__ << endl;
_initIterators();
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
Bool limitReached = False;
while (True) {
@@ -2029,7 +2034,7 @@
unsortedAry.clear();
return False;
}
- if (! dataProvider.null()) {
+ if (dataProvider) {
++(*dataProvider);
if (dataProvider->atEnd()) {
dataProvider->finalize();
@@ -2923,12 +2928,12 @@
AccumType mymin, AccumType mymax, Int64 minpos, Int64 maxpos, uInt
dataStride,
const Int64& currentDataset
) {
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
if (maxpos >= 0) {
_getStatsData().maxpos.first = currentDataset;
_getStatsData().maxpos.second = maxpos * dataStride;
- if (! dataProvider.null()) {
+ if (dataProvider) {
dataProvider->updateMaxPos(_getStatsData().maxpos);
}
_getStatsData().max = new AccumType(mymax);
@@ -2936,7 +2941,7 @@
if (minpos >= 0) {
_getStatsData().minpos.first = currentDataset;
_getStatsData().minpos.second = minpos * dataStride;
- if (! dataProvider.null()) {
+ if (dataProvider) {
dataProvider->updateMinPos(_getStatsData().minpos);
}
_getStatsData().min = new AccumType(mymin);
@@ -3074,6 +3079,7 @@
? (uInt64)_getStatsData().npts
: knownNpts.null()
? 0 : *knownNpts;
+ ThrowIf(myNpts == 0, "No valid data found");
if (myArray.empty()) {
if (myNpts > 0) {
// we have already computed npts
@@ -3089,7 +3095,7 @@
}
else {
// we have to calculate the number of good points
- if (this->_getDataProvider().null()) {
+ if (! this->_getDataProvider()) {
// we first get an upper limit by adding up the counts
uInt nr = 0;
const vector<Int64>& counts = this->_getCounts();
@@ -3238,6 +3244,5 @@
}

}
-

#endif
=======================================
--- /branches/nov14/scimath/Mathematics/FitToHalfStatistics.tcc Mon Feb 16
07:05:15 2015 UTC
+++ /branches/nov14/scimath/Mathematics/FitToHalfStatistics.tcc Tue Mar 3
12:22:11 2015 UTC
@@ -474,7 +474,7 @@
AccumType mymin, AccumType mymax, Int64 minpos, Int64 maxpos, uInt
dataStride,
const Int64& currentDataset
) {
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
dataProvider
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *dataProvider
= this->_getDataProvider();
if (maxpos >= 0) {
_realMax = new AccumType(mymax);
@@ -483,7 +483,7 @@
_getStatsData().maxpos.second = maxpos * dataStride;
_getStatsData().minpos.first = -1;
_getStatsData().minpos.second = -1;
- if (! dataProvider.null()) {
+ if (dataProvider) {
dataProvider->updateMaxPos(_getStatsData().maxpos);
}
_getStatsData().max = new AccumType(mymax);
@@ -497,7 +497,7 @@
_getStatsData().minpos.second = minpos * dataStride;
_getStatsData().maxpos.first = -1;
_getStatsData().maxpos.second = -1;
- if (! dataProvider.null()) {
+ if (dataProvider) {
dataProvider->updateMinPos(_getStatsData().minpos);
}
_getStatsData().min = new AccumType(mymin);
=======================================
--- /branches/nov14/scimath/Mathematics/StatisticsAlgorithm.h Mon Feb 16
07:05:15 2015 UTC
+++ /branches/nov14/scimath/Mathematics/StatisticsAlgorithm.h Tue Mar 3
12:22:11 2015 UTC
@@ -283,8 +283,9 @@
// instead of settng and adding data "by hand", set the data provider
that will provide
// all the data sets. Calling this method will clear any other data sets
that have
// previously been set or added.
- virtual void setDataProvider(CountedPtr<StatsDataProvider<AccumType,
InputIterator, MaskIterator> > dataProvider) {
- ThrowIf(dataProvider.null(), "Logic Error: data provider cannot be
NULL");
+ virtual void setDataProvider(StatsDataProvider<AccumType, InputIterator,
MaskIterator> *dataProvider) {
+ ThrowIf(! dataProvider, "Logic Error: data provider cannot be NULL");
+ _clearData();
_dataProvider = dataProvider;
}

@@ -310,7 +311,7 @@

const vector<InputIterator>& _getData() const { return _data; }

- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
_getDataProvider() {
+ StatsDataProvider<AccumType, InputIterator, MaskIterator>*
_getDataProvider() {
return _dataProvider;
}

@@ -371,7 +372,7 @@
std::map<uInt, DataRanges> _dataRanges;
vector<AccumType> _sortedArray;
std::set<StatisticsData::STATS> _statsToCalculate, _unsupportedStats;
- CountedPtr<StatsDataProvider<AccumType, InputIterator, MaskIterator> >
_dataProvider;
+ StatsDataProvider<AccumType, InputIterator, MaskIterator> *_dataProvider;

void _throwIfDataProviderDefined() const;
};
=======================================
--- /branches/nov14/scimath/Mathematics/StatisticsAlgorithm.tcc Mon Feb 16
07:05:15 2015 UTC
+++ /branches/nov14/scimath/Mathematics/StatisticsAlgorithm.tcc Tue Mar 3
12:22:11 2015 UTC
@@ -26,9 +26,9 @@

#ifndef SCIMATH_STATISTICSALGORITHM_TCC
#define SCIMATH_STATISTICSALGORITHM_TCC
-
+
#include <casacore/scimath/Mathematics/StatisticsAlgorithm.h>
-
+
#include <casacore/casa/BasicSL/STLIO.h>

namespace casacore {
@@ -37,7 +37,7 @@
StatisticsAlgorithm<AccumType, InputIterator,
MaskIterator>::StatisticsAlgorithm()
: _data(), _weights(), _masks(), _counts(), _dataStrides(), _maskStrides(),
_isIncludeRanges(), _dataRanges(), _sortedArray(), _statsToCalculate(),
- _unsupportedStats(), _dataProvider() {}
+ _unsupportedStats(), _dataProvider(NULL) {}

template <class AccumType, class InputIterator, class MaskIterator>
StatisticsAlgorithm<AccumType, InputIterator, MaskIterator>&
@@ -370,32 +370,15 @@
return indexToValuesMap;
}

-/*
-template <class AccumType, class InputIterator, class MaskIterator>
-std::map<Double, uInt64> StatisticsAlgorithm<AccumType, InputIterator,
MaskIterator>::_indicesFromQuantiles(
- uInt64 npts, const std::set<Double>& quantiles
-) {
- std::map<Double, uInt64> quantileToIndexMap;
- std::set<Double>::const_iterator qiter = quantiles.begin();
- std::set<Double>::const_iterator qend = quantiles.end();
- while (qiter != qend) {
- quantileToIndexMap[*qiter] = ((uInt64)ceil(*qiter * (Double)npts) - 1);
- ++qiter;
- }
- return quantileToIndexMap;
-}
-*/
-
template <class AccumType, class InputIterator, class MaskIterator>
void StatisticsAlgorithm<AccumType, InputIterator,
MaskIterator>::_throwIfDataProviderDefined() const {
ThrowIf(
- ! _dataProvider.null(),
+ _dataProvider,
"Logic Error: Cannot add data after a data provider has been set. Call
setData() to clear "
"the existing data provider and to add this new data set"
);
}

}
-

#endif