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 0 : StatWtColConfig::StatWtColConfig(
37 : casacore::MeasurementSet* ms, casacore::MeasurementSet* selMS, Bool preview,
38 : const String& dataColumn, const variant& chanbin
39 0 : ) : _ms(ms), _selMS(selMS), _preview(preview), _dataColumn(dataColumn) {
40 0 : ThrowIf(_dataColumn.empty(), "data column cannot be empty");
41 0 : _dataColumn.downcase();
42 0 : 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 0 : _possiblyWriteSigma = _dataColumn.startsWith("d")
52 0 : || _dataColumn.startsWith("residual_");
53 0 : auto chanBinType = chanbin.type();
54 0 : ThrowIf(
55 : (chanBinType == variant::BOOLVEC && ! chanbin.toBoolVec().empty())
56 : && chanBinType != variant::STRING && chanBinType != variant::INT,
57 : "Unsupported data type for chanbin"
58 : );
59 0 : _doChanBin = (
60 0 : chanBinType == variant::STRING && chanbin.toString() != "spw"
61 0 : ) || chanBinType == variant::INT;
62 0 : _determineFlags();
63 0 : _initSpecColsIfNecessary();
64 0 : }
65 :
66 0 : StatWtColConfig::~StatWtColConfig() {}
67 :
68 0 : void StatWtColConfig::getColWriteFlags(
69 : casacore::Bool& mustWriteWt, casacore::Bool& mustWriteWtSp,
70 : casacore::Bool& mustWriteSig, casacore::Bool& mustWriteSigSp
71 : ) const {
72 0 : mustWriteWt = _mustWriteWt;
73 0 : mustWriteWtSp = _mustWriteWtSp;
74 0 : mustWriteSig = _mustWriteSig;
75 0 : mustWriteSigSp = _mustWriteSigSp;
76 0 : }
77 :
78 0 : void StatWtColConfig::_determineFlags() {
79 0 : _mustWriteSig = _possiblyWriteSigma && ! _preview;
80 0 : Bool hasSigSp = False;
81 0 : Bool sigSpIsInitialized = False;
82 0 : _hasSpectrumIsSpectrumInitialized(
83 : hasSigSp, sigSpIsInitialized, MS::SIGMA_SPECTRUM
84 : );
85 0 : _mustWriteSigSp = _mustWriteSig && (sigSpIsInitialized || _doChanBin);
86 0 : _mustWriteWt = ! _preview
87 0 : && (
88 0 : ! _mustWriteSig
89 0 : || (
90 0 : _mustWriteSig
91 0 : && ! _ms->isColumn(MSMainEnums::CORRECTED_DATA)
92 : )
93 : );
94 0 : Bool hasWtSp = False;
95 0 : Bool wtSpIsInitialized = False;
96 0 : _hasSpectrumIsSpectrumInitialized(
97 : hasWtSp, wtSpIsInitialized, MS::WEIGHT_SPECTRUM
98 : );
99 0 : _mustWriteWtSp = _mustWriteWt && (wtSpIsInitialized || _doChanBin);
100 0 : static const auto colNameWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
101 : static const auto descWtSp = "weight spectrum";
102 : static const auto mgrWtSp = "TiledWgtSpectrum";
103 0 : _dealWithSpectrumColumn(
104 0 : hasWtSp, _mustWriteWtSp, _mustInitWtSp, _mustWriteWt,
105 : colNameWtSp, descWtSp, wtSpIsInitialized, mgrWtSp
106 : );
107 0 : static const auto colNameSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
108 : static const auto descSigSp = "sigma spectrum";
109 : static const auto mgrSigSp = "TiledSigSpectrum";
110 0 : _dealWithSpectrumColumn(
111 0 : hasSigSp, _mustWriteSigSp, _mustInitSigSp, _mustWriteSig,
112 : colNameSigSp, descSigSp, sigSpIsInitialized, mgrSigSp
113 : );
114 0 : LogIO log(LogOrigin("StatWtColConfig", __func__));
115 0 : if (_mustWriteWt) {
116 0 : 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 0 : << 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 0 : << "in the DATA column." << LogIO::POST;
128 : }
129 : }
130 0 : 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 0 : << "CORRECTED_DATA column." << LogIO::POST;
135 : }
136 0 : }
137 :
138 0 : void StatWtColConfig::_initSpecColsIfNecessary() {
139 0 : if (! _mustInitWtSp && ! _mustInitSigSp) {
140 0 : return;
141 : }
142 0 : LogIO log(LogOrigin("StatWtColConfig", __func__));
143 0 : if (_mustInitWtSp) {
144 : log << LogIO::NORMAL
145 : << "Fully initializing WEIGHT_SPECTRUM column"
146 0 : << LogIO::POST;
147 : }
148 0 : if (_mustInitWtSp) {
149 : log << LogIO::NORMAL
150 : << "Fully initializing SIGMA_SPECTRUM column"
151 0 : << LogIO::POST;
152 : }
153 0 : std::vector<Int> scs;
154 0 : scs.push_back(MS::ARRAY_ID);
155 0 : scs.push_back(MS::DATA_DESC_ID);
156 0 : scs.push_back(MS::TIME);
157 0 : Block<int> sort(scs.size());
158 0 : uInt i = 0;
159 0 : for (const auto& col: scs) {
160 0 : sort[i] = col;
161 0 : ++i;
162 : }
163 0 : vi::SortColumns sc(sort, False);
164 0 : vi::IteratingParameters ipar;
165 0 : for (uInt i=0; i<2; ++i) {
166 0 : 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 0 : MeasurementSet* ms = i == 0 ? _ms : _selMS;
172 0 : vi::VisIterImpl2LayerFactory mydata(ms, ipar, True);
173 0 : Vector<vi::ViiLayerFactory*> facts(1);
174 0 : facts[0] = &mydata;
175 0 : vi::VisibilityIterator2 vi(facts);
176 0 : vi::VisBuffer2 *vb = vi.getVisBuffer();
177 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
178 0 : for (vi.origin(); vi.more(); vi.next()) {
179 0 : auto nrow = vb->nRows();
180 0 : auto nchan = vb->nChannels();
181 0 : auto ncor = vb->nCorrelations();
182 0 : if (_mustInitWtSp) {
183 0 : Cube<Float> newsp(ncor, nchan, nrow, 0);
184 0 : _setEqual(newsp, vb->weight());
185 0 : vb->initWeightSpectrum(newsp);
186 : }
187 0 : if (_mustInitSigSp) {
188 0 : Cube<Float> newsp(ncor, nchan, nrow, 0);
189 0 : _setEqual(newsp, vb->sigma());
190 0 : vb->initSigmaSpectrum(newsp);
191 : }
192 0 : vb->writeChangesBack();
193 0 : if (_mustInitWtSp) {
194 0 : AlwaysAssert(vb->weightSpectrum().size() > 0, AipsError);
195 : }
196 0 : if (_mustInitSigSp) {
197 0 : AlwaysAssert(vb->sigmaSpectrum().size() > 0, AipsError);
198 : }
199 : }
200 : }
201 : }
202 : }
203 :
204 0 : void StatWtColConfig::_setEqual(
205 : Cube<Float>& newsp, const Matrix<Float>& col
206 : ) {
207 0 : IPosition start(3, 0);
208 0 : auto shape = newsp.shape();
209 0 : auto end = shape - 1;
210 0 : for (Int corr=0; corr<shape[0]; ++corr) {
211 0 : for (Int row=0; row<shape[2]; ++row) {
212 0 : start[0] = corr;
213 0 : end[0] = corr;
214 0 : start[2] = row;
215 0 : end[2] = row;
216 0 : newsp(start, end) = col(corr, row);
217 : }
218 : }
219 0 : }
220 :
221 0 : void StatWtColConfig::_hasSpectrumIsSpectrumInitialized(
222 : bool& hasSpectrum, bool& spectrumIsInitialzed,
223 : MS::PredefinedColumns col
224 : ) const {
225 0 : hasSpectrum = _ms->isColumn(col);
226 0 : if (! hasSpectrum) {
227 : // no column, so it is obviously not initialized
228 0 : spectrumIsInitialzed = False;
229 0 : return;
230 : }
231 0 : ArrayColumn<Float> column(*_ms, MS::columnName(col));
232 : try {
233 0 : column.get(0);
234 : // we were able to get a row, so its initialized.
235 0 : spectrumIsInitialzed = True;
236 : }
237 0 : catch (const AipsError& x) {
238 : // attempt to get first row failed, its not initialized.
239 0 : spectrumIsInitialzed = False;
240 : }
241 : }
242 :
243 0 : 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 0 : if (! hasSpec) {
252 0 : if (mustWriteSpec) {
253 : // we must create spectrum column
254 0 : hasSpec = True;
255 0 : mustInitSpec = True;
256 : // from Calibrater.cc
257 : // Nominal default tile shape
258 0 : IPosition dts(3, 4, 32, 1024);
259 : // Discern DATA's default tile shape and use it
260 0 : const auto dminfo = _ms->dataManagerInfo();
261 0 : for (uInt i=0; i<dminfo.nfields(); ++i) {
262 0 : Record col = dminfo.asRecord(i);
263 0 : if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
264 0 : dts = IPosition(
265 0 : col.asRecord("SPEC").asArrayInt("DEFAULTTILESHAPE")
266 0 : );
267 0 : break;
268 : }
269 : }
270 : // Add the column
271 0 : String colWtSp = colName;
272 0 : TableDesc tdWtSp;
273 0 : tdWtSp.addColumn(ArrayColumnDesc<Float>(colWtSp, descName, 2));
274 0 : TiledShapeStMan wtSpStMan(mgrName, dts);
275 0 : _ms->addColumn(tdWtSp, wtSpStMan);
276 : }
277 : }
278 0 : else if (mustWriteNonSpec) {
279 0 : 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 0 : 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 0 : }
291 :
292 : }
|
