Line data Source code
1 : //# SingleDishSkyCal.cc: implements SingleDishSkyCal
2 : //# Copyright (C) 2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library is distributed in the hope that it will be useful, but WITHOUT
11 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed as follows:
20 : //# Internet email: casa-feedback@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //# $Id$
27 :
28 : //# Includes
29 : #include <iostream>
30 : #include <sstream>
31 : #include <iomanip>
32 : #include <map>
33 : #include <type_traits>
34 :
35 : #include <casacore/casa/OS/File.h>
36 : #include <casacore/casa/Arrays/MaskedArray.h>
37 : #include <casacore/casa/Arrays/MaskArrMath.h>
38 : #include <casacore/tables/TaQL/TableParse.h>
39 : #include <casacore/tables/Tables/Table.h>
40 : #include <casacore/tables/Tables/ScalarColumn.h>
41 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
42 : #include <casacore/ms/MeasurementSets/MSState.h>
43 : #include <casacore/ms/MeasurementSets/MSSpectralWindow.h>
44 : #include <casacore/ms/MeasurementSets/MSIter.h>
45 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
46 : #include <synthesis/MeasurementComponents/SingleDishSkyCal.h>
47 : #include <synthesis/CalTables/CTGlobals.h>
48 : #include <synthesis/CalTables/CTMainColumns.h>
49 : #include <synthesis/CalTables/CTInterface.h>
50 : #include <casacore/ms/MSSel/MSSelection.h>
51 : #include <casacore/ms/MSSel/MSSelectionTools.h>
52 : #include <synthesis/Utilities/PointingDirectionCalculator.h>
53 : #include <synthesis/Utilities/PointingDirectionProjector.h>
54 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
55 : #include <casa_sakura/SakuraAlignedArray.h>
56 : #include <libsakura/sakura.h>
57 : #include <cassert>
58 :
59 : // Debug Message Handling
60 : // if SDCALSKY_DEBUG is defined, the macro debuglog and
61 : // debugpost point standard output stream (std::cout and
62 : // std::endl so that debug messages are sent to standard
63 : // output. Otherwise, these macros basically does nothing.
64 : // "Do nothing" behavior is implemented in NullLogger
65 : // and its associating << operator below.
66 : //
67 : // Usage:
68 : // Similar to standard output stream.
69 : //
70 : // debuglog << "Any message" << any_value << debugpost;
71 : //
72 : //#define SDCALSKY_DEBUG
73 :
74 : using namespace casacore;
75 :
76 : namespace {
77 : struct NullLogger {};
78 :
79 : template<class T>
80 0 : inline NullLogger &operator<<(NullLogger &logger, T /*value*/) {
81 0 : return logger;
82 : }
83 :
84 : #ifndef SDCALSKY_DEBUG
85 : NullLogger nulllogger;
86 : #endif
87 :
88 : // Helper class to format integers with thousand separators
89 : struct CommaSeparatedThousands : std::numpunct<char> {
90 0 : char_type do_thousands_sep() const override { return ','; }
91 0 : string_type do_grouping() const override { return "\3"; }
92 : };
93 :
94 : }
95 :
96 : #ifdef SDCALSKY_DEBUG
97 : #define debuglog std::cout
98 : #define debugpost std::endl
99 : #else
100 : #define debuglog nulllogger
101 : #define debugpost 0
102 : #endif
103 :
104 : // Local Functions
105 : namespace {
106 : inline Vector<rownr_t> getOffStateIdList(MeasurementSet const &ms) {
107 : String taql("SELECT FLAG_ROW FROM $1 WHERE UPPER(OBS_MODE) ~ m/^OBSERVE_TARGET#OFF_SOURCE/");
108 : Table stateSel = tableCommand(taql, ms.state()).table();
109 : Vector<rownr_t> stateIdList = stateSel.rowNumbers();
110 : debuglog << "stateIdList[" << stateIdList.nelements() << "]=";
111 : for (size_t i = 0; i < stateIdList.nelements(); ++i) {
112 : debuglog << stateIdList[i] << " ";
113 : }
114 : debuglog << debugpost;
115 : return stateIdList;
116 : }
117 :
118 : template<class T>
119 0 : inline std::string toString(casacore::Vector<T> const &v) {
120 0 : std::ostringstream oss;
121 0 : oss << "[";
122 0 : std::string delimiter = "";
123 0 : for (size_t i = 0; i < v.nelements(); ++i) {
124 0 : oss << delimiter << v[i];
125 0 : delimiter = ",";
126 : }
127 0 : oss << "]";
128 0 : return oss.str();
129 0 : }
130 :
131 : // unused
132 : /*
133 : inline casacore::String configureTaqlString(casacore::String const &msName, casacore::Vector<casacore::uInt> stateIdList) {
134 : std::ostringstream oss;
135 : oss << "SELECT FROM " << msName << " WHERE ANTENNA1 == ANTENNA2 && STATE_ID IN "
136 : << toString(stateIdList)
137 : << " ORDER BY FIELD_ID, ANTENNA1, FEED1, DATA_DESC_ID, TIME";
138 : return casacore::String(oss);
139 : }
140 : */
141 :
142 0 : inline void fillNChanParList(casacore::MeasurementSet const &ms, casacore::Vector<casacore::Int> &nChanParList) {
143 0 : casacore::MSSpectralWindow const &msspw = ms.spectralWindow();
144 0 : casacore::ScalarColumn<casacore::Int> nchanCol(msspw, "NUM_CHAN");
145 0 : debuglog << "nchanCol=" << toString(nchanCol.getColumn()) << debugpost;
146 0 : nChanParList = nchanCol.getColumn()(casacore::Slice(0,nChanParList.nelements()));
147 0 : debuglog << "nChanParList=" << nChanParList << debugpost;
148 0 : }
149 :
150 0 : inline void updateWeight(casa::NewCalTable &ct) {
151 0 : casa::CTMainColumns ctmc(ct);
152 0 : casacore::ArrayColumn<casacore::Double> chanWidthCol(ct.spectralWindow(), "CHAN_WIDTH");
153 0 : casacore::Vector<casacore::Int> chanWidth(chanWidthCol.nrow());
154 0 : for (size_t irow = 0; irow < chanWidthCol.nrow(); ++irow) {
155 0 : casacore::Double chanWidthVal = chanWidthCol(irow)(casacore::IPosition(1,0));
156 0 : chanWidth[irow] = abs(chanWidthVal);
157 : }
158 0 : for (size_t irow = 0; irow < ct.nrow(); ++irow) {
159 0 : casacore::Int spwid = ctmc.spwId()(irow);
160 0 : casacore::Double width = chanWidth[spwid];
161 0 : casacore::Double interval = ctmc.interval()(irow);
162 0 : casacore::Float weight = width * interval;
163 0 : casacore::Matrix<casacore::Float> weightMat(ctmc.fparam().shape(irow), weight);
164 0 : ctmc.weight().put(irow, weightMat);
165 0 : }
166 0 : }
167 :
168 : class DataColumnAccessor {
169 : public:
170 0 : DataColumnAccessor(casacore::Table const &ms,
171 : casacore::String const colName="DATA")
172 0 : : dataCol_(ms, colName) {
173 0 : }
174 0 : casacore::Matrix<casacore::Float> operator()(size_t irow) {
175 0 : return casacore::real(dataCol_(irow));
176 : }
177 : private:
178 : DataColumnAccessor() {}
179 : casacore::ArrayColumn<casacore::Complex> dataCol_;
180 : };
181 :
182 : class FloatDataColumnAccessor {
183 : public:
184 0 : FloatDataColumnAccessor(casacore::Table const &ms)
185 0 : : dataCol_(ms, "FLOAT_DATA") {
186 0 : }
187 0 : casacore::Matrix<casacore::Float> operator()(size_t irow) {
188 0 : return dataCol_(irow);
189 : }
190 : private:
191 : FloatDataColumnAccessor() {}
192 : casacore::ArrayColumn<casacore::Float> dataCol_;
193 : };
194 :
195 0 : inline casacore::Int nominalDataDesc(casacore::MeasurementSet const &ms, casacore::Int const ant)
196 : {
197 0 : casacore::Int goodDataDesc = -1;
198 0 : casacore::ScalarColumn<casacore::Int> col(ms.spectralWindow(), "NUM_CHAN");
199 0 : casacore::Vector<casacore::Int> nchanList = col.getColumn();
200 0 : size_t numSpw = col.nrow();
201 0 : casacore::Vector<casacore::Int> spwMap(numSpw);
202 0 : col.attach(ms.dataDescription(), "SPECTRAL_WINDOW_ID");
203 0 : for (size_t i = 0; i < col.nrow(); ++i) {
204 0 : spwMap[col(i)] = i;
205 : }
206 0 : casacore::ScalarColumn<casacore::String> obsModeCol(ms.state(), "OBS_MODE");
207 0 : casacore::Regex regex("^OBSERVE_TARGET#ON_SOURCE.*");
208 0 : for (size_t ispw = 0; ispw < numSpw; ++ispw) {
209 0 : if (nchanList[ispw] == 4) {
210 : // this should be WVR
211 0 : continue;
212 : }
213 : else {
214 0 : std::ostringstream oss;
215 0 : oss << "SELECT FROM $1 WHERE ANTENNA1 == " << ant
216 0 : << " && ANTENNA2 == " << ant << " && DATA_DESC_ID == " << spwMap[ispw]
217 0 : << " ORDER BY STATE_ID";
218 0 : casacore::MeasurementSet msSel(casacore::tableCommand(oss.str(), ms).table());
219 0 : col.attach(msSel, "STATE_ID");
220 0 : casacore::Vector<casacore::Int> stateIdList = col.getColumn();
221 0 : casacore::Int previousStateId = -1;
222 0 : for (size_t i = 0; i < msSel.nrow(); ++i) {
223 0 : casacore::Int stateId = stateIdList[i];
224 0 : if (stateId != previousStateId) {
225 0 : casacore::String obsmode = obsModeCol(stateId);
226 0 : if (regex.match(obsmode.c_str(), obsmode.size()) != casacore::String::npos) {
227 0 : goodDataDesc = spwMap[ispw];
228 0 : break;
229 : }
230 0 : }
231 0 : previousStateId = stateId;
232 : }
233 0 : }
234 :
235 0 : if (goodDataDesc >= 0)
236 0 : break;
237 : }
238 0 : return goodDataDesc;
239 0 : }
240 :
241 0 : inline casacore::Vector<size_t> detectGap(casacore::Vector<casacore::Double> timeList)
242 : {
243 0 : size_t n = timeList.size();
244 0 : casacore::Vector<casacore::Double> timeInterval = timeList(casacore::Slice(1, n-1)) - timeList(casacore::Slice(0, n-1));
245 0 : casacore::Double medianTime = casacore::median(timeInterval);
246 0 : casacore::Double const threshold = medianTime * 5.0;
247 0 : casacore::Vector<size_t> gapIndexList(casacore::IPosition(1, n/2 + 2), new size_t[n/2+2], casacore::TAKE_OVER);
248 0 : gapIndexList[0] = 0;
249 0 : size_t gapIndexCount = 1;
250 0 : for (size_t i = 0; i < timeInterval.size(); ++i) {
251 0 : if (timeInterval[i] > threshold) {
252 0 : gapIndexList[gapIndexCount] = i + 1;
253 0 : gapIndexCount++;
254 : }
255 : }
256 0 : if (gapIndexList[gapIndexCount] != n) {
257 0 : gapIndexList[gapIndexCount] = n;
258 0 : gapIndexCount++;
259 : }
260 0 : debuglog << "Detected " << gapIndexCount << " gaps." << debugpost;
261 0 : casacore::Vector<size_t> ret(casacore::IPosition(1, gapIndexCount), gapIndexList.data(), casacore::COPY);
262 0 : debuglog << "gapList=" << toString(ret) << debugpost;
263 0 : return ret;
264 0 : }
265 :
266 : struct DefaultRasterEdgeDetector
267 : {
268 0 : static size_t N(size_t numData, casacore::Float const /*fraction*/, casacore::Int const /*num*/)
269 : {
270 0 : debuglog << "DefaultRasterEdgeDetector" << debugpost;
271 0 : return max((size_t)1, static_cast<size_t>(sqrt(numData + 1) - 1));
272 : }
273 : };
274 :
275 : struct FixedNumberRasterEdgeDetector
276 : {
277 0 : static size_t N(size_t /*numData*/, casacore::Float const /*fraction*/, casacore::Int const num)
278 : {
279 0 : debuglog << "FixedNumberRasterEdgeDetector" << debugpost;
280 0 : if (num < 0) {
281 0 : throw casacore::AipsError("Negative number of edge points.");
282 : }
283 0 : return (size_t)num;
284 : }
285 : };
286 :
287 : struct FixedFractionRasterEdgeDetector
288 : {
289 0 : static casacore::Int N(size_t numData, casacore::Float const fraction, casacore::Int const /*num*/)
290 : {
291 0 : debuglog << "FixedFractionRasterEdgeDetector" << debugpost;
292 0 : return max((size_t)1, static_cast<size_t>(numData * fraction));
293 : }
294 : };
295 :
296 : template<class Detector>
297 0 : inline casacore::Vector<casacore::Double> detectEdge(casacore::Vector<casacore::Double> timeList, casacore::Double const interval, casacore::Float const fraction, casacore::Int const num)
298 : {
299 : // storage for time range for edges (at head and tail)
300 : // [(start of head edge), (end of head edge),
301 : // (start of tail edge), (end of tail edge)]
302 0 : casacore::Vector<casacore::Double> edgeList(4);
303 0 : size_t numList = timeList.size();
304 0 : size_t numEdge = Detector::N(numList, fraction, num);
305 0 : debuglog << "numEdge = " << numEdge << debugpost;
306 0 : if (numEdge == 0) {
307 0 : throw casacore::AipsError("Zero edge points.");
308 : }
309 0 : else if (timeList.size() > numEdge * 2) {
310 0 : edgeList[0] = timeList[0] - 0.5 * interval;
311 0 : edgeList[1] = timeList[numEdge-1] + 0.5 * interval;
312 0 : edgeList[2] = timeList[numList-numEdge] - 0.5 * interval;
313 0 : edgeList[3] = timeList[numList-1] + 0.5 * interval;
314 : }
315 : else {
316 0 : std::ostringstream oss;
317 0 : oss << "Too many edge points (" << 2.0 * numEdge << " out of "
318 0 : << timeList.size() << " points)";
319 0 : throw casacore::AipsError(oss.str());
320 : // edgeList[0] = timeList[0] - 0.5 * interval;
321 : // edgeList[1] = timeList[numList-1] + 0.5 * interval;
322 : // edgeList[2] = edgeList[0];
323 : // edgeList[3] = edgeList[2];
324 0 : }
325 0 : return edgeList;
326 0 : }
327 :
328 0 : inline casacore::Vector<casacore::String> detectRaster(casacore::MeasurementSet const &ms,
329 : casacore::Int const ant,
330 : casacore::Float const fraction,
331 : casacore::Int const num)
332 : {
333 0 : casacore::Int dataDesc = nominalDataDesc(ms, ant);
334 0 : debuglog << "nominal DATA_DESC_ID=" << dataDesc << debugpost;
335 0 : assert(dataDesc >= 0);
336 0 : if (dataDesc < 0) {
337 0 : return casacore::Vector<casacore::String>();
338 : }
339 :
340 0 : std::ostringstream oss;
341 0 : oss << "SELECT FROM $1 WHERE ANTENNA1 == " << ant
342 0 : << " && ANTENNA2 == " << ant << " && FEED1 == 0 && FEED2 == 0"
343 0 : << " && DATA_DESC_ID == " << dataDesc
344 0 : << " ORDER BY TIME";
345 0 : debuglog << "detectRaster: taql=" << oss.str() << debugpost;
346 0 : casacore::MeasurementSet msSel(casacore::tableCommand(oss.str(), ms).table());
347 0 : casacore::ScalarColumn<casacore::Double> timeCol(msSel, "TIME");
348 0 : casacore::ScalarColumn<casacore::Double> intervalCol(msSel, "INTERVAL");
349 0 : casacore::Vector<casacore::Double> timeList = timeCol.getColumn();
350 0 : casacore::Double interval = casacore::min(intervalCol.getColumn());
351 0 : casacore::Vector<size_t> gapList = detectGap(timeList);
352 0 : casacore::Vector<casacore::String> edgeAsTimeRange(gapList.size() * 2);
353 : typedef casacore::Vector<casacore::Double> (*DetectorFunc)(casacore::Vector<casacore::Double>, casacore::Double const, casacore::Float const, casacore::Int const);
354 0 : DetectorFunc detect = NULL;
355 0 : if (num > 0) {
356 0 : detect = detectEdge<FixedNumberRasterEdgeDetector>;
357 : }
358 0 : else if (fraction > 0.0) {
359 0 : detect = detectEdge<FixedFractionRasterEdgeDetector>;
360 : }
361 : else {
362 0 : detect = detectEdge<DefaultRasterEdgeDetector>;
363 : }
364 0 : for (size_t i = 0; i < gapList.size()-1; ++i) {
365 0 : size_t startRow = gapList[i];
366 0 : size_t endRow = gapList[i+1];
367 0 : size_t len = endRow - startRow;
368 0 : debuglog << "startRow=" << startRow << ", endRow=" << endRow << debugpost;
369 0 : casacore::Vector<casacore::Double> oneRow = timeList(casacore::Slice(startRow, len));
370 0 : casacore::Vector<casacore::Double> edgeList = detect(oneRow, interval, fraction, num);
371 0 : std::ostringstream s;
372 0 : s << std::setprecision(16) << "TIME BETWEEN " << edgeList[0] << " AND " << edgeList[1];
373 0 : edgeAsTimeRange[2*i] = s.str();
374 0 : s.str("");
375 0 : s << std::setprecision(16) << "TIME BETWEEN " << edgeList[2] << " AND " << edgeList[3];
376 0 : edgeAsTimeRange[2*i+1] = s.str();
377 0 : debuglog << "Resulting selection: (" << edgeAsTimeRange[2*i] << ") || ("
378 0 : << edgeAsTimeRange[2*i+1] << ")" << debugpost;
379 0 : }
380 0 : return edgeAsTimeRange;
381 0 : }
382 :
383 : // Formula for weight scaling factor, WF
384 : // 1. only one OFF spectrum is used (i.e. no interpolation)
385 : //
386 : // sigma = sqrt(sigma_ON^2 + sigma_OFF^2)
387 : // = sigma_ON * sqrt(1 + tau_ON / tau_OFF)
388 : //
389 : // weight = 1 / sigma^2
390 : // = 1 / sigma_ON^2 * tau_OFF / (tau_ON + tau_OFF)
391 : // = weight_ON * tau_OFF / (tau_ON + tau_OFF)
392 : //
393 : // WF = tau_OFF / (tau_ON + tau_OFF)
394 : //
395 : struct SimpleWeightScalingScheme
396 : {
397 0 : static casacore::Float SimpleScale(casacore::Double exposureOn, casacore::Double exposureOff)
398 : {
399 0 : return exposureOff / (exposureOn + exposureOff);
400 : }
401 : };
402 : // 2. two OFF spectrum is used (linear interpolation)
403 : //
404 : // sigma_OFF = {(t_OFF2 - t_ON)^2 * sigma_OFF1^2
405 : // + (t_ON - t_OFF1)^2 * sigma_OFF2^2}
406 : // / (t_OFF2 - t_OFF1)^2
407 : //
408 : // sigma = sqrt(sigma_ON^2 + sigma_OFF^2)
409 : // = sigma_ON * sqrt(1 + tau_ON / (t_OFF2 - t_OFF1)^2
410 : // * {(t_OFF2 - t_ON)^2 / tau_OFF1
411 : // + (t_ON - t_OFF1)^2 / tau_OFF2})
412 : //
413 : // weight = weight_ON / (1 + tau_ON / (t_OFF2 - t_OFF1)^2
414 : // * {(t_OFF2 - t_ON)^2 / tau_OFF1
415 : // + (t_ON - t_OFF1)^2 / tau_OFF2})
416 : //
417 : // WF = 1.0 / (1 + tau_ON / (t_OFF2 - t_OFF1)^2
418 : // * {(t_OFF2 - t_ON)^2 / tau_OFF1
419 : // + (t_ON - t_OFF1)^2 / tau_OFF2})
420 : //
421 : struct LinearWeightScalingScheme : public SimpleWeightScalingScheme
422 : {
423 0 : static casacore::Float InterpolateScale(casacore::Double timeOn, casacore::Double exposureOn,
424 : casacore::Double timeOff1, casacore::Double exposureOff1,
425 : casacore::Double timeOff2, casacore::Double exposureOff2)
426 : {
427 0 : casacore::Double dt = timeOff2 - timeOff1;
428 0 : casacore::Double dt1 = timeOn - timeOff1;
429 0 : casacore::Double dt2 = timeOff2 - timeOn;
430 0 : return 1.0f / (1.0f + exposureOn / (dt * dt)
431 0 : * (dt2 * dt2 / exposureOff1 + dt1 * dt1 / exposureOff2));
432 : }
433 : };
434 :
435 : // 3. two OFF spectrum is used (nearest interpolation)
436 : //
437 : // formulation is same as case 1.
438 : //
439 : struct NearestWeightScalingScheme : public SimpleWeightScalingScheme
440 : {
441 0 : static casacore::Float InterpolateScale(casacore::Double timeOn, casacore::Double exposureOn,
442 : casacore::Double timeOff1, casacore::Double exposureOff1,
443 : casacore::Double timeOff2, casacore::Double exposureOff2)
444 : {
445 0 : casacore::Double dt1 = abs(timeOn - timeOff1);
446 0 : casacore::Double dt2 = abs(timeOff2 - timeOn);
447 0 : return (dt1 <= dt2) ?
448 0 : SimpleScale(exposureOn, exposureOff1)
449 0 : : SimpleScale(exposureOn, exposureOff2);
450 : }
451 : };
452 :
453 0 : inline bool isEphemeris(casacore::String const &name) {
454 : // Check if given name is included in MDirection types
455 : casacore::Int nall, nextra;
456 : const casacore::uInt *typ;
457 0 : auto *types = casacore::MDirection::allMyTypes(nall, nextra, typ);
458 0 : auto start_extra = nall - nextra;
459 0 : auto capital_name = name;
460 0 : capital_name.upcase();
461 :
462 0 : for (auto i = start_extra; i < nall; ++i) {
463 0 : if (capital_name == types[i]) {
464 0 : return true;
465 : }
466 : }
467 :
468 0 : return false;
469 0 : }
470 :
471 : // Set number of correlations per spw
472 0 : inline void setNumberOfCorrelationsPerSpw(casacore::MeasurementSet const &ms,
473 : casacore::Vector<casacore::Int> &nCorr) {
474 0 : casacore::ScalarColumn<casacore::Int> spwIdCol(ms.dataDescription(), "SPECTRAL_WINDOW_ID");
475 0 : casacore::ScalarColumn<casacore::Int> polIdCol(ms.dataDescription(), "POLARIZATION_ID");
476 0 : casacore::ScalarColumn<casacore::Int> numCorrCol(ms.polarization(), "NUM_CORR");
477 0 : auto numSpw = ms.spectralWindow().nrow();
478 0 : auto const spwIdList = spwIdCol.getColumn();
479 0 : auto const polIdList = polIdCol.getColumn();
480 0 : auto const numCorrList = numCorrCol.getColumn();
481 0 : if (nCorr.size() != numSpw) {
482 0 : nCorr.resize(numSpw);
483 : }
484 : // initialize number of correlations with 0
485 0 : nCorr = 0;
486 0 : for (auto i = 0u; i < spwIdList.size(); ++i) {
487 0 : auto const spwId = spwIdList[i];
488 0 : auto const polId = polIdList[i];
489 0 : auto const numCorr = numCorrList[polId];
490 0 : nCorr[spwId] = numCorr;
491 : }
492 0 : }
493 : }
494 :
495 : namespace casa { //# NAMESPACE CASA - BEGIN
496 : //
497 : // SingleDishSkyCal
498 : //
499 :
500 : // Constructor
501 0 : SingleDishSkyCal::SingleDishSkyCal(VisSet& vs)
502 : : VisCal(vs),
503 : SolvableVisCal(vs),
504 0 : currAnt_(-1),
505 0 : engineC_(vs.numberSpw(), NULL),
506 0 : engineF_(vs.numberSpw(), NULL),
507 0 : currentSky_(vs.numberSpw(), NULL),
508 0 : currentSkyOK_(vs.numberSpw(), NULL),
509 0 : nCorr_(nSpw())
510 : {
511 0 : debuglog << "SingleDishSkyCal::SingleDishSkyCal(VisSet& vs)" << debugpost;
512 0 : append() = false;
513 :
514 0 : initializeSky();
515 0 : initializeCorr();
516 0 : }
517 :
518 0 : SingleDishSkyCal::SingleDishSkyCal(const MSMetaInfoForCal& msmc)
519 : : VisCal(msmc),
520 : SolvableVisCal(msmc),
521 0 : currAnt_(-1),
522 0 : engineC_(msmc.nSpw(), NULL),
523 0 : engineF_(msmc.nSpw(), NULL),
524 0 : currentSky_(msmc.nSpw(), NULL),
525 0 : currentSkyOK_(msmc.nSpw(), NULL),
526 0 : nCorr_(nSpw())
527 : {
528 0 : debuglog << "SingleDishSkyCal::SingleDishSkyCal(const MSMetaInfoForCal& msmc)" << debugpost;
529 0 : append() = False;
530 :
531 0 : initializeSky();
532 0 : initializeCorr();
533 0 : }
534 :
535 0 : SingleDishSkyCal::SingleDishSkyCal(const Int& nAnt)
536 : : VisCal(nAnt),
537 : SolvableVisCal(nAnt),
538 0 : currAnt_(-1),
539 0 : engineC_(1, NULL),
540 0 : engineF_(1, NULL),
541 0 : currentSky_(1, NULL),
542 0 : currentSkyOK_(1, NULL),
543 0 : nCorr_(nSpw())
544 : {
545 0 : debuglog << "SingleDishSkyCal::SingleDishSkyCal(const Int& nAnt)" << debugpost;
546 0 : append() = false;
547 :
548 0 : initializeSky();
549 0 : initializeCorr();
550 0 : }
551 :
552 : // Destructor
553 0 : SingleDishSkyCal::~SingleDishSkyCal()
554 : {
555 0 : debuglog << "SingleDishSkyCal::~SingleDishSkyCal()" << debugpost;
556 :
557 0 : finalizeSky();
558 0 : }
559 :
560 0 : void SingleDishSkyCal::guessPar(VisBuffer& /*vb*/)
561 : {
562 0 : }
563 :
564 0 : void SingleDishSkyCal::differentiate( CalVisBuffer & /*cvb*/)
565 : {
566 0 : }
567 :
568 0 : void SingleDishSkyCal::differentiate(VisBuffer& /*vb*/, Cube<Complex>& /*V*/,
569 : Array<Complex>& /*dV*/, Matrix<Bool>& /*Vflg*/)
570 : {
571 0 : }
572 :
573 0 : void SingleDishSkyCal::accumulate(SolvableVisCal* /*incr*/,
574 : const Vector<Int>& /*fields*/)
575 : {
576 0 : }
577 :
578 0 : void SingleDishSkyCal::diffSrc(VisBuffer& /*vb*/, Array<Complex>& /*dV*/)
579 : {
580 0 : }
581 :
582 0 : void SingleDishSkyCal::fluxscale(const String& /*outfile*/,
583 : const Vector<Int>& /*refFieldIn*/,
584 : const Vector<Int>& /*tranFieldIn*/,
585 : const Vector<Int>& /*inRefSpwMap*/,
586 : const Vector<String>& /*fldNames*/,
587 : const Float& /*inGainThres*/,
588 : const String& /*antSel*/,
589 : const String& /*timerangeSel*/,
590 : const String& /*scanSel*/,
591 : fluxScaleStruct& /*oFluxScaleStruct*/,
592 : const String& /*oListFile*/,
593 : const Bool& /*incremental*/,
594 : const Int& /*fitorder*/,
595 : const Bool& /*display*/)
596 : {
597 0 : }
598 :
599 0 : void SingleDishSkyCal::listCal(const Vector<Int> /*ufldids*/, const Vector<Int> /*uantids*/,
600 : const Matrix<Int> /*uchanids*/,
601 : //const Int& /*spw*/, const Int& /*chan*/,
602 : const String& /*listfile*/, const Int& /*pagerows*/)
603 : {
604 0 : }
605 :
606 0 : void SingleDishSkyCal::setApply(const Record& apply)
607 : {
608 : // Override interp
609 : // default frequency interpolation option is 'linearflag'
610 0 : Record applyCopy(apply);
611 0 : if (applyCopy.isDefined("interp")) {
612 0 : String interp = applyCopy.asString("interp");
613 0 : if (!interp.contains(',')) {
614 : //fInterpType() = "linearflag";
615 0 : String newInterp = interp + ",linearflag";
616 0 : applyCopy.define("interp", newInterp);
617 0 : }
618 0 : }
619 : else {
620 0 : applyCopy.define("interp", "linear,linearflag");
621 : }
622 :
623 : // call parent method
624 0 : SolvableVisCal::setApply(applyCopy);
625 0 : }
626 :
627 0 : void SingleDishSkyCal::fillCalibrationTable(casacore::MeasurementSet const &reference_data){
628 0 : String dataColName = (reference_data.tableDesc().isColumn("FLOAT_DATA")) ? "FLOAT_DATA" : "DATA";
629 :
630 0 : if ( dataColName == "FLOAT_DATA")
631 0 : fillCalibrationTable<FloatDataColumnAccessor>(reference_data);
632 : else
633 0 : fillCalibrationTable<DataColumnAccessor>(reference_data);
634 0 : }
635 :
636 : template<class DataRealComponentAccessor>
637 0 : void SingleDishSkyCal::fillCalibrationTable(MeasurementSet const &ms) {
638 0 : Int cols[] = {MS::FIELD_ID, MS::ANTENNA1, MS::FEED1,
639 : MS::DATA_DESC_ID};
640 0 : Int *colsp = cols;
641 0 : Block<Int> sortCols(4, colsp, false);
642 0 : MSIter msIter(ms, sortCols, 0.0, false, false);
643 0 : for (msIter.origin(); msIter.more(); msIter++) {
644 0 : Table const current = msIter.table();
645 0 : uInt nrow = current.nrow();
646 0 : debuglog << "nrow = " << nrow << debugpost;
647 0 : if (nrow == 0) {
648 0 : debuglog << "Skip" << debugpost;
649 0 : continue;
650 : }
651 0 : Int ispw = msIter.spectralWindowId();
652 0 : if (nChanParList()[ispw] == 4) {
653 : // Skip WVR
654 : debuglog << "Skip " << ispw
655 0 : << "(nchan " << nChanParList()[ispw] << ")"
656 0 : << debugpost;
657 0 : continue;
658 : }
659 : debuglog << "Process " << ispw
660 0 : << "(nchan " << nChanParList()[ispw] << ")"
661 0 : << debugpost;
662 :
663 0 : Int ifield = msIter.fieldId();
664 0 : ScalarColumn<Int> antennaCol(current, "ANTENNA1");
665 0 : currAnt_ = antennaCol(0);
666 0 : ScalarColumn<Int> feedCol(current, "FEED1");
667 : debuglog << "FIELD_ID " << msIter.fieldId()
668 : << " ANTENNA1 " << currAnt_
669 : << " FEED1 " << feedCol(0)
670 : << " DATA_DESC_ID " << msIter.dataDescriptionId()
671 0 : << debugpost;
672 0 : ScalarColumn<Double> timeCol(current, "TIME");
673 0 : ScalarColumn<Double> exposureCol(current, "EXPOSURE");
674 0 : ScalarColumn<Double> intervalCol(current, "INTERVAL");
675 0 : DataRealComponentAccessor dataCol(current);
676 0 : ArrayColumn<Bool> flagCol(current, "FLAG");
677 0 : ScalarColumn<Bool> flagRowCol(current, "FLAG_ROW");
678 0 : Vector<Double> timeList = timeCol.getColumn();
679 0 : Vector<Double> exposure = exposureCol.getColumn();
680 0 : Vector<Double> interval = intervalCol.getColumn();
681 0 : Vector<Double> timeInterval(timeList.nelements());
682 0 : Slice slice1(0, nrow - 1);
683 0 : Slice slice2(1, nrow - 1);
684 0 : timeInterval(slice1) = timeList(slice2) - timeList(slice1);
685 0 : timeInterval[nrow-1] = DBL_MAX;
686 0 : IPosition cellShape = flagCol.shape(0);
687 0 : size_t nelem = cellShape.product();
688 0 : Matrix<Float> dataSum(cellShape, new Float[nelem], TAKE_OVER);
689 0 : Matrix<Float> weightSum(cellShape, new Float[nelem], TAKE_OVER);
690 0 : dataSum = 0.0f;
691 0 : weightSum = 0.0f;
692 0 : Matrix<Bool> resultMask(cellShape, new Bool[nelem], TAKE_OVER);
693 0 : resultMask = true;
694 0 : Vector<Bool> flagRow = flagRowCol.getColumn();
695 0 : Double threshold = 1.1;
696 0 : Double timeCentroid = 0.0;
697 0 : size_t numSpectra = 0;
698 0 : Double effectiveExposure = 0.0;
699 0 : for (uInt i = 0; i < nrow; ++i) {
700 0 : if (flagRow(i)) {
701 0 : continue;
702 : }
703 :
704 0 : numSpectra++;
705 0 : timeCentroid += timeList[i];
706 0 : effectiveExposure += exposure[i];
707 :
708 0 : Matrix<Bool> mask = !flagCol(i);
709 0 : MaskedArray<Float> mdata(dataCol(i), mask);
710 0 : MaskedArray<Float> weight(Matrix<Float>(mdata.shape(), exposure[i]), mask);
711 0 : dataSum += mdata * weight;
712 0 : weightSum += weight;
713 :
714 0 : Double gap = 2.0 * timeInterval[i] /
715 0 : (interval[i] + interval[(i < nrow-1)?i+1:i]);
716 0 : if (gap > threshold) {
717 0 : debuglog << "flush accumulated data at row " << i << debugpost;
718 : // Here we can safely use data() since internal storeage
719 : // is contiguous
720 0 : Float *data_ = dataSum.data();
721 0 : const Float *weight_ = weightSum.data();
722 0 : Bool *flag_ = resultMask.data();
723 0 : for (size_t j = 0; j < dataSum.nelements(); ++j) {
724 0 : if (weight_[j] == 0.0f) {
725 0 : data_[j] = 0.0;
726 0 : flag_[j] = false;
727 : }
728 : else {
729 0 : data_[j] /= weight_[j];
730 : }
731 : }
732 :
733 0 : currSpw() = ispw;
734 0 : currField() = ifield;
735 0 : timeCentroid /= (Double)numSpectra;
736 0 : refTime() = timeCentroid;
737 0 : interval_ = effectiveExposure;
738 :
739 0 : debuglog << "spw " << ispw << ": solveAllRPar.shape=" << solveAllRPar().shape() << " nPar=" << nPar() << " nChanPar=" << nChanPar() << " nElem=" << nElem() << debugpost;
740 :
741 0 : size_t const nCorr = dataSum.shape()[0];
742 0 : Cube<Float> const rpar = dataSum.addDegenerate(1);
743 0 : Cube<Bool> const parOK = resultMask.addDegenerate(1);
744 0 : for (size_t iCorr = 0; iCorr < nCorr; ++iCorr) {
745 0 : solveAllRPar().yzPlane(iCorr) = rpar.yzPlane(iCorr);
746 0 : solveAllParOK().yzPlane(iCorr) = parOK.yzPlane(iCorr);
747 0 : solveAllParErr().yzPlane(iCorr) = 0.1; // TODO: this is tentative
748 0 : solveAllParSNR().yzPlane(iCorr) = 1.0; // TODO: this is tentative
749 : }
750 :
751 0 : keepNCT();
752 :
753 0 : dataSum = 0.0f;
754 0 : weightSum = 0.0f;
755 0 : resultMask = true;
756 0 : timeCentroid = 0.0;
757 0 : numSpectra = 0;
758 0 : effectiveExposure = 0.0;
759 0 : }
760 : }
761 : }
762 0 : }
763 :
764 0 : void SingleDishSkyCal::keepNCT() {
765 : // Call parent to do general stuff
766 : // This adds nElem() rows
767 0 : SolvableVisCal::keepNCT();
768 :
769 : // We are adding to the most-recently added rows
770 0 : RefRows rows(ct_->nrow()-nElem(),ct_->nrow()-1,1);
771 0 : Vector<Int> ant(nElem(), currAnt_);
772 :
773 : // update ANTENNA1 and ANTENNA2 with appropriate value
774 0 : CTMainColumns ncmc(*ct_);
775 0 : ncmc.antenna1().putColumnCells(rows,ant);
776 0 : ncmc.antenna2().putColumnCells(rows,ant);
777 :
778 : // update INTERVAL
779 0 : ncmc.interval().putColumnCells(rows,interval_);
780 0 : }
781 :
782 0 : void SingleDishSkyCal::initSolvePar()
783 : {
784 0 : debuglog << "SingleDishSkyCal::initSolvePar()" << debugpost;
785 0 : for (Int ispw=0;ispw<nSpw();++ispw) {
786 :
787 0 : currSpw()=ispw;
788 :
789 0 : switch(parType()) {
790 0 : case VisCalEnum::REAL: {
791 0 : solveAllRPar().resize(nPar(),nChanPar(),nElem());
792 0 : solveAllRPar()=0.0;
793 0 : solveRPar().reference(solveAllRPar());
794 0 : break;
795 : }
796 0 : default: {
797 0 : throw(AipsError("Internal error(Calibrater Module): Unsupported parameter type "
798 0 : "found in SingleDishSkyCal::initSolvePar()"));
799 : break;
800 : }
801 : }//switch
802 :
803 0 : solveAllParOK().resize(solveAllRPar().shape());
804 0 : solveAllParErr().resize(solveAllRPar().shape());
805 0 : solveAllParSNR().resize(solveAllRPar().shape());
806 0 : solveAllParOK()=false;
807 0 : solveAllParErr()=0.0;
808 0 : solveAllParSNR()=0.0;
809 0 : solveParOK().reference(solveAllParOK());
810 0 : solveParErr().reference(solveAllParErr());
811 0 : solveParSNR().reference(solveAllParSNR());
812 : }
813 0 : currSpw()=0;
814 0 : currAnt_ = 0;
815 :
816 0 : interval_.resize(nElem());
817 0 : }
818 :
819 0 : void SingleDishSkyCal::syncMeta2(const vi::VisBuffer2& vb)
820 : {
821 : // call method in parent class
822 0 : VisCal::syncMeta2(vb);
823 :
824 : // fill interval array with exposure
825 0 : interval_.reference(vb.exposure());
826 0 : debuglog << "SingleDishSkyCal::syncMeta2 interval_= " << interval_ << debugpost;
827 :
828 0 : setNumberOfCorrelationsPerSpw(vb.ms(), nCorr_);
829 0 : debuglog << "nCorr_ = " << nCorr_ << debugpost;
830 0 : debuglog << "currSpw() = " << currSpw() << debugpost;
831 0 : debuglog << "nPar() = " << nPar() << debugpost;
832 0 : }
833 :
834 0 : void SingleDishSkyCal::syncCalMat(const Bool &/*doInv*/)
835 : {
836 0 : debuglog << "SingleDishSkyCal::syncCalMat" << debugpost;
837 0 : debuglog << "nAnt()=" << nAnt() << ", nElem()=" << nElem() << ", nBln()=" << nBln() << debugpost;
838 0 : debuglog << "Spw " << currSpw() << "nchanmat, ncalmat" << nChanMat() << ", " << nCalMat() << debugpost;
839 0 : debuglog << "nChanPar = " << nChanPar() << debugpost;
840 0 : currentSky().resize(2, nChanMat(), nCalMat());
841 0 : currentSky().unique();
842 0 : currentSkyOK().resize(currentSky().shape());
843 0 : currentSkyOK().unique();
844 0 : debuglog << "currentSkyOK.shape()=" << currentSkyOK().shape() << debugpost;
845 :
846 : // sky data from caltable
847 0 : debuglog << "currRPar().shape()=" << currRPar().shape() << debugpost;
848 0 : if (currRPar().shape().product() > 0) {
849 0 : debuglog << "currRPar() = " << currRPar().xzPlane(0) << debugpost;
850 : }
851 :
852 0 : convertArray(currentSky(), currRPar());
853 0 : currentSkyOK() = currParOK();
854 0 : debuglog << "currentTime() = " << setprecision(16) << currTime() << debugpost;
855 0 : debuglog << "currentSky() = " << currentSky().xzPlane(0) << debugpost;
856 0 : debuglog << "currentSkyOK() = " << currentSkyOK().xzPlane(0) << debugpost;
857 :
858 : // weight calibration
859 0 : if (calWt())
860 0 : syncWtScale();
861 :
862 0 : debuglog << "SingleDishSkyCal::syncCalMat DONE" << debugpost;
863 0 : }
864 :
865 0 : void SingleDishSkyCal::syncDiffMat()
866 : {
867 0 : debuglog << "SingleDishSkyCal::syncDiffMat()" << debugpost;
868 0 : }
869 :
870 0 : void SingleDishSkyCal::syncWtScale()
871 : {
872 0 : debuglog << "syncWtScale" << debugpost;
873 :
874 : // allocate necessary memory
875 0 : currWtScale().resize(currentSky().shape());
876 0 : currWtScale() = 1.0;
877 :
878 : // Calculate the weight scaling
879 0 : if (tInterpType() == "nearest") {
880 0 : calcWtScale<NearestWeightScalingScheme>();
881 : }
882 : else {
883 0 : calcWtScale<LinearWeightScalingScheme>();
884 : }
885 :
886 0 : debuglog << "syncWtScale DONE" << debugpost;
887 0 : }
888 :
889 : template<class ScalingScheme>
890 0 : void SingleDishSkyCal::calcWtScale()
891 : {
892 0 : debuglog << "calcWtScale<ScalingScheme>" << debugpost;
893 0 : auto const key = std::make_pair(currObs(), currSpw());
894 0 : debuglog << "for obs " << currObs() << " and spw " << currSpw() << debugpost;
895 0 : decltype(wtScaleData_)::iterator iter = wtScaleData_.find(key);
896 0 : map<Int, Matrix<Double> > wtMap;
897 0 : if (iter == wtScaleData_.end()) {
898 0 : debuglog << "construct weight scaling data for obs " << currObs() << " spw " << currSpw() << debugpost;
899 0 : debuglog << "number of antennas = " << nAnt() << debugpost;
900 0 : for (Int iAnt = 0; iAnt < nAnt(); ++iAnt) {
901 0 : CTInterface cti(*ct_);
902 0 : MSSelection mss;
903 0 : mss.setObservationExpr(String::toString(currObs()));
904 0 : mss.setSpwExpr(String::toString(currSpw()));
905 0 : mss.setAntennaExpr(String::toString(iAnt) + "&&&");
906 0 : TableExprNode ten = mss.toTableExprNode(&cti);
907 0 : NewCalTable temp;
908 : try {
909 0 : getSelectedTable(temp, *ct_, ten, "");
910 0 : } catch (AipsError x) {
911 0 : continue;
912 : }
913 0 : temp = temp.sort("TIME");
914 0 : wtMap.emplace(std::piecewise_construct,
915 : std::forward_as_tuple(iAnt),
916 0 : std::forward_as_tuple(2, temp.nrow()));
917 0 : Matrix<Double> &arr = wtMap.at(iAnt);
918 0 : debuglog << "ant " << iAnt << " arr.shape = " << arr.shape() << debugpost;
919 0 : ScalarColumn<Double> col(temp, "TIME");
920 0 : auto timeCol = arr.row(0);
921 0 : col.getColumn(timeCol, False);
922 0 : col.attach(temp, "INTERVAL");
923 0 : auto intervalCol = arr.row(1);
924 0 : debuglog << "timeCol.size() == " << timeCol.size() << " intervalCol.size() = " << intervalCol.size() << debugpost;
925 0 : col.getColumn(intervalCol, False);
926 : }
927 0 : wtScaleData_.insert(std::make_pair(key, wtMap));
928 : } else {
929 0 : wtMap = iter->second;
930 : }
931 :
932 : {
933 : // for debugging
934 0 : debuglog << "wtMap keys: ";
935 0 : for (decltype(wtMap)::iterator i = wtMap.begin(); i != wtMap.end(); ++i) {
936 0 : debuglog << i->first << " ";
937 : }
938 0 : debuglog << debugpost;
939 : }
940 :
941 0 : for (Int iAnt = 0; iAnt < nAnt(); ++iAnt) {
942 0 : decltype(wtMap)::iterator i = wtMap.find(iAnt);
943 0 : if (i == wtMap.end()) {
944 0 : continue;
945 : }
946 0 : auto const &mat = i->second;
947 0 : debuglog << "matrix shape " << mat.shape() << debugpost;
948 0 : Vector<Double> const &timeCol = mat.row(0);
949 0 : Vector<Double> const &intervalCol = mat.row(1);
950 0 : size_t nrow = timeCol.size();
951 0 : debuglog << "timeCol = " << timeCol << debugpost;
952 0 : debuglog << "intervalCol = " << intervalCol << debugpost;
953 0 : debuglog << "iAnt " << iAnt << " nrow=" << nrow << debugpost;
954 0 : if (currTime() < timeCol[0]) {
955 0 : debuglog << "Use nearest OFF weight (0)" << debugpost;
956 0 : currWtScale().xyPlane(iAnt) = ScalingScheme::SimpleScale(interval_[iAnt], intervalCol[0]);
957 : }
958 0 : else if (currTime() > timeCol[nrow-1]) {
959 0 : debuglog << "Use nearest OFF weight (" << nrow-1 << ")" << debugpost;
960 0 : currWtScale().xyPlane(iAnt) = ScalingScheme::SimpleScale(interval_[iAnt], intervalCol[nrow-1]);
961 : }
962 : else {
963 0 : debuglog << "Use interpolated OFF weight" << debugpost;
964 0 : for (size_t irow = 0; irow < nrow ; ++irow) {
965 0 : if (currTime() == timeCol[irow]) {
966 0 : currWtScale().xyPlane(iAnt) = ScalingScheme::SimpleScale(interval_[iAnt], intervalCol[irow]);
967 0 : break;
968 : }
969 0 : else if (currTime() < timeCol[irow]) {
970 0 : currWtScale().xyPlane(iAnt) = ScalingScheme::InterpolateScale(currTime(), interval_[iAnt],
971 : timeCol[irow-1], intervalCol[irow-1],
972 : timeCol[irow], intervalCol[irow]);
973 0 : break;
974 : }
975 : }
976 : }
977 : }
978 0 : debuglog << "currWtScale() = " << currWtScale().xzPlane(0) << debugpost;
979 :
980 0 : debuglog << "calcWtScale<ScalingScheme> DONE" << debugpost;
981 0 : }
982 :
983 0 : Float SingleDishSkyCal::calcPowerNorm(Array<Float>& /*amp*/, const Array<Bool>& /*ok*/)
984 : {
985 0 : return 0.0f;
986 : }
987 :
988 0 : void SingleDishSkyCal::applyCal(VisBuffer& /*vb*/, Cube<Complex>& /*Vout*/, Bool /*trial*/)
989 : {
990 0 : throw AipsError("Single dish calibration doesn't support applyCal. Please use applyCal2");
991 : }
992 :
993 0 : void SingleDishSkyCal::applyCal2(vi::VisBuffer2 &vb, Cube<Complex> &Vout, Cube<Float> &Wout,
994 : Bool trial)
995 : {
996 0 : debuglog << "SingleDishSkyCal::applycal2" << debugpost;
997 0 : debuglog << "nrow, nchan=" << vb.nRows() << "," << vb.nChannels() << debugpost;
998 0 : debuglog << "antenna1: " << vb.antenna1() << debugpost;
999 0 : debuglog << "antenna2: " << vb.antenna2() << debugpost;
1000 0 : debuglog << "spw: " << vb.spectralWindows() << debugpost;
1001 0 : debuglog << "field: " << vb.fieldId() << debugpost;
1002 :
1003 : // References to VB2's contents' _data_
1004 0 : Vector<Bool> flagRv(vb.flagRow());
1005 0 : Vector<Int> a1v(vb.antenna1());
1006 0 : Vector<Int> a2v(vb.antenna2());
1007 0 : Cube<Bool> flagCube(vb.flagCube());
1008 0 : Cube<Complex> visCube(Vout);
1009 0 : ArrayIterator<Float> wt(Wout,2);
1010 0 : Matrix<Float> wtmat;
1011 :
1012 : // Data info/indices
1013 : Int* dataChan;
1014 0 : Bool* flagR=&flagRv(0);
1015 0 : Int* a1=&a1v(0);
1016 0 : Int* a2=&a2v(0);
1017 :
1018 : // iterate rows
1019 0 : Int nRow=vb.nRows();
1020 0 : Int nChanDat=vb.nChannels();
1021 0 : Vector<Int> dataChanv(vb.getChannelNumbers(0)); // All rows have same chans
1022 : // cout << currSpw() << " startChan() = " << startChan() << " nChanMat() = " << nChanMat() << " nChanDat="<<nChanDat <<endl;
1023 :
1024 : // setup calibration engine
1025 0 : engineC().setNumChannel(nChanDat);
1026 0 : engineC().setNumPolarization(vb.nCorrelations());
1027 :
1028 0 : debuglog << "typesize=" << engineC().typesize() << debugpost;
1029 :
1030 : // Matrix slice of visCube
1031 : // TODO: storage must be aligned for future use
1032 0 : Matrix<Complex> visCubeSlice;
1033 0 : Matrix<Bool> flagCubeSlice;
1034 :
1035 0 : for (Int row=0; row<nRow; row++,flagR++,a1++,a2++) {
1036 0 : debuglog << "spw: " << currSpw() << " antenna: " << *a1 << debugpost;
1037 0 : assert(*a1 == *a2);
1038 :
1039 : // Solution channel registration
1040 0 : Int solCh0(0);
1041 0 : dataChan=&dataChanv(0);
1042 :
1043 : // If cal _parameters_ are not freqDep (e.g., a delay)
1044 : // the startChan() should be the same as the first data channel
1045 0 : if (freqDepMat() && !freqDepPar())
1046 0 : startChan()=(*dataChan);
1047 :
1048 : // Solution and data array registration
1049 0 : engineC().sync(currentSky()(0,solCh0,*a1), currentSkyOK()(0,solCh0,*a1));
1050 0 : visCubeSlice.reference(visCube.xyPlane(row));
1051 0 : flagCubeSlice.reference(flagCube.xyPlane(row));
1052 :
1053 0 : if (trial) {
1054 : // only update flag info
1055 0 : engineC().flag(flagCubeSlice);
1056 : }
1057 : else {
1058 : // apply calibration
1059 0 : engineC().apply(visCubeSlice, flagCubeSlice);
1060 : }
1061 :
1062 : // If requested, update the weights
1063 0 : if (!trial && calWt()) {
1064 0 : wtmat.reference(wt.array());
1065 0 : updateWt2(wtmat,*a1);
1066 : }
1067 :
1068 0 : if (!trial)
1069 0 : wt.next();
1070 :
1071 : }
1072 0 : }
1073 :
1074 0 : void SingleDishSkyCal::selfGatherAndSolve(VisSet& vs, VisEquation& /*ve*/)
1075 : {
1076 0 : debuglog << "SingleDishSkyCal::selfGatherAndSolve()" << debugpost;
1077 :
1078 0 : MeasurementSet const &user_selection = vs.iter().ms();
1079 :
1080 0 : debuglog << "nspw = " << nSpw() << debugpost;
1081 : // Get and store the number of channels per spectral window
1082 0 : fillNChanParList(MeasurementSet(msName()), nChanParList());
1083 :
1084 : // Get and store the number of correlations per spectral window
1085 0 : setNumberOfCorrelationsPerSpw(user_selection, nCorr_);
1086 0 : debuglog << "nCorr_ = " << nCorr_ << debugpost;
1087 :
1088 : // Create a new in-memory calibration table to fill up
1089 0 : createMemCalTable();
1090 :
1091 : // Setup shape of solveAllRPar
1092 0 : nElem() = 1;
1093 0 : initSolvePar();
1094 :
1095 : // Select from user selection reference data associated with science target
1096 0 : MeasurementSet reference_data = selectReferenceData(user_selection);
1097 0 : logSink().origin(LogOrigin("SingleDishSkyCal","selfGatherAndSolve"));
1098 : {
1099 : // Log row numbers in a readable way
1100 0 : std::ostringstream msg;
1101 0 : auto * us_num_fmt = new CommaSeparatedThousands();
1102 0 : std::locale us_like_locale(std::locale::classic(), us_num_fmt);
1103 : // us_like_locale is responsible for deleting us_num_fmt from its destructor
1104 : // Ref: https://en.cppreference.com/w/cpp/locale/locale/locale
1105 : // Constructor (7) + Notes section
1106 : // So we must NOT delete us_num_fmt.
1107 0 : msg.imbue(us_like_locale);
1108 0 : msg << "Selected: " << std::right << std::setw(10) << reference_data.nrow() << " rows of reference data" << '\n'
1109 0 : << "out of : " << std::right << std::setw(10) << user_selection.nrow() << " rows of user-selected data";
1110 0 : logSink() << msg.str() << LogIO::POST;
1111 0 : }
1112 0 : logSink().origin(LogOrigin());
1113 0 : if (reference_data.nrow() == 0)
1114 0 : throw AipsError("No reference integration found in user-selected data. Please double-check your data selection criteria.");
1115 :
1116 : // Fill observing-mode-dependent columns of the calibration table
1117 : // Implementation is observing-mode-specific
1118 0 : fillCalibrationTable(reference_data);
1119 :
1120 : // Fill remaining columns of the calibration table,
1121 : // which are computed the same way for all observing modes
1122 : // ---- 1) FIELD_ID, SCAN_NUMBER, OBSERVATION_ID columns
1123 0 : assignCTScanField(*ct_, msName());
1124 :
1125 : // ---- 2) WEIGHT column
1126 : // update weight without Tsys
1127 : // formula is chanWidth [Hz] * interval [sec]
1128 0 : updateWeight(*ct_);
1129 :
1130 : // Store calibration table on disk
1131 0 : storeNCT();
1132 0 : }
1133 :
1134 0 : void SingleDishSkyCal::initializeSky()
1135 : {
1136 0 : debuglog << "SingleDishSkyCal::initializeSky()" << debugpost;
1137 0 : for (Int ispw=0;ispw<nSpw(); ispw++) {
1138 0 : currentSky_[ispw] = new Cube<Complex>();
1139 0 : currentSkyOK_[ispw] = new Cube<Bool>();
1140 0 : engineC_[ispw] = new SkyCal<Complex, Complex>();
1141 : }
1142 0 : }
1143 :
1144 0 : void SingleDishSkyCal::finalizeSky()
1145 : {
1146 0 : for (Int ispw=0;ispw<nSpw(); ispw++) {
1147 0 : if (currentSky_[ispw]) delete currentSky_[ispw];
1148 0 : if (currentSkyOK_[ispw]) delete currentSkyOK_[ispw];
1149 0 : if (engineC_[ispw]) delete engineC_[ispw];
1150 0 : if (engineF_[ispw]) delete engineF_[ispw];
1151 : }
1152 :
1153 0 : }
1154 :
1155 0 : void SingleDishSkyCal::initializeCorr()
1156 : {
1157 0 : File msPath(msName());
1158 0 : if (msPath.exists()) {
1159 0 : setNumberOfCorrelationsPerSpw(MeasurementSet(msName()), nCorr_);
1160 : } else {
1161 0 : nCorr_ = 1;
1162 : }
1163 0 : }
1164 :
1165 0 : void SingleDishSkyCal::updateWt2(Matrix<Float> &weight, const Int &antenna1)
1166 : {
1167 : // apply weight scaling factor
1168 0 : Matrix<Float> const factor = currWtScale().xyPlane(antenna1);
1169 0 : debuglog << "factor.shape() = " << factor.shape() << debugpost;
1170 0 : debuglog << "weight.shape() = " << weight.shape() << debugpost;
1171 0 : debuglog << "weight = " << weight << debugpost;
1172 :
1173 0 : auto const wtShape = weight.shape();
1174 0 : size_t const nCorr = wtShape[0];
1175 0 : size_t const nChan = wtShape[1];
1176 : // for each correlation
1177 0 : for (size_t iCorr = 0; iCorr < nCorr; ++iCorr) {
1178 0 : auto wSlice = weight.row(iCorr);
1179 0 : auto const fSlice = factor.row(iCorr);
1180 0 : if (fSlice.nelements() == nChan) {
1181 0 : wSlice *= fSlice;
1182 0 : } else if (nChan == 1) {
1183 : // take mean of spectral weight factor to apply it to scalar weight
1184 0 : wSlice *= mean(fSlice);
1185 : } else {
1186 0 : throw AipsError("Shape mismatch between input weight and weight scaling factor");
1187 : }
1188 0 : }
1189 0 : }
1190 :
1191 : //
1192 : // SingleDishPositionSwitchCal
1193 : //
1194 :
1195 : // Constructor
1196 0 : SingleDishPositionSwitchCal::SingleDishPositionSwitchCal(VisSet& vs)
1197 : : VisCal(vs),
1198 0 : SingleDishSkyCal(vs)
1199 : {
1200 0 : debuglog << "SingleDishPositionSwitchCal::SingleDishPositionSwitchCal(VisSet& vs)" << debugpost;
1201 0 : }
1202 :
1203 0 : SingleDishPositionSwitchCal::SingleDishPositionSwitchCal(const MSMetaInfoForCal& msmc)
1204 : : VisCal(msmc),
1205 0 : SingleDishSkyCal(msmc)
1206 : {
1207 0 : debuglog << "SingleDishPositionSwitchCal::SingleDishPositionSwitchCal(const MSMetaInfoForCal& msmc)" << debugpost;
1208 0 : }
1209 :
1210 0 : SingleDishPositionSwitchCal::SingleDishPositionSwitchCal(const Int& nAnt)
1211 : : VisCal(nAnt),
1212 0 : SingleDishSkyCal(nAnt)
1213 : {
1214 0 : debuglog << "SingleDishPositionSwitchCal::SingleDishPositionSwitchCal(const Int& nAnt)" << debugpost;
1215 0 : }
1216 :
1217 : // Destructor
1218 0 : SingleDishPositionSwitchCal::~SingleDishPositionSwitchCal()
1219 : {
1220 0 : debuglog << "SingleDishPositionSwitchCal::~SingleDishPositionSwitchCal()" << debugpost;
1221 0 : }
1222 :
1223 0 : MeasurementSet SingleDishPositionSwitchCal::selectReferenceData(MeasurementSet const &user_selection)
1224 : {
1225 0 : std::ostringstream qry;
1226 0 : constexpr auto eol = '\n';
1227 0 : qry << "using style python" << eol;
1228 : qry << "with [" << eol
1229 : << "select" << eol
1230 0 : << " [select TELESCOPE_NAME from ::OBSERVATION][OBSERVATION_ID] as TELESCOPE_NAME," << eol;
1231 : // Purposively not using TAQL's default mscal UDF library alias for derivedmscal
1232 : // to workaround a bug in casacore UDFBase::createUDF
1233 : qry << " derivedmscal.spwcol('NUM_CHAN') as NUM_CHAN" << eol
1234 : << "from" << eol
1235 : << " $1" << eol
1236 : << "] as metadata" << eol
1237 : << "select * from $1 , metadata" << eol
1238 0 : << "where " << eol;
1239 : // Row contains single-dish auto-correlation data,
1240 0 : qry << " ( ANTENNA1 == ANTENNA2 ) and" << eol ;
1241 0 : qry << " ( FEED1 == FEED2 ) and" << eol ;
1242 :
1243 : // has not been marked as invalid,
1244 0 : qry << " ( not(FLAG_ROW) ) and " << eol ;
1245 : // holds at least 1 data marked as valid,
1246 0 : qry << " ( not(all(FLAG)) ) and " << eol;
1247 : // ---- Note: both conditions above are required because FLAG and FLAG_ROW are not synchronized:
1248 : // a valid row (FLAG_ROW=False) may contain only invalid data: all(FLAG)=True
1249 :
1250 : // has observational intent: OBSERVE_TARGET#OFF_SOURCE,
1251 : qry << " ( STATE_ID in [ " << eol
1252 : << " select rowid() " << eol
1253 : << " from ::STATE" << eol
1254 : << " where " << eol
1255 : << " upper(OBS_MODE) ~ m/^OBSERVE_TARGET#OFF_SOURCE/ " << eol
1256 : << " ]" << eol
1257 0 : << " ) and" << eol;
1258 : // excluding - for ALMA - rows from Water Vapor Radiometers spectral windows, which have 4 channels
1259 : qry << " (" << eol
1260 : << " ( metadata.TELESCOPE_NAME != 'ALMA' ) or" << eol
1261 : << " (" << eol
1262 : << " ( metadata.TELESCOPE_NAME == 'ALMA' ) and" << eol
1263 : << " ( metadata.NUM_CHAN != 4 )" << eol
1264 : << " )" << eol
1265 0 : << " )" << eol;
1266 : debuglog << "SingleDishPositionSwitchCal::selectReferenceData(): taql query:" << eol
1267 0 : << qry.str() << debugpost;
1268 0 : return MeasurementSet(tableCommand(qry.str(), user_selection).table());
1269 0 : }
1270 :
1271 0 : void SingleDishPositionSwitchCal::fillCalibrationTable(casacore::MeasurementSet const &reference_data){
1272 0 : String dataColName = (reference_data.tableDesc().isColumn("FLOAT_DATA")) ? "FLOAT_DATA" : "DATA";
1273 :
1274 0 : if ( dataColName == "FLOAT_DATA")
1275 0 : fillCalibrationTable<FloatDataColumnAccessor>(reference_data);
1276 : else
1277 0 : fillCalibrationTable<DataColumnAccessor>(reference_data);
1278 0 : }
1279 :
1280 : template<class DataRealComponentAccessor>
1281 0 : void SingleDishPositionSwitchCal::fillCalibrationTable(casacore::MeasurementSet const &reference_data)
1282 : {
1283 0 : debuglog << "SingleDishPositionSwitchCal::fillCalibrationTable()" << debugpost;
1284 :
1285 : // Sort columns define the granularity at which we average data obtained
1286 : // by observing the reference position associated with the science target
1287 0 : constexpr size_t nSortColumns = 8;
1288 0 : Int columns[nSortColumns] = {
1289 : MS::OBSERVATION_ID,
1290 : MS::PROCESSOR_ID,
1291 : MS::FIELD_ID,
1292 : MS::ANTENNA1,
1293 : MS::FEED1,
1294 : MS::DATA_DESC_ID,
1295 : MS::SCAN_NUMBER,
1296 : MS::STATE_ID
1297 : };
1298 0 : Int *columnsPtr = columns;
1299 0 : Block<Int> sortColumns(nSortColumns, columnsPtr, false);
1300 :
1301 : // Iterator
1302 0 : constexpr Bool doGroupAllTimesTogether = true;
1303 0 : constexpr Double groupAllTimesTogether = doGroupAllTimesTogether ? 0.0 : -1.0;
1304 :
1305 0 : constexpr Bool doAddDefaultSortColumns = false;
1306 0 : constexpr Bool doStoreSortedTableOnDisk = false;
1307 :
1308 0 : MSIter msIter(reference_data, sortColumns,
1309 : groupAllTimesTogether, doAddDefaultSortColumns, doStoreSortedTableOnDisk);
1310 :
1311 : // Main loop: compute values of calibration table's columns
1312 0 : for (msIter.origin(); msIter.more(); msIter++) {
1313 0 : const auto iterTable = msIter.table();
1314 0 : const auto iterRows = iterTable.nrow();
1315 0 : if (iterRows == 0) continue;
1316 :
1317 : // TIME column of calibration table: mean of selected MAIN.TIME
1318 0 : ScalarColumn<Double> timeCol(iterTable, "TIME");
1319 0 : refTime() = mean(timeCol.getColumn());
1320 :
1321 : // FIELD_ID column of calibration table
1322 0 : currSpw() = msIter.spectralWindowId();
1323 :
1324 : // SPECTRAL_WINDOW_ID column of calibration table
1325 0 : currField() = msIter.fieldId();
1326 :
1327 : // ANTENNA1, ANTENNA2 columns of calibration table
1328 0 : ScalarColumn<Int> antenna1Col(iterTable, "ANTENNA1");
1329 0 : currAnt_ = antenna1Col(0);
1330 :
1331 : // INTERVAL column of calibration table: sum of selected MAIN.EXPOSURE
1332 0 : ScalarColumn<Double> exposureCol(iterTable, "EXPOSURE");
1333 0 : const auto exposure = exposureCol.getColumn();
1334 0 : interval_ = sum(exposure);
1335 :
1336 : // SCAN_NUMBER, OBSERVATION_ID columns of calibration table
1337 : // Not computed/updated here
1338 : #ifdef SDCALSKY_DEBUG
1339 : ScalarColumn<Int> scanNumberCol(iterTable, "SCAN_NUMBER");
1340 : const auto scan_number = scanNumberCol(0);
1341 : ScalarColumn<Int> stateIdCol(iterTable, "STATE_ID");
1342 : const auto state_id = stateIdCol(0);
1343 : cout << "field=" << currField()
1344 : << " ant=" << currAnt_
1345 : << " ddid=" << msIter.dataDescriptionId()
1346 : << " spw=" << currSpw()
1347 : << " scan_number=" << scan_number
1348 : << " state_id=" << state_id
1349 : << " nrows=" << iterRows
1350 : << endl;
1351 : #endif
1352 :
1353 : // FPARAM column of calibration table: weighted mean of valid data, weight = exposure
1354 : // + PARAMERR, FLAG, SNR
1355 : // ---- Get data shape from FLAG column
1356 0 : ArrayColumn<Bool> flagCol(iterTable, "FLAG");
1357 0 : const auto dataShape = flagCol.shape(0);
1358 0 : const auto nCorr = dataShape[0];
1359 0 : const auto nChannels = dataShape[1];
1360 : // ---- Initialize accumulators
1361 0 : Matrix<Float> weightedMean(nCorr, nChannels, 0.0f);
1362 0 : Matrix<Float> weightsSums(nCorr, nChannels, 0.0f);
1363 : // ---- Compute weighted mean of valid data
1364 0 : DataRealComponentAccessor dataAccessor(iterTable);
1365 0 : for (std::remove_const<decltype(iterRows)>::type iterRow = 0; iterRow < iterRows ; iterRow++){
1366 0 : Matrix<Bool> dataIsValid = not flagCol(iterRow);
1367 0 : MaskedArray<Float> validData(dataAccessor(iterRow), dataIsValid);
1368 0 : const auto rowExposure = static_cast<Float>(exposure[iterRow]);
1369 0 : weightedMean += validData * rowExposure;
1370 0 : MaskedArray<Float> validWeight(Matrix<Float>(validData.shape(), rowExposure), dataIsValid);
1371 0 : weightsSums += validWeight;
1372 : }
1373 0 : const Matrix<Bool> weightsSumsIsNonZero = ( weightsSums != 0.0f );
1374 0 : weightedMean /= MaskedArray<Float>(weightsSums,weightsSumsIsNonZero);
1375 : // ---- Update solveAll*() members
1376 0 : const Cube<Float> realParam = weightedMean.addDegenerate(1);
1377 0 : const Cube<Bool> realParamIsValid = weightsSumsIsNonZero.addDegenerate(1);
1378 0 : for (std::remove_const<decltype(nCorr)>::type corr = 0; corr < nCorr; corr++) {
1379 0 : solveAllRPar().yzPlane(corr) = realParam.yzPlane(corr); // FPARAM
1380 0 : solveAllParOK().yzPlane(corr) = realParamIsValid.yzPlane(corr); // not FLAG
1381 0 : solveAllParErr().yzPlane(corr) = 0.1; // PARAMERR. TODO: this is tentative
1382 0 : solveAllParSNR().yzPlane(corr) = 1.0; // SNR. TODO: this is tentative
1383 : }
1384 :
1385 : // WEIGHT column of calibration table
1386 : //
1387 :
1388 : // Update in-memory calibration table
1389 0 : keepNCT();
1390 : }
1391 0 : }
1392 :
1393 : //
1394 : // SingleDishRasterCal
1395 : //
1396 :
1397 : // Constructor
1398 0 : SingleDishRasterCal::SingleDishRasterCal(VisSet& vs)
1399 : : VisCal(vs),
1400 : SingleDishSkyCal(vs),
1401 0 : fraction_(0.1),
1402 0 : numEdge_(-1)
1403 : {
1404 0 : debuglog << "SingleDishRasterCal::SingleDishRasterCal(VisSet& vs)" << debugpost;
1405 0 : }
1406 :
1407 0 : SingleDishRasterCal::SingleDishRasterCal(const MSMetaInfoForCal& msmc)
1408 : : VisCal(msmc),
1409 : SingleDishSkyCal(msmc),
1410 0 : fraction_(0.1),
1411 0 : numEdge_(-1)
1412 : {
1413 0 : debuglog << "SingleDishRasterCal::SingleDishRasterCal(const MSMetaInfoForCal& msmc)" << debugpost;
1414 0 : }
1415 :
1416 0 : SingleDishRasterCal::SingleDishRasterCal(const Int& nAnt)
1417 : : VisCal(nAnt),
1418 0 : SingleDishSkyCal(nAnt)
1419 : {
1420 0 : debuglog << "SingleDishRasterCal::SingleDishRasterCal(const Int& nAnt)" << debugpost;
1421 0 : }
1422 :
1423 : // Destructor
1424 0 : SingleDishRasterCal::~SingleDishRasterCal()
1425 : {
1426 0 : debuglog << "SingleDishRasterCal::~SingleDishRasterCal()" << debugpost;
1427 0 : }
1428 :
1429 0 : void SingleDishRasterCal::setSolve(const Record& solve)
1430 : {
1431 : // edge detection parameter for otfraster mode
1432 0 : if (solve.isDefined("fraction")) {
1433 0 : fraction_ = solve.asFloat("fraction");
1434 : }
1435 0 : if (solve.isDefined("numedge")) {
1436 0 : numEdge_ = solve.asInt("numedge");
1437 : }
1438 :
1439 0 : logSink() << "fraction=" << fraction_ << endl
1440 0 : << "numedge=" << numEdge_ << LogIO::POST;
1441 :
1442 : // call parent setSolve
1443 0 : SolvableVisCal::setSolve(solve);
1444 0 : }
1445 :
1446 0 : MeasurementSet SingleDishRasterCal::selectReferenceData(MeasurementSet const &ms)
1447 : {
1448 0 : debuglog << "SingleDishRasterCal::selectReferenceData" << debugpost;
1449 0 : const Record specify;
1450 0 : std::ostringstream oss;
1451 0 : oss << "SELECT FROM $1 WHERE ";
1452 0 : String delimiter = "";
1453 : // for (Int iant = 0; iant < nAnt(); ++iant) {
1454 : // Vector<String> timeRangeList = detectRaster(msName(), iant, fraction_, numEdge_);
1455 : // debuglog << "timeRangeList=" << ::toString(timeRangeList) << debugpost;
1456 : // oss << delimiter;
1457 : // oss << "(ANTENNA1 == " << iant << " && ANTENNA2 == " << iant << " && (";
1458 : // String separator = "";
1459 : // for (size_t i = 0; i < timeRangeList.size(); ++i) {
1460 : // if (timeRangeList[i].size() > 0) {
1461 : // oss << separator << "(" << timeRangeList[i] << ")";
1462 : // separator = " || ";
1463 : // }
1464 : // }
1465 : // oss << "))";
1466 : // debuglog << "oss.str()=" << oss.str() << debugpost;
1467 : // delimiter = " || ";
1468 : // }
1469 : // use ANTENNA 0 for reference antenna
1470 0 : Vector<String> timeRangeList = detectRaster(ms, 0, fraction_, numEdge_);
1471 0 : debuglog << "timeRangeList=" << ::toString(timeRangeList) << debugpost;
1472 0 : oss << delimiter;
1473 0 : oss << "(ANTENNA1 == ANTENNA2 && (";
1474 0 : String separator = "";
1475 0 : for (size_t i = 0; i < timeRangeList.size(); ++i) {
1476 0 : if (timeRangeList[i].size() > 0) {
1477 0 : oss << separator << "(" << timeRangeList[i] << ")";
1478 0 : separator = " || ";
1479 : }
1480 : }
1481 0 : oss << "))";
1482 0 : debuglog << "oss.str()=" << oss.str() << debugpost;
1483 :
1484 : oss //<< ")"
1485 0 : << " ORDER BY FIELD_ID, ANTENNA1, FEED1, DATA_DESC_ID, TIME";
1486 0 : return MeasurementSet(tableCommand(oss.str(), ms).table());
1487 0 : }
1488 :
1489 : //
1490 : // SingleDishOtfCal
1491 : //
1492 :
1493 : // Constructor
1494 0 : SingleDishOtfCal::SingleDishOtfCal(VisSet& vs)
1495 : : VisCal(vs),
1496 : SingleDishSkyCal(vs),
1497 0 : fraction_(0.1),
1498 0 : pixel_scale_(0.5),
1499 0 : msSel_(vs.iter().ms())
1500 : {
1501 0 : debuglog << "SingleDishOtfCal::SingleDishOtfCal(VisSet& vs)" << debugpost;
1502 0 : }
1503 :
1504 : /*
1505 : SingleDishOtfCal::SingleDishOtfCal(const MSMetaInfoForCal& msmc)
1506 : : VisCal(msmc),
1507 : SingleDishSkyCal(msmc),
1508 : fraction_(0.1),
1509 : pixel_scale_(0.5),
1510 : msSel_(vs.iter().ms()) ************need MS!
1511 : {
1512 : debuglog << "SingleDishOtfCal::SingleDishOtfCal(VisSet& vs)" << debugpost;
1513 : }
1514 : */
1515 0 : void SingleDishOtfCal::setSolve(const Record& solve)
1516 : {
1517 : // edge detection parameter for otfraster mode
1518 0 : if (solve.isDefined("fraction")) {
1519 0 : fraction_ = solve.asFloat("fraction");
1520 : }
1521 :
1522 0 : logSink() << "fraction=" << fraction_ << LogIO::POST;
1523 :
1524 : // call parent setSolve
1525 0 : SolvableVisCal::setSolve(solve);
1526 0 : }
1527 :
1528 : /*
1529 : SingleDishOtfCal::SingleDishOtfCal(const Int& nAnt)
1530 : : VisCal(nAnt),
1531 : SingleDishSkyCal(nAnt)
1532 : {
1533 : debuglog << "SingleDishOtfCal::SingleDishOtfCal(const Int& nAnt)" << debugpost;
1534 : }
1535 : */
1536 :
1537 : // Destructor
1538 0 : SingleDishOtfCal::~SingleDishOtfCal()
1539 : {
1540 0 : debuglog << "SingleDishOtfCal::~SingleDishOtfCal()" << debugpost;
1541 0 : }
1542 :
1543 0 : MeasurementSet SingleDishOtfCal::selectReferenceData(MeasurementSet const &ms)
1544 : {
1545 0 : PointingDirectionCalculator calc(ms);
1546 0 : calc.setDirectionListMatrixShape(PointingDirectionCalculator::ROW_MAJOR);
1547 :
1548 : // Check the coordinates system type used to store the pointing measurements
1549 0 : const MSPointing& tbl_pointing = ms.pointing();
1550 0 : MSPointingColumns pointing_cols(tbl_pointing);
1551 0 : const ROArrayMeasColumn< MDirection >& direction_cols = pointing_cols.directionMeasCol();
1552 0 : const MeasRef<MDirection>& direction_ref_frame = direction_cols.getMeasRef();
1553 0 : uInt ref_frame_type = direction_ref_frame.getType();
1554 :
1555 : // If non-celestial coordinates (AZEL*) are used, convert to celestial ones
1556 0 : switch (ref_frame_type) {
1557 0 : case MDirection::AZEL : // Fall through
1558 : case MDirection::AZELSW :
1559 : case MDirection::AZELGEO :
1560 : case MDirection::AZELSWGEO : {
1561 0 : const String& ref_frame_name = MDirection::showType(ref_frame_type);
1562 0 : debuglog << "Reference frame of pointings coordinates is non-celestial: " << ref_frame_name << debugpost;
1563 0 : String j2000(MDirection::showType(MDirection::J2000));
1564 0 : debuglog << "Pointings coordinates will be converted to: " << j2000 << debugpost;
1565 0 : calc.setFrame(j2000);
1566 0 : }
1567 : }
1568 : // Extract edge pointings for each (field_id,antenna,spectral window) triple
1569 : // MeasurementSet 2 specification / FIELD table:
1570 : // . FIELD_ID column is removed
1571 : // . FIELD table is directly indexed using the FIELD_ID value in MAIN
1572 0 : const MSField& tbl_field = ms.field();
1573 0 : const String &field_col_name_str = tbl_field.columnName(MSField::MSFieldEnums::SOURCE_ID);
1574 0 : ScalarColumn<Int> source_id_col(tbl_field, field_col_name_str);
1575 0 : const MSAntenna& tbl_antenna = ms.antenna();
1576 0 : const String &col_name_str = tbl_antenna.columnName(MSAntenna::MSAntennaEnums::NAME);
1577 0 : ScalarColumn<String> antenna_name(tbl_antenna,col_name_str);
1578 0 : const MSSpectralWindow& tbl_spectral_window = ms.spectralWindow();
1579 :
1580 : // make a map between SOURCE_ID and source NAME
1581 0 : const MSSource &tbl_source = ms.source();
1582 0 : ScalarColumn<Int> id_col(tbl_source, tbl_source.columnName(MSSource::MSSourceEnums::SOURCE_ID));
1583 0 : ScalarColumn<String> name_col(tbl_source, tbl_source.columnName(MSSource::MSSourceEnums::NAME));
1584 0 : std::map<Int, String> source_map;
1585 0 : for (uInt irow = 0; irow < tbl_source.nrow(); ++irow) {
1586 0 : auto source_id = id_col(irow);
1587 0 : if (source_map.find(source_id) == source_map.end()) {
1588 0 : source_map[source_id] = name_col(irow);
1589 : }
1590 : }
1591 :
1592 0 : Vector<casacore::rownr_t> rowList;
1593 :
1594 0 : for (uInt field_id=0; field_id < tbl_field.nrow(); ++field_id){
1595 0 : String field_sel(casacore::String::toString<uInt>(field_id));
1596 0 : String const source_name = source_map.at(source_id_col(field_id));
1597 :
1598 : // Set ephemeris flag if source name is the one recognized as a moving source
1599 0 : if (isEphemeris(source_name)) {
1600 0 : calc.setMovingSource(source_name);
1601 : }
1602 : else {
1603 0 : calc.unsetMovingSource();
1604 : }
1605 :
1606 0 : for (uInt ant_id=0; ant_id < tbl_antenna.nrow(); ++ant_id){
1607 0 : String ant_sel(antenna_name(ant_id) + "&&&");
1608 0 : for (uInt spw_id=0; spw_id < tbl_spectral_window.nrow(); ++spw_id){
1609 0 : String spw_sel(casacore::String::toString<uInt>(spw_id));
1610 : // Filter user selection by (field_id,antenna,spectral window) triple
1611 : try {
1612 0 : calc.selectData(ant_sel,spw_sel,field_sel);
1613 : }
1614 0 : catch (AipsError& e) { // Empty selection
1615 : // Note: when the underlying MSSelection is empty
1616 : // MSSelection internally catches an MSSelectionError error
1617 : // but does not re-throw it. It throws instead an AipsError
1618 : // copy-constructed from the MSSelectionError
1619 0 : continue;
1620 0 : }
1621 : debuglog << "field_id: " << field_id
1622 : << " ant_id: " << ant_id
1623 : << " spw: " << spw_id
1624 0 : << " ==> selection rows: " << calc.getNrowForSelectedMS() << debugpost;
1625 : // Get time-interpolated celestial pointing directions for the filtered user selection
1626 0 : Matrix<Double> pointings_dirs = calc.getDirection();
1627 : // Project directions onto image plane
1628 : // pixel_scale_ :
1629 : // . hard-coded to 0.5 in constructor
1630 : // . is applied to the median separation of consecutive pointing directions by the projector
1631 : // . projector pixel size = 0.5*directions_median
1632 0 : debuglog << "pixel_scale:" << pixel_scale_ << debugpost;
1633 0 : OrthographicProjector p(pixel_scale_);
1634 0 : p.setDirection(pointings_dirs);
1635 0 : const Matrix<Double> &pointings_coords = p.project();
1636 : // Extract edges of the observed region for the (field_id,antenna,spectral window) triple
1637 0 : SakuraAlignedArray<Double> pointings_x(pointings_coords.row(0));
1638 0 : SakuraAlignedArray<Double> pointings_y(pointings_coords.row(1));
1639 0 : SakuraAlignedArray<Bool> is_edge_storage(pointings_coords.ncolumn());
1640 0 : Vector<Bool> is_edge = is_edge_storage.casaVector();
1641 0 : is_edge = false;
1642 0 : double pixel_size = 0.0;
1643 : {
1644 : // CAS-9956
1645 : // Mitigation of memory usage due to unexpectedly large number of pixels.
1646 : // Currently CreateMaskNearEdgeDouble requires 2*sizeof(size_t)*num_pixels bytes
1647 : // of memory. If this value exceeds 2GB, the mitigation will be activated.
1648 0 : Double const num_pixels = p.p_size()[0] * p.p_size()[1];
1649 0 : auto const estimated_memory = num_pixels * 2 * sizeof(size_t);
1650 0 : constexpr Double kMaxMemory = 2.0e9;
1651 0 : if (estimated_memory >= kMaxMemory && pixel_scale_ < 1.0) {
1652 0 : LogIO os;
1653 0 : os << LogOrigin("PointingDirectionProjector", "scale_and_center", WHERE);
1654 0 : os << LogIO::WARN << "Estimated Memory: " << estimated_memory << LogIO::POST;
1655 0 : os << LogIO::WARN << "Mitigation of memory usage is activated. pixel scale is set to 1.0" << LogIO::POST;
1656 : // pixel_size can be set to 2.0 since projection grid spacing is estimated from half of median separation
1657 : // between neighboring pixels so that pixel_width will become about 1.0 if pixel_size is 0.
1658 0 : pixel_size = 2.0;
1659 0 : os << LogIO::WARN << "pixel_size is set to " << pixel_size << LogIO::POST;
1660 0 : }
1661 : }
1662 : // libsakura 2.0: setting pixel_size=0.0 means that CreateMaskNearEdgeDouble will
1663 : // . compute the median separation of consecutive pointing coordinates
1664 : // . use an "edge detection pixel size" = 0.5*coordinates_median (pixel scale hard-coded to 0.5)
1665 0 : debuglog << "sakura library function call: parameters info:" << debugpost;
1666 0 : debuglog << "in: fraction: " << fraction_ << debugpost;
1667 0 : debuglog << "in: pixel size: " << pixel_size << debugpost;
1668 0 : debuglog << "in: pixels count: (nx = " << p.p_size()[0] << " , ny = " << p.p_size()[1] << debugpost;
1669 0 : debuglog << "in: pointings_coords.ncolumn(): " << pointings_coords.ncolumn() << debugpost;
1670 0 : LIBSAKURA_SYMBOL(Status) status = LIBSAKURA_SYMBOL(CreateMaskNearEdgeDouble)(
1671 : fraction_, pixel_size,
1672 0 : pointings_coords.ncolumn(), pointings_x.data(), pointings_y.data(),
1673 : nullptr /* blc_x */, nullptr /* blc_y */,
1674 : nullptr /* trc_x */, nullptr /* trc_y */,
1675 : is_edge.data());
1676 0 : bool edges_detection_ok = ( status == LIBSAKURA_SYMBOL(Status_kOK) );
1677 0 : if ( ! edges_detection_ok ) {
1678 0 : debuglog << "sakura error: status=" << status << debugpost;
1679 : }
1680 0 : AlwaysAssert(edges_detection_ok,AipsError);
1681 : // Compute ROW ids of detected edges. ROW "ids" are ROW ids in the MS filtered by user selection.
1682 0 : auto index_2_rowid = calc.getRowIdForOriginalMS();
1683 0 : size_t edges_count = ntrue(is_edge);
1684 0 : size_t rowListIndex = rowList.size();
1685 0 : rowList.resize(rowList.size() + edges_count, True);
1686 0 : for (size_t i = 0; i < is_edge.size(); ++i){
1687 0 : if ( is_edge[i] ) {
1688 0 : rowList[rowListIndex] = index_2_rowid[i]; // i;
1689 0 : ++rowListIndex;
1690 : }
1691 : }
1692 0 : debuglog << "edges_count=" << edges_count << debugpost;
1693 0 : AlwaysAssert(edges_count > 0, AipsError);
1694 : #ifdef SDCALSKY_DEBUG
1695 : stringstream fname;
1696 : fname << calTableName().c_str() << ".edges."
1697 : << field_id << "_" << ant_id << "_" << spw_id
1698 : << ".csv" ;
1699 : debuglog << "Save pointing directions and coordinates to:" << debugpost;
1700 : debuglog << fname.str() << debugpost;
1701 : ofstream ofs(fname.str());
1702 : AlwaysAssert(ofs.good(), AipsError);
1703 : ofs << "row_id,field_id,ant_id,spw_id,triple_key,dir_0,dir_1,coord_0,coord_1,edge_0,edge_1,is_edge" << endl;
1704 : const auto &d0 = pointings_dirs.row(0);
1705 : const auto &d1 = pointings_dirs.row(1);
1706 : const auto &c0 = pointings_coords.row(0);
1707 : const auto &c1 = pointings_coords.row(1);
1708 : for (uInt j=0; j<d0.size(); j++) {
1709 : ofs << index_2_rowid[j] << ","
1710 : << field_id << "," << ant_id << "," << spw_id << ","
1711 : << field_id << "_" << ant_id << "_" << spw_id << ","
1712 : << d0(j) << "," << d1(j) << ","
1713 : << c0(j) << "," << c1(j) << "," ;
1714 : if ( is_edge[j] ) ofs << c0(j) << "," << c1(j) << "," << 1 << endl;
1715 : else ofs << ",," << 0 << endl;
1716 : }
1717 : #endif
1718 0 : }
1719 0 : }
1720 0 : }
1721 0 : Bool have_off_spectra = (rowList.size() > 0);
1722 0 : AlwaysAssert(have_off_spectra, AipsError);
1723 0 : MeasurementSet msSel = ms(rowList);
1724 0 : debuglog << "rowList = " << rowList << debugpost;
1725 0 : return msSel;
1726 0 : }
1727 :
1728 : } //# NAMESPACE CASA - END
1729 :
|
