Line data Source code
1 : //#
2 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
3 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
4 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
5 : //#
6 : //# This library is free software; you can redistribute it and/or
7 : //# modify it under the terms of the GNU Lesser General Public
8 : //# License as published by the Free software Foundation; either
9 : //# version 2.1 of the License, or (at your option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful,
12 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
13 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 : //# Lesser General Public License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Lesser General Public
17 : //# License along with this library; if not, write to the Free Software
18 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
19 : //# MA 02111-1307 USA
20 :
21 : #include <mstransform/MSTransform/StatWtColConfig.h>
22 :
23 : #include <casacore/tables/DataMan/TiledShapeStMan.h>
24 : #include <casacore/tables/Tables/ArrColDesc.h>
25 : #include <casacore/tables/Tables/ArrayColumn.h>
26 :
27 : #include <msvis/MSVis/IteratingParameters.h>
28 : #include <msvis/MSVis/LayeredVi2Factory.h>
29 : #include <stdcasa/variant.h>
30 :
31 : using namespace casacore;
32 : using namespace casac;
33 :
34 : namespace casa {
35 :
36 48 : StatWtColConfig::StatWtColConfig(
37 : casacore::MeasurementSet* ms, casacore::MeasurementSet* selMS, Bool preview,
38 : const String& dataColumn, const variant& chanbin
39 48 : ) : _ms(ms), _selMS(selMS), _preview(preview), _dataColumn(dataColumn) {
40 48 : ThrowIf(_dataColumn.empty(), "data column cannot be empty");
41 48 : _dataColumn.downcase();
42 48 : ThrowIf (
43 : ! (
44 : _dataColumn.startsWith("c")
45 : || _dataColumn.startsWith("d")
46 : || _dataColumn.startsWith("residual")
47 : || _dataColumn.startsWith("residual_")
48 : ),
49 : "Unsupported value for data column: " + dataColumn
50 : );
51 96 : _possiblyWriteSigma = _dataColumn.startsWith("d")
52 96 : || _dataColumn.startsWith("residual_");
53 48 : auto chanBinType = chanbin.type();
54 48 : ThrowIf(
55 : (chanBinType == variant::BOOLVEC && ! chanbin.toBoolVec().empty())
56 : && chanBinType != variant::STRING && chanBinType != variant::INT,
57 : "Unsupported data type for chanbin"
58 : );
59 48 : _doChanBin = (
60 89 : chanBinType == variant::STRING && chanbin.toString() != "spw"
61 89 : ) || chanBinType == variant::INT;
62 48 : _determineFlags();
63 48 : _initSpecColsIfNecessary();
64 48 : }
65 :
66 48 : StatWtColConfig::~StatWtColConfig() {}
67 :
68 47 : void StatWtColConfig::getColWriteFlags(
69 : casacore::Bool& mustWriteWt, casacore::Bool& mustWriteWtSp,
70 : casacore::Bool& mustWriteSig, casacore::Bool& mustWriteSigSp
71 : ) const {
72 47 : mustWriteWt = _mustWriteWt;
73 47 : mustWriteWtSp = _mustWriteWtSp;
74 47 : mustWriteSig = _mustWriteSig;
75 47 : mustWriteSigSp = _mustWriteSigSp;
76 47 : }
77 :
78 48 : void StatWtColConfig::_determineFlags() {
79 48 : _mustWriteSig = _possiblyWriteSigma && ! _preview;
80 48 : Bool hasSigSp = False;
81 48 : Bool sigSpIsInitialized = False;
82 48 : _hasSpectrumIsSpectrumInitialized(
83 : hasSigSp, sigSpIsInitialized, MS::SIGMA_SPECTRUM
84 : );
85 48 : _mustWriteSigSp = _mustWriteSig && (sigSpIsInitialized || _doChanBin);
86 96 : _mustWriteWt = ! _preview
87 95 : && (
88 47 : ! _mustWriteSig
89 8 : || (
90 8 : _mustWriteSig
91 8 : && ! _ms->isColumn(MSMainEnums::CORRECTED_DATA)
92 : )
93 : );
94 48 : Bool hasWtSp = False;
95 48 : Bool wtSpIsInitialized = False;
96 48 : _hasSpectrumIsSpectrumInitialized(
97 : hasWtSp, wtSpIsInitialized, MS::WEIGHT_SPECTRUM
98 : );
99 48 : _mustWriteWtSp = _mustWriteWt && (wtSpIsInitialized || _doChanBin);
100 48 : static const auto colNameWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
101 : static const auto descWtSp = "weight spectrum";
102 : static const auto mgrWtSp = "TiledWgtSpectrum";
103 96 : _dealWithSpectrumColumn(
104 48 : hasWtSp, _mustWriteWtSp, _mustInitWtSp, _mustWriteWt,
105 : colNameWtSp, descWtSp, wtSpIsInitialized, mgrWtSp
106 : );
107 48 : static const auto colNameSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
108 : static const auto descSigSp = "sigma spectrum";
109 : static const auto mgrSigSp = "TiledSigSpectrum";
110 96 : _dealWithSpectrumColumn(
111 48 : hasSigSp, _mustWriteSigSp, _mustInitSigSp, _mustWriteSig,
112 : colNameSigSp, descSigSp, sigSpIsInitialized, mgrSigSp
113 : );
114 96 : LogIO log(LogOrigin("StatWtColConfig", __func__));
115 48 : if (_mustWriteWt) {
116 45 : if (_mustWriteSig) {
117 : log << LogIO::NORMAL
118 : << "CORRECTED_DATA is not present. Updating the "
119 : << "SIGMA/SIGMA_SPECTRUM and WEIGHT/WEIGHT_SPECTRUM values "
120 : << "based on calculations using the DATA column."
121 6 : << LogIO::POST;
122 : }
123 : else {
124 : log << LogIO::NORMAL
125 : << "Updating the WEIGHT/WEIGHT_SPECTRUM values. SIGMA/SIGMA_SPECTRUM "
126 : << "values will not be recalculated as they are related to the values "
127 39 : << "in the DATA column." << LogIO::POST;
128 : }
129 : }
130 3 : else if (_mustWriteSig) {
131 : log << LogIO::NORMAL
132 : << "Updating the SIGMA/SIGMA_SPECTRUM values. WEIGHT/WEIGHT_SPECTRUM will "
133 : << "not be recalculated as they are related to the values in the "
134 2 : << "CORRECTED_DATA column." << LogIO::POST;
135 : }
136 48 : }
137 :
138 48 : void StatWtColConfig::_initSpecColsIfNecessary() {
139 48 : if (! _mustInitWtSp && ! _mustInitSigSp) {
140 43 : return;
141 : }
142 10 : LogIO log(LogOrigin("StatWtColConfig", __func__));
143 5 : if (_mustInitWtSp) {
144 : log << LogIO::NORMAL
145 : << "Fully initializing WEIGHT_SPECTRUM column"
146 5 : << LogIO::POST;
147 : }
148 5 : if (_mustInitWtSp) {
149 : log << LogIO::NORMAL
150 : << "Fully initializing SIGMA_SPECTRUM column"
151 5 : << LogIO::POST;
152 : }
153 5 : std::vector<Int> scs;
154 5 : scs.push_back(MS::ARRAY_ID);
155 5 : scs.push_back(MS::DATA_DESC_ID);
156 5 : scs.push_back(MS::TIME);
157 5 : Block<int> sort(scs.size());
158 5 : uInt i = 0;
159 20 : for (const auto& col: scs) {
160 15 : sort[i] = col;
161 15 : ++i;
162 : }
163 5 : vi::SortColumns sc(sort, False);
164 5 : vi::IteratingParameters ipar;
165 15 : for (uInt i=0; i<2; ++i) {
166 10 : if (i == 1 && _ms == _selMS) {
167 0 : break;
168 : }
169 : // FIXME, in some cases, this still doesn't add the needed columns
170 : // to _selMS
171 10 : MeasurementSet* ms = i == 0 ? _ms : _selMS;
172 10 : vi::VisIterImpl2LayerFactory mydata(ms, ipar, True);
173 10 : Vector<vi::ViiLayerFactory*> facts(1);
174 10 : facts[0] = &mydata;
175 10 : vi::VisibilityIterator2 vi(facts);
176 10 : vi::VisBuffer2 *vb = vi.getVisBuffer();
177 602 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
178 2100 : for (vi.origin(); vi.more(); vi.next()) {
179 1508 : auto nrow = vb->nRows();
180 1508 : auto nchan = vb->nChannels();
181 1508 : auto ncor = vb->nCorrelations();
182 1508 : if (_mustInitWtSp) {
183 1508 : Cube<Float> newsp(ncor, nchan, nrow, 0);
184 1508 : _setEqual(newsp, vb->weight());
185 1508 : vb->initWeightSpectrum(newsp);
186 1508 : }
187 1508 : if (_mustInitSigSp) {
188 1148 : Cube<Float> newsp(ncor, nchan, nrow, 0);
189 1148 : _setEqual(newsp, vb->sigma());
190 1148 : vb->initSigmaSpectrum(newsp);
191 1148 : }
192 1508 : vb->writeChangesBack();
193 1508 : if (_mustInitWtSp) {
194 1508 : AlwaysAssert(vb->weightSpectrum().size() > 0, AipsError);
195 : }
196 1508 : if (_mustInitSigSp) {
197 1148 : AlwaysAssert(vb->sigmaSpectrum().size() > 0, AipsError);
198 : }
199 : }
200 : }
201 10 : }
202 5 : }
203 :
204 2656 : void StatWtColConfig::_setEqual(
205 : Cube<Float>& newsp, const Matrix<Float>& col
206 : ) {
207 2656 : IPosition start(3, 0);
208 2656 : auto shape = newsp.shape();
209 2656 : auto end = shape - 1;
210 12560 : for (Int corr=0; corr<shape[0]; ++corr) {
211 104320 : for (Int row=0; row<shape[2]; ++row) {
212 94416 : start[0] = corr;
213 94416 : end[0] = corr;
214 94416 : start[2] = row;
215 94416 : end[2] = row;
216 94416 : newsp(start, end) = col(corr, row);
217 : }
218 : }
219 2656 : }
220 :
221 96 : void StatWtColConfig::_hasSpectrumIsSpectrumInitialized(
222 : bool& hasSpectrum, bool& spectrumIsInitialzed,
223 : MS::PredefinedColumns col
224 : ) const {
225 96 : hasSpectrum = _ms->isColumn(col);
226 96 : if (! hasSpectrum) {
227 : // no column, so it is obviously not initialized
228 11 : spectrumIsInitialzed = False;
229 11 : return;
230 : }
231 85 : ArrayColumn<Float> column(*_ms, MS::columnName(col));
232 : try {
233 85 : column.get(0);
234 : // we were able to get a row, so its initialized.
235 85 : spectrumIsInitialzed = True;
236 : }
237 0 : catch (const AipsError& x) {
238 : // attempt to get first row failed, its not initialized.
239 0 : spectrumIsInitialzed = False;
240 0 : }
241 85 : }
242 :
243 96 : void StatWtColConfig::_dealWithSpectrumColumn(
244 : Bool& hasSpec, Bool& mustWriteSpec, Bool& mustInitSpec,
245 : Bool mustWriteNonSpec, const String& colName, const String& descName,
246 : Bool specIsInitialized, const String& mgrName
247 : ) {
248 : // this conditional structure supports the
249 : // case of ! hasSpec && ! mustWriteSpec, in which case,
250 : // nothing need be done
251 96 : if (! hasSpec) {
252 11 : if (mustWriteSpec) {
253 : // we must create spectrum column
254 7 : hasSpec = True;
255 7 : mustInitSpec = True;
256 : // from Calibrater.cc
257 : // Nominal default tile shape
258 7 : IPosition dts(3, 4, 32, 1024);
259 : // Discern DATA's default tile shape and use it
260 7 : const auto dminfo = _ms->dataManagerInfo();
261 119 : for (uInt i=0; i<dminfo.nfields(); ++i) {
262 119 : Record col = dminfo.asRecord(i);
263 119 : if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
264 14 : dts = IPosition(
265 14 : col.asRecord("SPEC").asArrayInt("DEFAULTTILESHAPE")
266 7 : );
267 7 : break;
268 : }
269 119 : }
270 : // Add the column
271 7 : String colWtSp = colName;
272 7 : TableDesc tdWtSp;
273 7 : tdWtSp.addColumn(ArrayColumnDesc<Float>(colWtSp, descName, 2));
274 7 : TiledShapeStMan wtSpStMan(mgrName, dts);
275 7 : _ms->addColumn(tdWtSp, wtSpStMan);
276 7 : }
277 : }
278 85 : else if (mustWriteNonSpec) {
279 42 : if (specIsInitialized) {
280 : // it's initialized, so even if we are using the full
281 : // spw for binning, we still need to update *_SPECTRUM
282 42 : mustWriteSpec = True;
283 : }
284 : else {
285 : // it's not initialized, so we aren't going to write to it unless
286 : // chanbin has been specified to be less than the spw width
287 0 : mustInitSpec = mustWriteSpec;
288 : }
289 : }
290 96 : }
291 :
292 : }
|
