Line data Source code
1 : //# SDMSManager.cc: this defines single dish MS transform manager
2 : //# inheriting MSTransformManager.
3 : //#
4 : //# Copyright (C) 2015
5 : //# National Astronomical Observatory of Japan
6 : //#
7 : //# This library is free software; you can redistribute it and/or modify it
8 : //# under the terms of the GNU Library General Public License as published by
9 : //# the Free Software Foundation; either version 2 of the License, or (at your
10 : //# option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful, but WITHOUT
13 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15 : //# License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Library General Public License
18 : //# along with this library; if not, write to the Free Software Foundation,
19 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
20 : //#
21 : //# Correspondence concerning AIPS++ should be addressed as follows:
22 : //# Internet email: aips2-request@nrao.edu.
23 : //# Postal address: AIPS++ Project Office
24 : //# National Radio Astronomy Observatory
25 : //# 520 Edgemont Road
26 : //# Charlottesville, VA 22903-2475 USA
27 : //#
28 : //# $Id$
29 : #include <algorithm>
30 : #include <array>
31 : #include <iostream>
32 : #include <iterator>
33 : #include <list>
34 :
35 : //#include <libsakura/sakura.h>
36 : //#include <libsakura/config.h>
37 :
38 : #include <casacore/casa/Logging/LogIO.h>
39 : #include <casacore/casa/Logging/LogOrigin.h>
40 : #include <casacore/casa/Utilities/Assert.h>
41 : #include <casacore/casa/Arrays/ArrayMath.h>
42 : #include <casacore/casa/Utilities/Sort.h>
43 :
44 : #include <casacore/ms/MSSel/MSSelectionTools.h>
45 : #include <msvis/MSVis/VisibilityIterator2.h>
46 : #include <msvis/MSVis/VisSetUtil.h>
47 :
48 : #include <casa_sakura/SakuraUtils.h>
49 : #include <singledish/SingleDish/SDMSManager.h>
50 :
51 :
52 : #define _ORIGIN LogOrigin("SDMSManager", __func__, WHERE)
53 :
54 : using namespace casacore;
55 :
56 : namespace casa {
57 :
58 0 : SDMSManager::SDMSManager() : doSmoothing_(false) {
59 0 : }
60 :
61 : // SDMSManager &SDMSManager::operator=(SDMSManager const &other)
62 : // {
63 : // return *this;
64 : // }
65 :
66 0 : SDMSManager::~SDMSManager() {
67 0 : LogIO os(_ORIGIN);
68 0 : }
69 :
70 : // -----------------------------------------------------------------------
71 : // Fill output MS with data from an input VisBuffer
72 : // -----------------------------------------------------------------------
73 0 : void SDMSManager::fillCubeToOutputMs(vi::VisBuffer2 *vb,Cube<Float> const &data_cube,
74 : Cube<Bool> const *flag_cube, Matrix<Float> const *weight_matrix) {
75 0 : setupBufferTransformations(vb);
76 :
77 : // make sure the shape of cube matches the number of rows to add.
78 0 : AlwaysAssert(data_cube.nplane() == nRowsToAdd_p, AipsError);
79 :
80 0 : if (bufferMode_p) {
81 0 : return;
82 : }
83 :
84 : // Create RowRef object to fill new rows
85 0 : uInt currentRows = outputMs_p->nrow();
86 0 : RefRows rowRef(currentRows, currentRows + nRowsToAdd_p / nspws_p - 1);
87 :
88 : // Add new rows to output MS
89 0 : outputMs_p->addRow(nRowsToAdd_p, false);
90 :
91 : // Fill new rows
92 0 : if (weight_matrix == nullptr) {
93 0 : weightSpectrumFlatFilled_p = false;
94 0 : weightSpectrumFromSigmaFilled_p = false;
95 0 : fillWeightCols(vb, rowRef);
96 : }
97 0 : fillCubeToDataCols(vb, rowRef, data_cube, flag_cube);
98 0 : fillIdCols(vb, rowRef);
99 0 : if (weight_matrix != nullptr) { // for update_weight=True (CAS-13161)
100 0 : outputMsCols_p->weight().putColumnCells(rowRef, *weight_matrix);
101 : }
102 : }
103 :
104 : // ----------------------------------------------------------------------------------------
105 : // Fill main (data) columns which have to be combined together to produce bigger SPWs
106 : // ----------------------------------------------------------------------------------------
107 0 : void SDMSManager::fillCubeToDataCols(vi::VisBuffer2 *vb, RefRows &rowRef,
108 : Cube<Float> const &data_cube, Cube<Bool> const *flag_cube) {
109 0 : ArrayColumn<Bool> *outputFlagCol = NULL;
110 0 : for (dataColMap::iterator iter = dataColMap_p.begin(); iter != dataColMap_p.end(); iter++) {
111 : // Get applicable *_SPECTRUM (copy constructor uses reference semantics)
112 : // If channel average or combine, otherwise no need to copy
113 0 : const Cube<Float> applicableSpectrum = getApplicableSpectrum(vb, iter->first);
114 :
115 : // Apply transformations
116 0 : switch (iter->first) {
117 0 : case MS::DATA:
118 : {
119 0 : if (mainColumn_p == MS::DATA) {
120 0 : outputFlagCol = &(outputMsCols_p->flag());
121 0 : setTileShape(rowRef, outputMsCols_p->flag());
122 : } else {
123 0 : outputFlagCol = NULL;
124 : }
125 0 : setTileShape(rowRef, outputMsCols_p->data());
126 0 : Cube<Complex> cdata_cube(data_cube.shape());
127 0 : convertArray(cdata_cube, data_cube);
128 0 : transformCubeOfData(vb, rowRef, cdata_cube, outputMsCols_p->data(),
129 : outputFlagCol, applicableSpectrum);
130 0 : break;
131 : }
132 0 : case MS::CORRECTED_DATA:
133 : {
134 0 : if (mainColumn_p == MS::CORRECTED_DATA) {
135 0 : outputFlagCol = &(outputMsCols_p->flag());
136 0 : setTileShape(rowRef, outputMsCols_p->flag());
137 : } else {
138 0 : outputFlagCol = NULL;
139 : }
140 0 : Cube<Complex> cdata_cube(data_cube.shape());
141 0 : convertArray(cdata_cube, data_cube);
142 0 : if (iter->second == MS::DATA) {
143 0 : setTileShape(rowRef, outputMsCols_p->data());
144 0 : transformCubeOfData(vb, rowRef, cdata_cube, outputMsCols_p->data(),
145 : outputFlagCol, applicableSpectrum);
146 : } else {
147 0 : setTileShape(rowRef, outputMsCols_p->correctedData());
148 0 : transformCubeOfData(vb, rowRef, cdata_cube, outputMsCols_p->correctedData(),
149 : outputFlagCol, applicableSpectrum);
150 : }
151 0 : break;
152 : }
153 0 : case MS::MODEL_DATA:
154 : {
155 0 : if (mainColumn_p == MS::MODEL_DATA) {
156 0 : outputFlagCol = &(outputMsCols_p->flag());
157 0 : setTileShape(rowRef, outputMsCols_p->flag());
158 : } else {
159 0 : outputFlagCol = NULL;
160 : }
161 :
162 0 : if (iter->second == MS::DATA) {
163 0 : setTileShape(rowRef, outputMsCols_p->data());
164 0 : transformCubeOfData(vb, rowRef, vb->visCubeModel(), outputMsCols_p->data(),
165 : outputFlagCol, applicableSpectrum);
166 : } else {
167 0 : setTileShape(rowRef, outputMsCols_p->modelData());
168 0 : transformCubeOfData(vb, rowRef, vb->visCubeModel(), outputMsCols_p->modelData(),
169 : outputFlagCol, applicableSpectrum);
170 : }
171 0 : break;
172 : }
173 0 : case MS::FLOAT_DATA:
174 : {
175 0 : if (mainColumn_p == MS::FLOAT_DATA) {
176 0 : outputFlagCol = &(outputMsCols_p->flag());
177 0 : setTileShape(rowRef, outputMsCols_p->flag());
178 : } else {
179 0 : outputFlagCol = NULL;
180 : }
181 0 : setTileShape(rowRef, outputMsCols_p->floatData());
182 0 : transformCubeOfData(vb, rowRef, data_cube, outputMsCols_p->floatData(),
183 : outputFlagCol, applicableSpectrum);
184 0 : break;
185 : }
186 0 : case MS::LAG_DATA:
187 : {
188 : // jagonzal: TODO
189 0 : break;
190 : }
191 0 : default:
192 : {
193 : // jagonzal: TODO
194 0 : break;
195 : }
196 : } // end switch
197 :
198 : //KS: THIS PART ASSUMES INROW==OUTROW
199 0 : if ((outputFlagCol != NULL) && (flag_cube != nullptr)) {
200 0 : setTileShape(rowRef, outputMsCols_p->flag());
201 0 : writeCube(*flag_cube, *outputFlagCol, rowRef);
202 : }
203 : } // end loop for dataColMap
204 :
205 : // Special case for flag category
206 0 : if (inputFlagCategoryAvailable_p) {
207 0 : if (spectrumReshape_p) {
208 0 : IPosition transformedCubeShape = getShape(); //[nC,nF,nR]
209 0 : IPosition inputFlagCategoryShape = vb->flagCategory().shape(); // [nC,nF,nCategories,nR]
210 : IPosition flagCategoryShape(4,
211 0 : inputFlagCategoryShape(1),
212 0 : transformedCubeShape(2),
213 0 : inputFlagCategoryShape(2),
214 0 : transformedCubeShape(2));
215 0 : Array<Bool> flagCategory(flagCategoryShape, false);
216 :
217 0 : outputMsCols_p->flagCategory().putColumnCells(rowRef, flagCategory);
218 : } else {
219 0 : outputMsCols_p->flagCategory().putColumnCells(rowRef, vb->flagCategory());
220 : }
221 : }
222 :
223 0 : return;
224 : }
225 :
226 : // -----------------------------------------------------------------------
227 : // Method to determine sort columns order
228 : // -----------------------------------------------------------------------
229 0 : void SDMSManager::setSortColumns(Block<Int> sortColumns,
230 : bool addDefaultSortCols,
231 : Double timebin) {
232 0 : size_t const num_in = sortColumns.nelements();
233 0 : LogIO os(_ORIGIN);
234 0 : os << "Setting user sort columns with " << num_in << " elements" << LogIO::POST;
235 0 : if (addDefaultSortCols) {
236 0 : Block<Int> defaultCols(6);
237 0 : defaultCols[0] = MS::OBSERVATION_ID;
238 0 : defaultCols[1] = MS::ARRAY_ID;
239 : // Assuming auto-correlation here
240 0 : defaultCols[2] = MS::ANTENNA1;
241 0 : defaultCols[3] = MS::FEED1;
242 0 : defaultCols[4] = MS::DATA_DESC_ID;
243 0 : defaultCols[5] = MS::TIME;
244 0 : os << "Adding default sort columns to user sort column" << LogIO::POST;
245 0 : if (num_in > 0) {
246 0 : size_t num_elem(num_in);
247 0 : sortColumns.resize(num_in + defaultCols.nelements(), false, true);
248 0 : for (size_t i = 0; i < defaultCols.nelements(); ++i) {
249 0 : bool addcol = true;
250 0 : if (getBlockId(sortColumns, defaultCols[i]) > -1) {
251 0 : addcol = false; // the columns is in sortColumns
252 : }
253 0 : if (addcol) {
254 0 : sortColumns[num_elem] = defaultCols[i];
255 0 : ++num_elem;
256 : }
257 : }
258 : // shrink block
259 0 : sortColumns.resize(num_elem, true, true);
260 0 : userSortCols_ = sortColumns;
261 : } else {
262 0 : userSortCols_ = defaultCols;
263 : }
264 : } else { // do not add
265 0 : if (num_in > 0) {
266 0 : userSortCols_ = sortColumns;
267 : } else {
268 0 : userSortCols_ = Block<Int>();
269 : }
270 : }
271 0 : os << "Defined user sort columns with " << userSortCols_.nelements() << " elements" << LogIO::POST;
272 0 : timeBin_p = timebin;
273 0 : os << "Time bin is " << timeBin_p << " sec" << LogIO::POST;
274 : // regenerate iterator if necessary
275 0 : if (visibilityIterator_p != NULL) {
276 0 : os << "Regenerating iterator" << LogIO::POST;
277 0 : setIterationApproach();
278 0 : generateIterator();
279 : }
280 :
281 0 : return;
282 : }
283 :
284 0 : int SDMSManager::getBlockId(Block<Int> const &data, Int const value) {
285 0 : for (size_t i = 0; i < data.nelements(); ++i) {
286 0 : if (data[i] == value) {
287 0 : return (int)i;
288 : }
289 : }
290 :
291 0 : return -1;
292 : }
293 :
294 0 : void SDMSManager::setIterationApproach() {
295 0 : if (userSortCols_.nelements() == 0) {
296 0 : sortColumns_p = Block<Int>();
297 :
298 0 : return;
299 : }
300 :
301 : // User column is set.
302 0 : logger_p.origin(_ORIGIN);
303 :
304 : using ColumnId = MSMainEnums::PredefinedColumns;
305 : struct CheckItem {
306 0 : CheckItem(ColumnId const &id, String const &name, String const &type, Bool const &remove)
307 0 : : columnId{id}, columnName{name}, paramType{type}, removeIt{remove}, addIt{!remove}
308 0 : {}
309 : ColumnId columnId;
310 : String columnName;
311 : String paramType;
312 : Bool removeIt;
313 : Bool addIt;
314 : };
315 0 : auto const userSortColExists = [&](ColumnId const &columnId) {
316 0 : return getBlockId(userSortCols_, columnId) > -1;
317 0 : };
318 :
319 : // copy userSortCols_ to std::list
320 0 : std::list<ColumnId> userSortColsList;
321 : std::transform(
322 : userSortCols_.begin(), userSortCols_.end(),
323 : std::back_inserter(userSortColsList),
324 0 : [](Int const &i) {return static_cast<ColumnId>(i);}
325 0 : );
326 :
327 : // list of columns that may be added/removed depending on the value of timespan_p
328 : std::array<CheckItem, 3> const checkList {{
329 0 : {MS::SCAN_NUMBER, "SCAN_NUMBER", "scan", timespan_p.contains("scan")},
330 0 : {MS::STATE_ID, "STATE_ID", "state", timespan_p.contains("state")},
331 0 : {MS::FIELD_ID, "FIELD_ID", "field", timespan_p.contains("field")},
332 0 : }};
333 :
334 : // remove columns from userSortColsList if necessary
335 0 : for (auto const &item: checkList) {
336 0 : if (item.removeIt && userSortColExists(item.columnId)) {
337 0 : logger_p << LogIO::NORMAL;
338 0 : logger_p << "Combining data through " << item.paramType << "s for time average. ";
339 0 : logger_p << "Removing " << item.columnName << " from user sort list." << LogIO::POST;
340 0 : userSortColsList.remove(item.columnId);
341 : }
342 : }
343 :
344 : // add columns to userSortColsList if necessary
345 0 : if (timeAverage_p) {
346 0 : for (auto const &item: checkList) {
347 0 : if (item.addIt && !userSortColExists(item.columnId)) {
348 0 : logger_p << LogIO::NORMAL
349 0 : << "Splitting data by " << item.paramType << "s for time average. "
350 0 : << "Adding " << item.columnName << " to user sort list." << LogIO::POST;
351 0 : userSortColsList.push_back(item.columnId);
352 : }
353 : }
354 : }
355 :
356 : // copy back userSortColsList to sortColumns_p
357 0 : constexpr bool forceSmaller = true;
358 0 : constexpr bool copyElements = false;
359 0 : sortColumns_p.resize(userSortColsList.size(), forceSmaller, copyElements);
360 0 : std::transform(
361 : userSortColsList.begin(), userSortColsList.end(), sortColumns_p.begin(),
362 0 : [](ColumnId const &i) {return static_cast<Int>(i);}
363 : );
364 :
365 0 : ostringstream oss;
366 0 : for (size_t i = 0; i < sortColumns_p.nelements(); ++i) {
367 0 : oss << sortColumns_p[i] << ", ";
368 : }
369 0 : logger_p << LogIO::DEBUG1 << "Final Sort order: "
370 0 : << oss.str() << LogIO::POST;
371 :
372 0 : return;
373 : }
374 :
375 0 : Record SDMSManager::getSelRec(string const &spw) {
376 0 : MeasurementSet myms = MeasurementSet(inpMsName_p);
377 0 : return myms.msseltoindex(toCasaString(spw));
378 : }
379 :
380 : /*
381 : MeasurementSet SDMSManager::getMS()
382 : {
383 : MeasurementSet myms = MeasurementSet(inpMsName_p);
384 : return myms;
385 : }
386 : */
387 :
388 0 : void SDMSManager::setSmoothing(string const &kernelType, float const &kernelWidth) {
389 : // kernel type
390 0 : VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernelType);
391 :
392 : // Fail if type is neither GAUSSIAN nor BOXCAR since other ones are yet to be supported
393 0 : if ((type != VectorKernel::GAUSSIAN) && (type != VectorKernel::BOXCAR)) {
394 0 : stringstream oss;
395 0 : oss << "Smoothing kernel type \"" << kernelType << "\" is not supported yet.";
396 0 : throw AipsError(oss.str());
397 : }
398 :
399 : // Fail if kernel width is zero or negative
400 0 : if (kernelWidth <= 0.0) {
401 0 : throw AipsError("Zero or negative kernel width is not allowed.");
402 : }
403 :
404 0 : doSmoothing_ = true;
405 0 : kernelType_ = type;
406 0 : kernelWidth_ = kernelWidth;
407 0 : }
408 :
409 0 : void SDMSManager::unsetSmoothing() {
410 0 : doSmoothing_ = false;
411 0 : }
412 :
413 0 : void SDMSManager::initializeSmoothing() {
414 0 : LogIO os(_ORIGIN);
415 0 : if (!doSmoothing_) {
416 0 : return;
417 : }
418 :
419 0 : Vector<Int> numChanList = inspectNumChan();
420 0 : for (size_t i = 0; i < numChanList.nelements(); ++i) {
421 0 : Int numChan = numChanList[i];
422 0 : Vector<Float> theKernel = VectorKernel::make(kernelType_, kernelWidth_, numChan, true, false);
423 : //shift 1 channel for boxcar kernel---(for CAS-7442, 2015/11/18 WK)
424 0 : if (kernelType_ == VectorKernel::BOXCAR) {
425 0 : for (size_t j = theKernel.nelements()-1; j >= 1; --j) {
426 0 : theKernel[j] = theKernel[j-1];
427 : }
428 0 : theKernel[0] = 0.0;
429 : }
430 : //------
431 0 : convolverPool_[numChan] = Convolver<Float>();
432 0 : convolverPool_[numChan].setPsf(theKernel, IPosition(1, numChan));
433 : }
434 :
435 : // set smoothing kernel for weight
436 0 : smoothBin_p = static_cast<uInt>(kernelWidth_ + 0.5f);
437 0 : smoothBin_p += (smoothBin_p % 2 == 0) ? 1 : 0; // to make smoothBin_p odd
438 0 : uInt numChanMinimum = min(numChanList);
439 0 : smoothBin_p = min(smoothBin_p, numChanMinimum - ((numChanMinimum % 2 == 0) ? 1 : 0)); // smoothBin_p < numChanMinimum
440 0 : uInt halfWidth = smoothBin_p / 2;
441 0 : Sort sort;
442 0 : Vector<Float> kernelForMinimumNumChan = convolverPool_[numChanMinimum].getPsf();
443 0 : sort.sortKey(kernelForMinimumNumChan.data(), TpFloat, 0, Sort::Descending);
444 0 : Vector<uInt> indexArray;
445 : //uInt indexArrayLength = sort.sort(indexArray, numChanMinimum);
446 0 : sort.sort(indexArray, numChanMinimum);
447 0 : uInt startChan = indexArray[0] - halfWidth;
448 0 : uInt endChan = startChan + smoothBin_p;
449 0 : smoothCoeff_p.resize(smoothBin_p, false);
450 0 : for (uInt i = startChan, j = 0; i < endChan; ++i, ++j) {
451 0 : smoothCoeff_p[j] = kernelForMinimumNumChan[i];
452 : }
453 : // normalize smoothCoeff_p
454 0 : smoothCoeff_p /= sum(smoothCoeff_p);
455 0 : os << LogIO::DEBUGGING << "smoothBin_p = " << smoothBin_p << LogIO::POST;
456 0 : os << LogIO::DEBUGGING << "smoothCoeff_p = " << smoothCoeff_p << LogIO::POST;
457 : }
458 :
459 0 : Vector<Int> SDMSManager::inspectNumChan() {
460 0 : LogIO os(_ORIGIN);
461 :
462 0 : if (selectedInputMs_p == NULL) {
463 0 : throw AipsError("Input MS is not opened yet.");
464 : }
465 :
466 0 : ScalarColumn<Int> col(*selectedInputMs_p, "DATA_DESC_ID");
467 0 : Vector<Int> ddIdList = col.getColumn();
468 0 : uInt numDDId = GenSort<Int>::sort(ddIdList, Sort::Ascending,
469 : Sort::QuickSort | Sort::NoDuplicates);
470 0 : col.attach(selectedInputMs_p->dataDescription(), "SPECTRAL_WINDOW_ID");
471 0 : Vector<Int> spwIdList(numDDId);
472 0 : for (uInt i = 0; i < numDDId; ++i) {
473 0 : spwIdList[i] = col(ddIdList[i]);
474 : }
475 0 : Vector<Int> numChanList(numDDId);
476 0 : col.attach(selectedInputMs_p->spectralWindow(), "NUM_CHAN");
477 0 : os << LogIO::DEBUGGING << "spwIdList = " << spwIdList << LogIO::POST;
478 0 : for (size_t i = 0; i < spwIdList.nelements(); ++i) {
479 0 : Int spwId = spwIdList[i];
480 0 : Int numChan = col(spwId);
481 0 : numChanList[i] = numChan;
482 0 : os << LogIO::DEBUGGING << "examine spw " << spwId << ": nchan = " << numChan << LogIO::POST;
483 0 : if (numChan == 1) {
484 0 : stringstream ss;
485 0 : ss << "smooth: Failed due to wrong spw " << i;
486 0 : throw AipsError(ss.str());
487 : }
488 : }
489 :
490 0 : uInt numNumChan = GenSort<Int>::sort(numChanList, Sort::Ascending,
491 : Sort::QuickSort | Sort::NoDuplicates);
492 0 : return numChanList(Slice(0, numNumChan));
493 : }
494 :
495 : } // End of casa namespace.
|
