Line data Source code
1 : #include <iostream>
2 : #include <set>
3 : #include <memory>
4 : #include <sstream>
5 :
6 : #include <mstransform/TVI/PointingInterpolationTVI.h>
7 :
8 : // Visibility Buffer
9 : #include <msvis/MSVis/VisBuffer2.h>
10 :
11 : // Visibility Iterator
12 : #include <msvis/MSVis/VisibilityIteratorImpl2.h>
13 : #include <msvis/MSVis/TransformingVi2.h>
14 :
15 : // Measures
16 : #include <casacore/measures/Measures/MDirection.h>
17 : #include <casacore/measures/Measures/MeasFrame.h>
18 : #include <casacore/tables/TaQL.h>
19 : #include <casacore/casa/Arrays/Cube.h>
20 : #include <casacore/casa/Exceptions/Error.h>
21 :
22 :
23 :
24 : // Messy
25 : // !!! This include is required and must be after PI_TVI.h
26 : // due to forward declarations
27 : // of the MeasurementSet class in files included by PI_TVI.h
28 : #include <casacore/ms/MeasurementSets.h>
29 : #include <casacore/ms/MeasurementSets/MSColumns.h>
30 :
31 : using namespace casa::vi;
32 : using namespace casacore;
33 : using std::cout;
34 : using std::endl;
35 :
36 : using PI_TVI = PointingInterpolationTVI;
37 : using PI_Vi2Factory = PointingInterpolationVi2Factory;
38 : using PI_TVILayerFactory = PointingInterpolationTVILayerFactory;
39 :
40 0 : PI_TVI::PointingInterpolationTVI(ViImplementation2 *inputVII) :
41 0 : TransformingVi2(inputVII) {
42 : // Initialize associated output VisBuffer2 object
43 0 : VisBufferOptions vbOptions {VbNoOptions};
44 0 : setVisBuffer(createAttachedVisBuffer(vbOptions));
45 :
46 : // Expecting a pointing table
47 0 : inputVII->ms().pointing().throwIfNull();
48 :
49 : // Input MS / pointing direction units and ref
50 0 : MSPointingColumns mspc(getVii()->ms().pointing());
51 0 : auto dirUnits = mspc.directionMeasCol().measDesc().getUnits();
52 0 : lonUnit_ = dirUnits[0];
53 0 : latUnit_ = dirUnits[1];
54 :
55 0 : MDirection::Ref ref(mspc.directionMeasCol().measDesc().getRefCode());
56 0 : dirRef_ = ref;
57 0 : toRef_.setType(MDirection::ICRS);
58 :
59 0 : setupInterpolator();
60 0 : }
61 :
62 0 : void PI_TVI::setupInterpolator(){
63 :
64 : try {
65 : // To clarify: behavior of TVIs iterating over multiple MSs = ?
66 0 : const auto & ms = this->getVii()->ms();
67 :
68 : // Total number of antennas
69 0 : auto nAnts = ms.antenna().nrow();
70 :
71 : // Selected, data-taking antennas
72 0 : Vector<Bool> isSelected(nAnts);
73 0 : MSColumns msCols(ms);
74 0 : cout << "Extracting active antennas ..." << endl;
75 0 : for (auto antId : msCols.antenna1().getColumn()) isSelected[antId] = true;
76 0 : for (auto antId : msCols.antenna2().getColumn()) isSelected[antId] = true;
77 0 : std::set<Int> selectedAntIds;
78 0 : for (uInt antId=0; antId<nAnts; antId++){
79 0 : if (isSelected[antId]) selectedAntIds.insert(antId);
80 : }
81 :
82 : //cout << "TVI input ms rows: " << ms.nrow() << endl;
83 : //cout << "TVI input ms antennas: " << nAnts << endl;
84 : //cout << "TVI input ms selected antennas: " << selectedAntIds.size() << endl;
85 :
86 : // Distinct pointings
87 0 : std::stringstream ss;
88 : ss << "select ANTENNA_ID, TIME from $1 "
89 0 : << "where ANTENNA_ID in " << taQLSet(selectedAntIds) << " "
90 0 : << "order by distinct ANTENNA_ID, TIME";
91 0 : String qryDistinctPointings(ss.str());
92 :
93 : //cout << "Query string:" << endl;
94 : //cout << qryDistinctPointings << endl;
95 :
96 0 : auto keys = tableCommand(qryDistinctPointings,ms.pointing()).table();
97 0 : auto distinctPointings = ms.pointing()(keys.rowNumbers());
98 :
99 : // Pointings per antenna
100 0 : auto distinctPointingsAntIds = ScalarColumn<Int>(distinctPointings,"ANTENNA_ID").getColumn();
101 0 : Vector<uInt> antsPointingsCount(nAnts,0);
102 0 : for ( auto antId : distinctPointingsAntIds ) ++antsPointingsCount[antId];
103 0 : Vector<Bool> hasPointings(nAnts,false);
104 0 : for (uInt antId=0; antId<nAnts; antId++){
105 0 : if (antsPointingsCount[antId] >0) hasPointings[antId] = true;
106 : }
107 :
108 : //for (uInt antId=0; antId < antsPointingsCount.nelements(); ++antId)
109 : // cout << "ant: " << antId << " pointings: " << antsPointingsCount[antId] << endl;
110 :
111 : // Extract pointings data from distinctPointings table
112 : using PointingTimes = Vector<Double>;
113 : using PointingDir = Vector<Double>;
114 : using PointingDirs = Vector<PointingDir>;
115 0 : Vector<PointingTimes> antsPointingTimes(nAnts);
116 0 : Vector<PointingDirs> antsPointingDirs(nAnts);
117 :
118 0 : auto distinctPointingTimes = ScalarColumn<Double>(distinctPointings,"TIME").getColumn();
119 0 : ArrayMeasColumn<MDirection> distinctPointingDirsArray(distinctPointings,"DIRECTION");
120 0 : const auto & dirArray = distinctPointingDirsArray;
121 :
122 0 : uInt antStartRow = 0; // In distinct pointings table
123 0 : for ( uInt antId=0; antId<nAnts; antId++ ) {
124 0 : if ( not isSelected[antId] or not hasPointings[antId] ) continue;
125 : // Time
126 0 : auto antPointingsCount = antsPointingsCount[antId];
127 0 : antsPointingTimes[antId] = distinctPointingTimes(Slice(antStartRow,antPointingsCount));
128 : // Direction
129 0 : auto & antPointingDirs = antsPointingDirs[antId];
130 0 : antPointingDirs.resize(antPointingsCount);
131 0 : for ( uInt pIndex = 0 ; pIndex < antPointingsCount ; pIndex++ ){
132 0 : auto pointingRow = antStartRow + pIndex;
133 0 : auto rowDirArray = dirArray(pointingRow);
134 0 : auto rowDirection = rowDirArray.begin()->getAngle().getValue();
135 0 : antPointingDirs[pIndex] = rowDirection;
136 : }
137 : // Next antenna in distinct pointings table
138 0 : antStartRow += antPointingsCount;
139 : }
140 :
141 0 : interpolator_.setData(antsPointingTimes,antsPointingDirs,isSelected);
142 :
143 : // cout << "myTVI.constructor: done" << endl;
144 0 : } catch (const std::exception& e) {
145 0 : cout << "Unexcepted exception caught in " << __func__ << endl;
146 0 : cout << "Error message: " << endl;
147 0 : cout << e.what() << endl;
148 0 : exit(1);
149 : }
150 0 : }
151 :
152 0 : PI_TVI::~PointingInterpolationTVI() {
153 : // cout << "myTVI: destructor: start" << endl;
154 : // delete interpolator_;
155 : // delete distinctPointingsCols_;
156 0 : inputVii_p = NULL; // Following TVI Development Guide
157 : //cout << "myTVI: destructor: done" << endl;
158 0 : }
159 :
160 0 : String PI_TVI::ViiType() const {
161 0 : return String("PointingInterpolation( ") + getVii()->ViiType() + " )";
162 : }
163 :
164 0 : void PI_TVI::origin(){
165 : // Drive underlying ViImplementation
166 0 : getVii()->origin();
167 :
168 : // Synchronize own output VisBuffer
169 0 : configureNewSubchunk();
170 0 : }
171 :
172 0 : void PI_TVI::next(){
173 : // Drive underlying ViImplementation
174 0 : getVii()->next();
175 :
176 : // Synchronize own output VisBuffer
177 0 : configureNewSubchunk();
178 0 : }
179 :
180 0 : PI_TVI::Interpolator &PI_TVI::getInterpolator() {
181 0 : return interpolator_;
182 : }
183 :
184 0 : void PI_TVI::setOutputDirectionFrame(MDirection::Types toRefType){
185 0 : toRef_.setType(toRefType);
186 0 : }
187 :
188 0 : std::pair<bool, MDirection> PI_TVI::getPointingAngle (int antId, double timeStamp) const {
189 0 : auto dirValue = interpolator_.pointingDir(antId,timeStamp);
190 0 : Quantity qLon(dirValue[0], lonUnit_);
191 0 : Quantity qLat(dirValue[1], latUnit_);
192 0 : MDirection fromDir(qLon,qLat,dirRef_);
193 :
194 0 : if (dirRef_.getType() == toRef_.getType())
195 0 : return std::make_pair(true,fromDir);
196 :
197 : using InterpMethod = PI_TVI::Interpolator::InterpMethod;
198 0 : double refTimeStamp = interpolator_.getInterpMethod() == InterpMethod::NEAREST ?
199 0 : interpolator_.nearestPointingTimeStamp()
200 0 : : timeStamp;
201 0 : MEpoch refEpoch(Quantity(refTimeStamp,"s"),MEpoch::UTC);
202 :
203 0 : auto antPosition = getVisBuffer()->subtableColumns().antenna().positionMeas()(antId);
204 0 : auto toDir = MDirection::Convert(fromDir,
205 0 : MDirection::Ref(toRef_.getType(),MeasFrame(refEpoch,antPosition)))();
206 :
207 0 : return std::make_pair(true,toDir);
208 : }
209 :
210 0 : String PI_TVI::taQLSet(const std::set<int> & inputSet){
211 0 : std::stringstream ss;
212 0 : ss << "[ ";
213 0 : if (not inputSet.empty()){
214 0 : auto it = inputSet.begin();
215 0 : ss << *it;
216 0 : while ( ++it != inputSet.end() ) ss << "," << *it ;
217 : }
218 0 : ss << " ]";
219 0 : return String(ss.str());
220 : }
221 :
222 : using PI_Interp = PI_TVI::Interpolator;
223 :
224 0 : PI_Interp::Interpolator()
225 : : interp_(InterpMethod::CUBIC_SPLINE),
226 0 : nearestPointingTimeStamp_(0.0)
227 : {
228 0 : }
229 :
230 0 : void PI_Interp::setInterpMethod(InterpMethod m){
231 0 : interp_ = m;
232 0 : }
233 :
234 0 : PI_Interp::InterpMethod PI_Interp::getInterpMethod() const {
235 0 : return interp_ ;
236 : }
237 :
238 0 : void PI_Interp::setData(const Vector<PointingTimes> &antsTimes,
239 : const Vector<PointingDirs> &antsDirs,
240 : const Vector<bool> &antSelected){
241 :
242 : // TODO: check sizes of input vectors
243 :
244 : // Initialize
245 0 : clear();
246 0 : init(antSelected);
247 :
248 : // Count pointings per antenna
249 0 : auto nAnt = antsTimes.nelements();
250 0 : Vector<uInt> pointingsCount(nAnt,0);
251 0 : for (decltype(nAnt) antId = 0; antId < nAnt; antId++) {
252 0 : pointingsCount[antId] = antsTimes[antId].nelements();
253 : }
254 :
255 : // Compute spline coefficients for selected antennas,
256 : // where possible
257 0 : for (decltype(nAnt) antId = 0; antId < nAnt; ++antId) {
258 0 : if (not antSelected[antId]) continue;
259 0 : if (pointingsCount[antId] < 4) continue;
260 :
261 : // Resize private member: antsTimes_
262 : // auto antPointings = pointingsCount[antId];
263 : // antsTimes_[antId].resize(antPointings);
264 :
265 : // Store pointing times: required for interpolation
266 0 : antsTimes_[antId] = antsTimes[antId];
267 :
268 : // No need
269 : //dirPointing(i).resize(nPointingData(i));
270 :
271 : // Resize private member: antsSplinesCoeffs_
272 : //splineCoeff(i).resize(nPointingData(i) - 1);
273 0 : auto antPointings = pointingsCount[antId];
274 0 : auto antSplineSegments = antPointings - 1;
275 0 : antsSplinesCoeffs_[antId].resize(antSplineSegments);
276 :
277 : // No need
278 : //for (uInt j = 0; j < dirPointing(i).nelements(); ++j) {
279 : //dirPointing(i)(j).resize(2);
280 : //}
281 :
282 : // Resize private member: antsSplinesCoeffs_: continued
283 : ///for (uInt j = 0; j < splineCoeff(i).nelements(); ++j) {
284 0 : auto & antSplinesCoeffs = antsSplinesCoeffs_[antId];
285 0 : for (uInt s = 0; s < antSplineSegments; ++s){
286 : // splineCoeff(i)(j).resize(2);
287 0 : antSplinesCoeffs[s].resize(2);
288 : // splineCoeff(i)(j)(0).resize(4); // x
289 : // splineCoeff(i)(j)(1).resize(4); // y
290 0 : antSplinesCoeffs[s][0].resize(4);
291 0 : antSplinesCoeffs[s][1].resize(4);
292 : }
293 :
294 : // Copy input data into private members
295 : //Int npoi = nPointingData(i);
296 : //for (Int j = 0; j < npoi; ++j) {
297 :
298 : // Copy time into private member: done above
299 : // timePointing(i)(j) = time(i)(j);
300 :
301 : // Not needed
302 : //for (Int k = 0; k < 2; ++k) {
303 : // dirPointing(i)(j)(k) = dir(i)(j)(k);
304 : //}
305 : //}
306 :
307 :
308 : // calcSplineCoeff(timePointing(i), dirPointing(i), splineCoeff(i));
309 0 : const auto & antPointingTimes = antsTimes[antId];
310 0 : const auto & antPointingDirs = antsDirs[antId];
311 0 : computeSplineCoeffs(antPointingTimes,antPointingDirs,antSplinesCoeffs);
312 0 : isInterpolated_[antId] = true;
313 : }
314 :
315 0 : }
316 :
317 0 : void PI_Interp::init(const Vector<bool> &antSelected){
318 0 : isSelected_ = antSelected;
319 0 : auto nAnt = antSelected.nelements();
320 0 : antsTimes_.resize(nAnt);
321 0 : antsSplinesCoeffs_.resize(nAnt);
322 0 : isInterpolated_.resize(nAnt);
323 0 : isInterpolated_ = false;
324 0 : }
325 :
326 0 : void PI_Interp::clear(){
327 0 : init();
328 0 : }
329 :
330 0 : void PI_Interp::computeSplineCoeffs(const PointingTimes& time,
331 : const PointingDirs& dir,
332 : SplineCoeffs& coeff) {
333 0 : Vector<Double> h, vx, vy;
334 0 : Vector<Double> a;
335 0 : Vector<Double> c;
336 0 : Vector<Double> alpha, beta, gamma;
337 0 : Vector<Double> wx, wy;
338 0 : Vector<Double> ux, uy;
339 :
340 0 : const auto num_data = time.nelements();
341 0 : h.resize(num_data-1);
342 0 : vx.resize(num_data-1);
343 0 : vy.resize(num_data-1);
344 0 : a.resize(num_data-1);
345 0 : c.resize(num_data-1);
346 0 : alpha.resize(num_data-1);
347 0 : beta.resize(num_data-1);
348 0 : gamma.resize(num_data-1);
349 0 : wx.resize(num_data-1);
350 0 : wy.resize(num_data-1);
351 0 : ux.resize(num_data);
352 0 : uy.resize(num_data);
353 :
354 0 : h[0] = time[1] - time[0];
355 0 : for (uInt i = 1; i < num_data-1; ++i) {
356 0 : h[i] = time[i+1] - time[i];
357 0 : vx[i] = 6.0*((dir[i+1][0]-dir[i][0])/h[i] - (dir[i][0]-dir[i-1][0])/h[i-1]);
358 0 : vy[i] = 6.0*((dir[i+1][1]-dir[i][1])/h[i] - (dir[i][1]-dir[i-1][1])/h[i-1]);
359 0 : a[i] = 2.0*(time[i+1] - time[i-1]);
360 0 : c[i] = h[i];
361 0 : gamma[i] = c[i];
362 : }
363 0 : alpha[2] = c[1]/a[1];
364 0 : for (uInt i = 3; i < num_data-1; ++i) {
365 0 : alpha[i] = c[i-1]/(a[i-1] - alpha[i-1]*c[i-2]);
366 : }
367 0 : beta[1] = a[1];
368 0 : for (uInt i = 2; i < num_data-2; ++i) {
369 0 : beta[i] = c[i]/alpha[i+1];
370 : }
371 0 : beta[num_data-2] = a[num_data-2] - alpha[num_data-2] * c[num_data-3];
372 0 : wx[0] = 0.0;
373 0 : wx[1] = vx[1];
374 0 : wy[0] = 0.0;
375 0 : wy[1] = vy[1];
376 0 : for (uInt i = 2; i < num_data-1; ++i) {
377 0 : wx[i] = vx[i] - alpha[i]*wx[i-1];
378 0 : wy[i] = vy[i] - alpha[i]*wy[i-1];
379 : }
380 0 : ux[num_data-1] = 0.0;
381 0 : uy[num_data-1] = 0.0;
382 0 : for (uInt i = num_data-2; i >= 1; --i) {
383 0 : ux[i] = (wx[i] - gamma[i]*ux[i+1])/beta[i];
384 0 : uy[i] = (wy[i] - gamma[i]*uy[i+1])/beta[i];
385 : }
386 0 : ux[0] = 0.0;
387 0 : uy[0] = 0.0;
388 :
389 0 : for (uInt i = 0; i < num_data-1; ++i) {
390 0 : coeff[i][0][0] = dir[i][0];
391 0 : coeff[i][1][0] = dir[i][1];
392 0 : const auto dt = time[i+1]-time[i];
393 0 : coeff[i][0][1] = (dir[i+1][0]-dir[i][0])/dt - dt*(2.0*ux[i]+ux[i+1])/6.0;
394 0 : coeff[i][1][1] = (dir[i+1][1]-dir[i][1])/dt - dt*(2.0*uy[i]+uy[i+1])/6.0;
395 0 : coeff[i][0][2] = ux[i]/2.0;
396 0 : coeff[i][1][2] = uy[i]/2.0;
397 0 : coeff[i][0][3] = (ux[i+1]-ux[i])/dt/6.0;
398 0 : coeff[i][1][3] = (uy[i+1]-uy[i])/dt/6.0;
399 : }
400 0 : }
401 :
402 :
403 0 : Vector<Double> PI_Interp::pointingDir(int antId, double timeStamp) const {
404 : //const MSPointingColumns& mspc,
405 : //const Double& time,
406 : //const Int& index,
407 : //const Int& antid) {
408 :
409 :
410 0 : if ( not isInterpolated_[antId] ) {
411 0 : Vector<Double> pole(2,0.0);
412 0 : return pole;
413 : }
414 : // Search the segment to which the timestamp belongs
415 : // TODO: std::lowerbound
416 :
417 : //Int lastIndex = timePointing(antid).nelements() - 1;
418 0 : auto antPointingsCount = antsTimes_[antId].nelements();
419 0 : Int lastIndex = antPointingsCount - 1;
420 :
421 0 : Int aindex = lastIndex;
422 : //for (uInt i = 0; i < timePointing(antId).nelements(); ++i) {
423 0 : for (uInt i = 0; i < antPointingsCount; ++i) {
424 0 : if (timeStamp < antsTimes_[antId][i]) {
425 0 : aindex = i-1;
426 0 : break;
427 : }
428 : }
429 0 : if (aindex < 0) aindex = 0;
430 0 : if (lastIndex <= aindex) aindex = lastIndex - 1;
431 0 : nearestPointingTimeStamp_ = antsTimes_[antId][aindex];
432 :
433 : //auto const &coeff = splineCoeff(antid)(aindex);
434 0 : const auto & coeff = antsSplinesCoeffs_[antId][aindex];
435 : // Double dt = time - timePointing(antid)(aindex);
436 0 : Double dt = interp_ == InterpMethod::NEAREST ? 0.0 :
437 0 : timeStamp - antsTimes_[antId][aindex];
438 :
439 0 : Vector<Double> newdir(2);
440 : //newdir(0) = coeff(0)(0) + coeff(0)(1)*dt + coeff(0)(2)*dt*dt + coeff(0)(3)*dt*dt*dt;
441 : //newdir(1) = coeff(1)(0) + coeff(1)(1)*dt + coeff(1)(2)*dt*dt + coeff(1)(3)*dt*dt*dt;
442 0 : newdir(0) = coeff(0)(0) + dt*(coeff(0)(1) + dt*(coeff(0)(2) + dt*coeff(0)(3)));
443 0 : newdir(1) = coeff(1)(0) + dt*(coeff(1)(1) + dt*(coeff(1)(2) + dt*coeff(1)(3)));
444 :
445 0 : return newdir;
446 :
447 : // Quantity rDirLon(newdir(0), "rad");
448 : // Quantity rDirLat(newdir(1), "rad");
449 : // auto const &directionMeasColumn = mspc.directionMeasCol();
450 : // MDirection::Ref rf(directionMeasColumn.measDesc().getRefCode());
451 : // if (directionMeasColumn.measDesc().isRefCodeVariable()) {
452 : // rf = mspc.directionMeas(index).getRef();
453 : // }
454 : //return MDirection(rDirLon, rDirLat, rf);
455 : }
456 :
457 0 : double PI_Interp::nearestPointingTimeStamp() const {
458 0 : return nearestPointingTimeStamp_;
459 : }
460 : // ------------------------- Factories ---------------------------------
461 :
462 : // ---- Vi2Factory
463 0 : PI_Vi2Factory::PointingInterpolationVi2Factory(Record const &configuration,
464 0 : ViImplementation2 *inputVII) :
465 : inputVII_p(inputVII),
466 0 : configuration_(configuration)
467 : {
468 :
469 0 : }
470 :
471 0 : PI_Vi2Factory::PointingInterpolationVi2Factory(const Record & /*configuration*/,
472 : const MeasurementSet *ms,
473 : const SortColumns &sortColumns,
474 : Double timeInterval,
475 0 : Bool isWritable) :
476 0 : inputVII_p(nullptr) {
477 :
478 0 : inputVII_p = new VisibilityIteratorImpl2(
479 0 : Block<const MeasurementSet *>(1, ms),
480 0 : sortColumns, timeInterval, isWritable);
481 0 : }
482 :
483 0 : PI_Vi2Factory::~PointingInterpolationVi2Factory() {
484 0 : }
485 :
486 0 : ViImplementation2 * PI_Vi2Factory::createVi() const {
487 0 : return new PI_TVI(inputVII_p);
488 : }
489 :
490 : // ---- TVILayerFactory
491 0 : PI_TVILayerFactory::PointingInterpolationTVILayerFactory(
492 0 : Record const &configuration) :
493 0 : ViiLayerFactory()
494 : {
495 0 : configuration_p = configuration;
496 0 : }
497 :
498 0 : PI_TVILayerFactory::~PointingInterpolationTVILayerFactory()
499 : {
500 :
501 0 : }
502 :
503 : ViImplementation2*
504 0 : PI_TVILayerFactory::createInstance(ViImplementation2* vii0) const {
505 : // Make the PI_TVI, using supplied ViImplementation2, and return it
506 0 : PointingInterpolationVi2Factory factory(configuration_p, vii0);
507 0 : ViImplementation2 *vii = nullptr;
508 : try {
509 0 : vii = factory.createVi();
510 0 : } catch (...) {
511 0 : if (vii0) {
512 0 : delete vii0;
513 : }
514 0 : throw;
515 : }
516 0 : return vii;
517 : }
518 :
519 :
520 :
521 :
|
