Line data Source code
1 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
2 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
3 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
4 : //#
5 : //# This library is free software; you can redistribute it and/or
6 : //# modify it under the terms of the GNU Lesser General Public
7 : //# License as published by the Free software Foundation; either
8 : //# version 2.1 of the License, or (at your option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful,
11 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
12 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 : //# Lesser General Public License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Lesser General Public
16 : //# License along with this library; if not, write to the Free Software
17 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
18 : //# MA 02111-1307 USA
19 :
20 : #ifndef STATWTTVI_H_
21 : #define STATWTTVI_H_
22 :
23 : #include <msvis/MSVis/TransformingVi2.h>
24 :
25 : #include <casacore/casa/aipstype.h>
26 : #include <casacore/casa/Arrays/ArrayIter.h>
27 : #include <casacore/casa/BasicSL/Complex.h>
28 : #include <casacore/ms/MSSel/MSSelection.h>
29 : #include <casacore/scimath/StatsFramework/StatisticsAlgorithmFactory.h>
30 :
31 : #include <msvis/MSVis/VisBuffer2.h>
32 : #include <msvis/MSVis/VisibilityIterator2.h>
33 : #include <mstransform/TVI/StatWtClassicalDataAggregator.h>
34 : #include <mstransform/TVI/StatWtTypes.h>
35 : #include <mstransform/TVI/UtilsTVI.h>
36 : #include <stdcasa/variant.h>
37 : #include <stdcasa/StdCasa/CasacSupport.h>
38 :
39 : #include <map>
40 : #include <vector>
41 : #include <utility>
42 : #include <set>
43 : #include <memory>
44 :
45 : namespace casa {
46 :
47 : class StatWtVarianceAndWeightCalculator;
48 :
49 : namespace vi {
50 :
51 : class StatWtTVI : public TransformingVi2 {
52 :
53 : public:
54 :
55 : static const casacore::String CHANBIN;
56 :
57 : // The following fields are supported in the input configuration record
58 : // combine String, if contains "corr", data will be aggregated
59 : // across correlations.
60 : // value in CHANBIN: Int or Quantity String, describes channel bin widths
61 : // in which to aggregate data within spectral windows
62 : // (spw boundaries are not crossed). If not supplied,
63 : // data for all channels in each spectral window are
64 : // aggregated.
65 : // minsamp: Int, minimum number of samples required in an
66 : // aggregated set, if less than that, stats are not
67 : // computed and the data in the sample are flagged. If
68 : // not supplied, 2 is used.
69 : // statalg String representing what statistics algorithm to use.
70 : // "cl", "ch", "f", "h".
71 : // maxiter Int max number of iterations for Chauvenet algorithm
72 : // zscore Double zscore for Chauvenet algorithm
73 : // center String center for FitToHalf algorithm, "mean",
74 : // "median", or "zero"
75 : // lside Bool side to use for FitToHalf algorithm, True means
76 : // <= center side.
77 : // fence Double fence value for HingesFences algorithm
78 : // wtrange Zero or two element Array<Double>. Specifies the range
79 : // of "good"
80 : // weight values. Data with weights computed to be outside
81 : // this range will be flagged. Both elements must be
82 : // non-negative. If zero length, all weights are
83 : // acceptable.
84 : // fitspw String. MSSelection string representing channels to
85 : // exclude from weight computation.
86 : // datacolumn String. Data column to use for computing weights.
87 : // Supports 'data' or 'corrected'. Minimum match, case
88 : // insensitive. If not provided. 'corrected' is used.
89 : // slidetimebin Bool. If true, use a sliding window for binning in
90 : // time.
91 : // timebin Double. Width of sliding time window. Not used if
92 : // doslidetime is not supplied or if doslidetime = false;
93 : StatWtTVI(
94 : ViImplementation2* inputVii, const casacore::Record &configuration
95 : );
96 :
97 : virtual ~StatWtTVI();
98 :
99 0 : virtual casacore::String ViiType() const {
100 0 : return casacore::String("StatWt( ") + getVii()->ViiType() + " )";
101 : };
102 :
103 : void initWeightSpectrum (const casacore::Cube<casacore::Float>& wtspec);
104 :
105 : void initSigmaSpectrum (const casacore::Cube<casacore::Float>& sigspec);
106 :
107 : void next();
108 :
109 : void origin();
110 :
111 : virtual void weightSpectrum(casacore::Cube<casacore::Float>& wtsp) const;
112 :
113 : virtual void sigmaSpectrum(casacore::Cube<casacore::Float>& sigmaSp) const;
114 :
115 : virtual void weight(casacore::Matrix<casacore::Float> & wtmat) const;
116 :
117 : virtual void sigma(casacore::Matrix<casacore::Float> & sigmaMat) const;
118 :
119 : virtual void flag(casacore::Cube<casacore::Bool>& flagCube) const;
120 :
121 : virtual void flagRow (casacore::Vector<casacore::Bool> & flagRow) const;
122 :
123 : void summarizeFlagging() const;
124 :
125 : void summarizeStats(
126 : casacore::Double& mean, casacore::Double& variance
127 : ) const;
128 :
129 : // Override unimplemented TransformingVi2 version
130 : void writeBackChanges(VisBuffer2* vb);
131 :
132 : // these are public so that class StatWt can call them. In general, other
133 : // clients shouldn't call them.
134 : static casacore::Double getTimeBinWidthInSec(
135 : const casacore::Quantity& binWidth
136 : );
137 :
138 : static void checkTimeBinWidth(casacore::Double binWidth);
139 :
140 : protected:
141 :
142 : void originChunks(casacore::Bool forceRewind);
143 :
144 : void nextChunk();
145 :
146 : private:
147 :
148 : mutable casacore::Bool _weightsComputed = false;
149 : mutable std::shared_ptr<casacore::Bool> _mustComputeWtSp {};
150 : mutable casacore::Cube<casacore::Float> _newWtSp {};
151 : mutable casacore::Matrix<casacore::Float> _newWt {};
152 : mutable casacore::Cube<casacore::Bool> _newFlag {};
153 : mutable casacore::Vector<casacore::Bool> _newFlagRow {};
154 : // The key refers to the spw, the value vector refers to the
155 : // channel numbers within that spw that are the first, last channel pair
156 : // in their respective bins
157 : std::map<casacore::Int, std::vector<StatWtTypes::ChanBin>> _chanBins {};
158 : casacore::Int _minSamp = 2;
159 : casacore::Bool _combineCorr {false};
160 : std::shared_ptr<
161 : casacore::StatisticsAlgorithm<
162 : casacore::Double, casacore::Array<casacore::Float>::const_iterator,
163 : casacore::Array<casacore::Bool>::const_iterator,
164 : casacore::Array<casacore::Double>::const_iterator
165 : >
166 : > _statAlg {} ;
167 : std::shared_ptr<std::pair<casacore::Double, casacore::Double>> _wtrange {};
168 : // The _chanSelFlags key is the spw. The value is a Cube for convenience
169 : // for subchunk computations that require the same shaped cube of flags to
170 : // be applied. The dimension that counts is the second (zero-based 1) as it
171 : // has length equal to the number of channels in the spw. A value of True
172 : // indicates that the channel is "flagged", ie should not be used.
173 : std::map<casacore::uInt, casacore::Cube<casacore::Bool>> _chanSelFlags {};
174 :
175 : mutable size_t _nTotalPts = 0;
176 : mutable size_t _nNewFlaggedPts = 0;
177 : mutable size_t _nOrigFlaggedPts = 0;
178 : mutable StatWtTypes::Column _column = StatWtTypes::CORRECTED;
179 : mutable std::shared_ptr<
180 : std::map<casacore::uInt, std::pair<casacore::uInt, casacore::uInt>>
181 : > _samples {
182 : new std::map<
183 : casacore::uInt, std::pair<casacore::uInt, casacore::uInt>
184 : >()
185 : };
186 : // mutable std::set<casacore::uInt> _processedRowIDs {};
187 : // mutable std::vector<std::vector<casacore::Double>> _timeWindowWts {};
188 : // mutable casacore::Cube<casacore::Double> _multiLoopWeights {};
189 : // if False, the a sliding time window is being used
190 : casacore::Bool _timeBlockProcessing = true;
191 : // we can process using classical VI/VB2 algorithm. Only happens if
192 : // we are not using a sliding time window and if we are not using an
193 : // integer number of time bins
194 : casacore::Bool _doClassicVIVB = true;
195 : // if defined means we are using a window width in seconds
196 : std::shared_ptr<casacore::Double> _binWidthInSeconds {};
197 : // if defined means we are using an integer number of timestamps for the
198 : // bin width
199 : std::shared_ptr<casacore::Int> _nTimeStampsInBin {};
200 :
201 : casacore::Bool _mustComputeSigma = casacore::False;
202 : casacore::Bool _updateWeight = casacore::True;
203 : casacore::Bool _noModel = casacore::False;
204 :
205 : std::shared_ptr<StatWtDataAggregator> _dataAggregator {};
206 :
207 : std::shared_ptr<
208 : casacore::ClassicalStatistics<casacore::Double,
209 : casacore::Array<casacore::Float>::const_iterator,
210 : casacore::Array<casacore::Bool>::const_iterator>
211 : > _wtStats {};
212 :
213 : void _computeWeightSpectrumAndFlags() const;
214 :
215 : // combines the flag cube with the channel selection flags (if any)
216 : casacore::Cube<casacore::Bool> _getResultantFlags(
217 : casacore::Cube<casacore::Bool>& chanSelFlagTemplate,
218 : casacore::Cube<casacore::Bool>& chanSelFlags,
219 : casacore::Bool& initChanSelFlags, casacore::Int spw,
220 : const casacore::Cube<casacore::Bool>& flagCube
221 : ) const;
222 :
223 : // CAS-12358
224 : void _logUsedChannels() const;
225 :
226 : casacore::Bool _parseConfiguration(const casacore::Record &configuration);
227 :
228 : std::pair<
229 : casacore::Cube<casacore::Float>, casacore::Cube<casacore::Bool>
230 : > _getLowerLayerWtSpFlags(size_t& nOrigFlagged) const;
231 :
232 : void _setChanBinMap(casacore::Int binWidth);
233 :
234 : void _setChanBinMap(const casacore::Quantity& binWidth);
235 :
236 : void _setDefaultChanBinMap();
237 :
238 : void _clearCache();
239 :
240 : void _configureStatAlg(const casacore::Record& config);
241 :
242 : };
243 :
244 : }
245 :
246 : }
247 :
248 : #endif
249 :
|
