Line data Source code
1 : //# CTPatchedInterp.cc: Implementation of CTPatchedInterp.h
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,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 adressed 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 : //#
27 :
28 : #include <synthesis/CalTables/CTPatchedInterp.h>
29 : #include <synthesis/CalTables/CTIter.h>
30 :
31 : #include <casacore/scimath/Mathematics/InterpolateArray1D.h>
32 : #include <casacore/casa/OS/Path.h>
33 : #include <casacore/casa/Utilities/GenSort.h>
34 : #include <casacore/casa/aips.h>
35 :
36 : #define CTPATCHEDINTERPVERB false
37 :
38 : #define NEAREST InterpolateArray1D<Double,Float>::nearestNeighbour
39 : #define LINEAR InterpolateArray1D<Double,Float>::linear
40 : #define CUBIC InterpolateArray1D<Double,Float>::cubic
41 : #define SPLINE InterpolateArray1D<Double,Float>::spline
42 :
43 : #include <casacore/ms/MSOper/MSMetaData.h>
44 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
45 : #include <casacore/casa/Logging/LogMessage.h>
46 : #include <casacore/casa/Logging/LogIO.h>
47 :
48 : using namespace casacore;
49 : namespace casa { //# NAMESPACE CASA - BEGIN
50 :
51 : // Ctor
52 0 : CTPatchedInterp::CTPatchedInterp(NewCalTable& ct,
53 : VisCalEnum::MatrixType mtype,
54 : Int nPar,
55 : const String& timetype,
56 : const String& freqtype,
57 : const String& fieldtype,
58 : Vector<Int> spwmap,
59 : Vector<Int> fldmap,
60 0 : const CTTIFactoryPtr cttifactoryptr) :
61 0 : ct_(ct),
62 0 : ctname_(),
63 0 : msmc_(NULL),
64 0 : mtype_(mtype),
65 0 : isCmplx_(false),
66 0 : nPar_(nPar),
67 0 : nFPar_(nPar),
68 0 : timeType_(timetype),
69 0 : freqTypeStr_(freqtype),
70 0 : relativeFreq_(freqtype.contains("rel")),
71 0 : freqInterpMethod0_(ftype(freqTypeStr_)),
72 0 : freqInterpMethod_(freqInterpMethod0_),
73 0 : freqInterpMethodVec_(),
74 0 : byObs_(timetype.contains("perobs")), // detect slicing by obs
75 0 : byScan_(timetype.contains("perscan")), // detect slicing by scan
76 0 : byField_(fieldtype=="nearest" || fieldtype=="map"),
77 0 : nChanIn_(),
78 0 : freqIn_(),
79 0 : nMSTimeSeg_(1), // byObs_?ct.observation().nrow():1), // assume CT shapes for MS shapes
80 0 : nMSFld_(ct.field().nrow()),
81 0 : nMSSpw_(ct.spectralWindow().nrow()),
82 0 : nMSAnt_(ct.antenna().nrow()),
83 0 : altFld_(),
84 0 : nCTTimeSeg_(1), // byObs_?ct.observation().nrow():1),
85 0 : nCTFld_(byField_?ct.field().nrow():1),
86 0 : nCTSpw_(ct.spectralWindow().nrow()),
87 0 : nCTAnt_(ct.antenna().nrow()),
88 0 : nCTElem_(0),
89 0 : spwInOK_(),
90 0 : fldMap_(),
91 0 : spwMap_(),
92 0 : antMap_(),
93 0 : elemMap_(),
94 0 : conjTab_(),
95 0 : result_(),
96 0 : resFlag_(),
97 0 : tI_(),
98 0 : tIdel_(),
99 0 : tIMissingLogged_(),
100 0 : lastFld_(ct.spectralWindow().nrow(),-1),
101 0 : lastObs_(ct.spectralWindow().nrow(),-1),
102 0 : lastScan_(ct.spectralWindow().nrow(),-1),
103 0 : cttifactoryptr_(cttifactoryptr)
104 : {
105 : if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::CTPatchedInterp(<no MS>)" << endl;
106 :
107 : // cout << "freqInterpMethod_ = " << freqInterpMethod_ << endl;
108 :
109 0 : freqInterpMethodVec_.resize(nMSSpw_);
110 0 : freqInterpMethodVec_.set(freqInterpMethod0_);
111 :
112 0 : switch(mtype_) {
113 0 : case VisCalEnum::GLOBAL: {
114 :
115 0 : throw(AipsError("CTPatchedInterp::ctor: No non-Mueller/Jones yet."));
116 :
117 : nCTElem_=1;
118 : nMSElem_=1;
119 : break;
120 : }
121 0 : case VisCalEnum::MUELLER: {
122 0 : nCTElem_=nCTAnt_*(nCTAnt_+1)/2;
123 0 : nMSElem_=nMSAnt_*(nMSAnt_+1)/2;
124 0 : break;
125 : }
126 0 : case VisCalEnum::JONES: {
127 0 : nCTElem_=nCTAnt_;
128 0 : nMSElem_=nMSAnt_;
129 0 : break;
130 : }
131 : }
132 :
133 : // How many _Float_ parameters?
134 0 : isCmplx_=(ct_.keywordSet().asString("ParType")=="Complex");
135 0 : if (isCmplx_) // Complex input
136 0 : nFPar_*=2; // interpolating 2X as many Float values
137 :
138 :
139 : // Set channel/freq info
140 0 : CTSpWindowColumns ctspw(ct_.spectralWindow());
141 0 : ctspw.numChan().getColumn(nChanIn_);
142 0 : freqIn_.resize(nCTSpw_);
143 0 : for (uInt iCTspw=0;iCTspw<ctspw.nrow();++iCTspw) {
144 0 : ctspw.chanFreq().get(iCTspw,freqIn_(iCTspw),true);
145 0 : if (relativeFreq_) {
146 0 : Vector<Double>& fIn(freqIn_(iCTspw));
147 0 : fIn-=mean(fIn); // assume mean freq is center, and apply offset
148 : // Flip LSB
149 0 : if (fIn.nelements()>1 && fIn(0)>fIn(1)) {
150 : //cout << " spw=" << iCTspw << " LSB!" << endl;
151 0 : fIn*=Double(-1.0);
152 : }
153 : }
154 : //cout << " freqIn_(" << iCTspw << ")=" << freqIn_(iCTspw) << endl;
155 : }
156 :
157 : // byScan_ supercedes byObs_, so turn OFF byObs_, if necessary
158 0 : if (byObs_ && byScan_) {
159 0 : byObs_=false;
160 0 : LogIO log;
161 0 : ostringstream msg;
162 0 : msg << "For " << Path(ct_.tableName()).baseName().before(".tempMemCal")
163 0 : << " 'perscan' interpolation supercedes 'perobs' interpolation. ";
164 0 : log << msg.str() << LogIO::WARN;
165 0 : }
166 :
167 : // Manage 'byObs_' carefully
168 0 : if (byObs_) {
169 :
170 0 : Int nMSObs=ct_.observation().nrow(); // assume CT shapes for MS shapes
171 0 : Int nCTObs=ct_.observation().nrow();
172 :
173 : // Count _available_ obsids in caltable
174 0 : ROCTMainColumns ctmc(ct_);
175 0 : Vector<Int> obsid;
176 0 : ctmc.obsId().getColumn(obsid);
177 0 : Int nctobsavail=genSort(obsid,Sort::Ascending,
178 0 : (Sort::QuickSort | Sort::NoDuplicates));
179 :
180 :
181 0 : LogIO log;
182 0 : ostringstream msg;
183 :
184 0 : if (nctobsavail==1) {
185 0 : byObs_=false;
186 0 : msg << "Only one ObsId found in "
187 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
188 0 : << "; ignoring 'perobs' interpolation.";
189 0 : log << msg.str() << LogIO::WARN;
190 : }
191 : else {
192 :
193 : // Verify consistency between CT and MS
194 0 : if (nctobsavail==nCTObs &&
195 : nctobsavail==nMSObs) {
196 : // Everything ok
197 0 : nCTTimeSeg_=nCTObs;
198 0 : nMSTimeSeg_=nMSObs;
199 : }
200 : else {
201 : // only 1 obs, or available nobs doesn't match MS
202 0 : byObs_=false;
203 0 : msg << "Multiple ObsIds found in "
204 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
205 : << ", but they do not match the MS ObsIds;"
206 0 : << " turning off 'perobs'.";
207 0 : log << msg.str() << LogIO::WARN;
208 : }
209 : }
210 :
211 0 : }
212 :
213 : // Manage 'byScan_' even more carefully
214 0 : if (byScan_) {
215 :
216 0 : LogIO log;
217 0 : ostringstream msg;
218 :
219 : // Count _availale_ scan numbers in caltable
220 0 : ROCTMainColumns ctmc(ct_);
221 0 : Vector<Int> scan;
222 : Int minScan, maxScan;
223 0 : ctmc.scanNo().getColumn(scan);
224 0 : minMax(minScan, maxScan, scan);
225 0 : if (minScan == -1) {
226 0 : byScan_=false;
227 0 : msg << "No scan numbers found in "
228 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
229 0 : << "; ignoring 'perscan' interpolation.";
230 0 : log << msg.str() << LogIO::WARN;
231 : } else {
232 0 : nCTTimeSeg_=maxScan+1;
233 0 : nMSTimeSeg_=maxScan+1; // assume CT shapes for MS shapes
234 : }
235 0 : }
236 :
237 : // Initialize caltable slices
238 0 : sliceTable();
239 :
240 : // Set spwmap
241 0 : setSpwMap(spwmap);
242 :
243 : // Set fldmap
244 0 : if (byField_) {
245 0 : if (fieldtype=="map")
246 0 : setFldMap(fldmap); // Use specified map
247 : else
248 0 : setFldMap(ct.field()); // Use CalTable's fields ('nearest')
249 : }
250 : else
251 0 : setDefFldMap();
252 :
253 : // Set defaultmaps
254 0 : setDefAntMap();
255 0 : setElemMap();
256 :
257 : // Resize working arrays
258 0 : result_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
259 0 : resFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
260 0 : timeResult_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
261 0 : timeResFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
262 0 : freqResult_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
263 0 : freqResFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
264 :
265 : // Figure out where we can duplicate field interpolators
266 0 : calcAltFld();
267 :
268 : // Setup mapped interpolators
269 : // TBD: defer this to later, so that spwmap, etc. can be revised
270 : // before committing to the interpolation engines
271 0 : makeInterpolators();
272 :
273 : // state();
274 0 : }
275 :
276 0 : CTPatchedInterp::CTPatchedInterp(NewCalTable& ct,
277 : VisCalEnum::MatrixType mtype,
278 : Int nPar,
279 : const String& timetype,
280 : const String& freqtype,
281 : const String& fieldtype,
282 : const MeasurementSet& ms,
283 : const MSMetaInfoForCal& msmc,
284 : Vector<Int> spwmap,
285 0 : const CTTIFactoryPtr cttifactoryptr) :
286 0 : ct_(ct),
287 0 : ctname_(),
288 0 : msmc_(&msmc),
289 0 : mtype_(mtype),
290 0 : isCmplx_(false),
291 0 : nPar_(nPar),
292 0 : nFPar_(nPar),
293 0 : timeType_(timetype),
294 0 : freqTypeStr_(freqtype),
295 0 : relativeFreq_(freqtype.contains("rel")),
296 0 : freqInterpMethod0_(ftype(freqTypeStr_)),
297 0 : freqInterpMethod_(freqInterpMethod0_),
298 0 : freqInterpMethodVec_(),
299 0 : byObs_(timetype.contains("perobs")), // detect slicing by obs
300 0 : byScan_(timetype.contains("perscan")), // detect slicing by scan
301 0 : byField_(fieldtype=="nearest"), // for now we are NOT slicing by field
302 0 : nChanIn_(),
303 0 : freqIn_(),
304 0 : nMSTimeSeg_(1), // byObs_?ms.observation().nrow():1),
305 0 : nMSFld_(ms.field().nrow()),
306 0 : nMSSpw_(ms.spectralWindow().nrow()),
307 0 : nMSAnt_(ms.antenna().nrow()),
308 0 : altFld_(),
309 0 : nCTTimeSeg_(1), // byObs_?ct.observation().nrow():1),
310 0 : nCTFld_(byField_?ct.field().nrow():1),
311 0 : nCTSpw_(ct.spectralWindow().nrow()),
312 0 : nCTAnt_(ct.antenna().nrow()),
313 0 : nCTElem_(0),
314 0 : spwInOK_(),
315 0 : fldMap_(),
316 0 : spwMap_(),
317 0 : antMap_(),
318 0 : elemMap_(),
319 0 : conjTab_(),
320 0 : result_(),
321 0 : resFlag_(),
322 0 : tI_(),
323 0 : tIdel_(),
324 0 : tIMissingLogged_(),
325 0 : lastFld_(ms.spectralWindow().nrow(),-1),
326 0 : lastObs_(ms.spectralWindow().nrow(),-1),
327 0 : lastScan_(ms.spectralWindow().nrow(),-1),
328 0 : cttifactoryptr_(cttifactoryptr)
329 : {
330 :
331 : if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::CTPatchedInterp(CT,MS)" << endl;
332 :
333 : //freqInterpMethod_=ftype(freqTypeStr_);
334 :
335 0 : freqInterpMethodVec_.resize(nMSSpw_);
336 0 : freqInterpMethodVec_.set(freqInterpMethod0_);
337 :
338 : // cout << "freqInterpMethod_ = " << freqInterpMethod_ << endl;
339 :
340 0 : switch(mtype_) {
341 0 : case VisCalEnum::GLOBAL: {
342 :
343 0 : throw(AipsError("CTPatchedInterp::ctor: No non-Mueller/Jones yet."));
344 :
345 : nCTElem_=1;
346 : nMSElem_=1;
347 : break;
348 : }
349 0 : case VisCalEnum::MUELLER: {
350 0 : nCTElem_=nCTAnt_*(nCTAnt_+1)/2;
351 0 : nMSElem_=nMSAnt_*(nMSAnt_+1)/2;
352 0 : break;
353 : }
354 0 : case VisCalEnum::JONES: {
355 0 : nCTElem_=nCTAnt_;
356 0 : nMSElem_=nMSAnt_;
357 0 : break;
358 : }
359 : }
360 :
361 : // How many _Float_ parameters?
362 0 : isCmplx_=(ct_.keywordSet().asString("ParType")=="Complex"); // Complex input
363 0 : if (isCmplx_)
364 0 : nFPar_*=2; // interpolating 2X as many Float values
365 :
366 : // Set channel/freq info
367 0 : CTSpWindowColumns ctspw(ct_.spectralWindow());
368 0 : ctspw.numChan().getColumn(nChanIn_);
369 0 : freqIn_.resize(nCTSpw_);
370 0 : for (uInt iCTspw=0;iCTspw<ctspw.nrow();++iCTspw) {
371 0 : ctspw.chanFreq().get(iCTspw,freqIn_(iCTspw),true);
372 0 : if (relativeFreq_) {
373 0 : Vector<Double>& fIn(freqIn_(iCTspw));
374 0 : fIn-=mean(fIn); // assume mean freq is center, and apply offset
375 : // Flip LSB
376 0 : if (fIn.nelements()>1 && fIn(0)>fIn(1)) {
377 : //cout << " spw=" << iCTspw << " LSB!" << endl;
378 0 : fIn*=Double(-1.0);
379 : }
380 : }
381 : //cout << " freqIn_(" << iCTspw << ")=" << freqIn_(iCTspw) << endl;
382 : }
383 :
384 : // byScan_ supercedes byObs_, so turn OFF byObs_, if necessary
385 0 : if (byObs_ && byScan_) {
386 0 : byObs_=false;
387 0 : LogIO log;
388 0 : ostringstream msg;
389 0 : msg << "For " << Path(ct_.tableName()).baseName().before(".tempMemCal")
390 0 : << " 'perscan' interpolation supercedes 'perobs' interpolation. ";
391 0 : log << msg.str() << LogIO::WARN;
392 0 : }
393 :
394 : // Manage 'byObs_' carefully
395 0 : if (byObs_) {
396 0 : Int nMSObs=ms.observation().nrow();
397 0 : Int nCTObs=ct_.observation().nrow();
398 :
399 : // Count _available_ obsids in caltable
400 0 : ROCTMainColumns ctmc(ct_);
401 0 : Vector<Int> obsid;
402 0 : ctmc.obsId().getColumn(obsid);
403 0 : Int nctobsavail=genSort(obsid,Sort::Ascending,
404 0 : (Sort::QuickSort | Sort::NoDuplicates));
405 :
406 :
407 0 : LogIO log;
408 0 : ostringstream msg;
409 :
410 0 : if (nctobsavail==1) {
411 0 : byObs_=false;
412 0 : msg << "Only one ObsId found in "
413 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
414 0 : << "; ignoring 'perobs' interpolation.";
415 0 : log << msg.str() << LogIO::WARN;
416 : }
417 : else {
418 :
419 : // Verify consistency between CT and MS
420 0 : if (nctobsavail==nCTObs &&
421 : nctobsavail==nMSObs) {
422 : // Everything ok
423 0 : nCTTimeSeg_=nCTObs;
424 0 : nMSTimeSeg_=nMSObs;
425 : }
426 : else {
427 : // only 1 obs, or available nobs doesn't match MS
428 0 : byObs_=false;
429 0 : msg << "Multiple ObsIds found in "
430 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
431 : << ", but they do not match the MS ObsIds;"
432 0 : << " turning off 'perobs'.";
433 0 : log << msg.str() << LogIO::WARN;
434 : }
435 : }
436 0 : }
437 :
438 : // Manage 'byScan_' even more carefully
439 0 : if (byScan_) {
440 :
441 0 : LogIO log;
442 0 : ostringstream msg;
443 :
444 : // Count _availale_ scan numbers in caltable
445 0 : ROCTMainColumns ctmc(ct_);
446 0 : Vector<Int> scan;
447 : Int minScan, maxScan;
448 0 : ctmc.scanNo().getColumn(scan);
449 0 : minMax(minScan, maxScan, scan);
450 0 : if (minScan == -1) {
451 0 : byScan_=false;
452 0 : msg << "No scan numbers found in "
453 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
454 0 : << "; ignoring 'perscan' interpolation.";
455 0 : log << msg.str() << LogIO::WARN;
456 : } else {
457 0 : std::set<Int> scanNumbers = msmc.msmd().getScanNumbers(0,0);
458 0 : nCTTimeSeg_=maxScan+1;
459 0 : nMSTimeSeg_=*scanNumbers.rbegin()+1;
460 :
461 0 : }
462 0 : }
463 :
464 : // Initialize caltable slices
465 0 : sliceTable();
466 :
467 : // Set spwmap
468 0 : setSpwMap(spwmap);
469 :
470 : // Set fldmap
471 0 : if (byField_)
472 0 : setFldMap(ms.field()); // on a trial basis
473 : else
474 0 : setDefFldMap();
475 :
476 : // Set defaultmaps
477 0 : setDefAntMap();
478 0 : setElemMap();
479 :
480 : // Resize working arrays
481 0 : result_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
482 0 : resFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
483 0 : timeResult_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
484 0 : timeResFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
485 0 : freqResult_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
486 0 : freqResFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
487 :
488 : // Figure out where we can duplicate field interpolators
489 0 : calcAltFld();
490 :
491 : // Setup mapped interpolators
492 : // TBD: defer this to later, so that spwmap, etc. can be revised
493 : // before committing to the interpolation engines
494 0 : makeInterpolators();
495 :
496 : // state();
497 :
498 0 : }
499 :
500 0 : CTPatchedInterp::CTPatchedInterp(NewCalTable& ct,
501 : VisCalEnum::MatrixType mtype,
502 : Int nPar,
503 : const String& timetype,
504 : const String& freqtype,
505 : const String& fieldtype,
506 : const MSColumns& mscol,
507 : Vector<Int> spwmap,
508 0 : const CTTIFactoryPtr cttifactoryptr) :
509 0 : ct_(ct),
510 0 : ctname_(),
511 0 : msmc_(NULL),
512 0 : mtype_(mtype),
513 0 : isCmplx_(false),
514 0 : nPar_(nPar),
515 0 : nFPar_(nPar),
516 0 : timeType_(timetype),
517 0 : freqTypeStr_(freqtype),
518 0 : relativeFreq_(freqtype.contains("rel")),
519 0 : freqInterpMethod0_(ftype(freqTypeStr_)),
520 0 : freqInterpMethod_(freqInterpMethod0_),
521 0 : freqInterpMethodVec_(),
522 0 : byObs_(false), // turn off for old-fashioned
523 0 : byScan_(false), // turn off for old-fashioned
524 0 : byField_(fieldtype=="nearest"), // for now we are NOT slicing by field
525 0 : nChanIn_(),
526 0 : freqIn_(),
527 0 : nMSFld_(mscol.field().nrow()),
528 0 : nMSSpw_(mscol.spectralWindow().nrow()),
529 0 : nMSAnt_(mscol.antenna().nrow()),
530 0 : altFld_(),
531 0 : nCTFld_(byField_?ct.field().nrow():1),
532 0 : nCTSpw_(ct.spectralWindow().nrow()),
533 0 : nCTAnt_(ct.antenna().nrow()),
534 0 : nCTElem_(0),
535 0 : spwInOK_(),
536 0 : fldMap_(),
537 0 : spwMap_(),
538 0 : antMap_(),
539 0 : elemMap_(),
540 0 : conjTab_(),
541 0 : result_(),
542 0 : resFlag_(),
543 0 : tI_(),
544 0 : tIdel_(),
545 0 : tIMissingLogged_(),
546 0 : lastFld_(mscol.spectralWindow().nrow(),-1),
547 0 : lastObs_(mscol.spectralWindow().nrow(),-1),
548 0 : lastScan_(mscol.spectralWindow().nrow(),-1),
549 0 : cttifactoryptr_(cttifactoryptr)
550 : {
551 : if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::CTPatchedInterp(mscol)" << endl;
552 :
553 : // cout << "freqInterpMethod_ = " << freqInterpMethod_ << endl;
554 0 : freqInterpMethodVec_.resize(nMSSpw_);
555 0 : freqInterpMethodVec_.set(freqInterpMethod0_);
556 :
557 0 : switch(mtype_) {
558 0 : case VisCalEnum::GLOBAL: {
559 :
560 0 : throw(AipsError("CTPatchedInterp::ctor: No non-Mueller/Jones yet."));
561 :
562 : nCTElem_=1;
563 : nMSElem_=1;
564 : break;
565 : }
566 0 : case VisCalEnum::MUELLER: {
567 0 : nCTElem_=nCTAnt_*(nCTAnt_+1)/2;
568 0 : nMSElem_=nMSAnt_*(nMSAnt_+1)/2;
569 0 : break;
570 : }
571 0 : case VisCalEnum::JONES: {
572 0 : nCTElem_=nCTAnt_;
573 0 : nMSElem_=nMSAnt_;
574 0 : break;
575 : }
576 : }
577 :
578 : // How many _Float_ parameters?
579 0 : isCmplx_=(ct_.keywordSet().asString("ParType")=="Complex");
580 0 : if (isCmplx_) // Complex input
581 0 : nFPar_*=2; // interpolating 2X as many Float values
582 :
583 : // Set channel/freq info
584 0 : CTSpWindowColumns ctspw(ct_.spectralWindow());
585 0 : ctspw.numChan().getColumn(nChanIn_);
586 0 : freqIn_.resize(nCTSpw_);
587 0 : for (uInt iCTspw=0;iCTspw<ctspw.nrow();++iCTspw) {
588 0 : ctspw.chanFreq().get(iCTspw,freqIn_(iCTspw),true);
589 0 : if (relativeFreq_) {
590 0 : Vector<Double>& fIn(freqIn_(iCTspw));
591 0 : fIn-=mean(fIn); // assume mean freq is center, and apply offset
592 : // Flip LSB
593 0 : if (fIn.nelements()>1 && fIn(0)>fIn(1)) {
594 : //cout << " spw=" << iCTspw << " LSB!" << endl;
595 0 : fIn*=Double(-1.0);
596 : }
597 : }
598 : // cout << " freqIn_(" << iCTspw << ")=" << freqIn_(iCTspw) << endl;
599 : }
600 :
601 : // Initialize caltable slices
602 0 : sliceTable();
603 :
604 : // Set spwmap
605 0 : setSpwMap(spwmap);
606 :
607 : // Set fldmap
608 0 : if (byField_)
609 0 : setFldMap(mscol.field()); // on a trial basis
610 : else
611 0 : setDefFldMap();
612 :
613 :
614 : // Set defaultmaps
615 0 : setDefAntMap();
616 0 : setElemMap();
617 :
618 : // Resize working arrays
619 0 : result_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
620 0 : resFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
621 0 : timeResult_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
622 0 : timeResFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
623 0 : freqResult_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
624 0 : freqResFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
625 :
626 : // Figure out where we can duplicate field interpolators
627 0 : calcAltFld();
628 :
629 : // Setup mapped interpolators
630 : // TBD: defer this to later, so that spwmap, etc. can be revised
631 : // before committing to the interpolation engines
632 0 : makeInterpolators();
633 :
634 : // state();
635 :
636 0 : }
637 :
638 0 : CTPatchedInterp::CTPatchedInterp(NewCalTable& ct,
639 : VisCalEnum::MatrixType mtype,
640 : Int nPar,
641 : const String& timetype,
642 : const String& freqtype,
643 : const String& fieldtype,
644 : const MeasurementSet& ms,
645 : Vector<Int> spwmap,
646 0 : const CTTIFactoryPtr cttifactoryptr) :
647 0 : ct_(ct),
648 0 : ctname_(),
649 0 : msmc_(NULL),
650 0 : mtype_(mtype),
651 0 : isCmplx_(false),
652 0 : nPar_(nPar),
653 0 : nFPar_(nPar),
654 0 : timeType_(timetype),
655 0 : freqTypeStr_(freqtype),
656 0 : relativeFreq_(freqtype.contains("rel")),
657 0 : freqInterpMethod0_(ftype(freqTypeStr_)),
658 0 : freqInterpMethod_(freqInterpMethod0_),
659 0 : freqInterpMethodVec_(),
660 0 : byObs_(timetype.contains("perobs")), // detect slicing by obs
661 0 : byScan_(timetype.contains("perscan")), // detect slicing by scan
662 0 : byField_(fieldtype=="nearest"), // for now we are NOT slicing by field
663 0 : nChanIn_(),
664 0 : freqIn_(),
665 0 : nMSTimeSeg_(1), // byObs_?ms.observation().nrow():1),
666 0 : nMSFld_(ms.field().nrow()),
667 0 : nMSSpw_(ms.spectralWindow().nrow()),
668 0 : nMSAnt_(ms.antenna().nrow()),
669 0 : altFld_(),
670 0 : nCTTimeSeg_(1), // byObs_?ct.observation().nrow():1),
671 0 : nCTFld_(byField_?ct.field().nrow():1),
672 0 : nCTSpw_(ct.spectralWindow().nrow()),
673 0 : nCTAnt_(ct.antenna().nrow()),
674 0 : nCTElem_(0),
675 0 : spwInOK_(),
676 0 : fldMap_(),
677 0 : spwMap_(),
678 0 : antMap_(),
679 0 : elemMap_(),
680 0 : conjTab_(),
681 0 : result_(),
682 0 : resFlag_(),
683 0 : tI_(),
684 0 : tIdel_(),
685 0 : tIMissingLogged_(),
686 0 : lastFld_(ms.spectralWindow().nrow(),-1),
687 0 : lastObs_(ms.spectralWindow().nrow(),-1),
688 0 : lastScan_(ms.spectralWindow().nrow(),-1),
689 0 : cttifactoryptr_(cttifactoryptr)
690 : {
691 :
692 : if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::CTPatchedInterp(CT,MS)" << endl;
693 :
694 : //freqInterpMethod_=ftype(freqTypeStr_);
695 :
696 0 : freqInterpMethodVec_.resize(nMSSpw_);
697 0 : freqInterpMethodVec_.set(freqInterpMethod0_);
698 :
699 : // cout << "freqInterpMethod_ = " << freqInterpMethod_ << endl;
700 :
701 0 : switch(mtype_) {
702 0 : case VisCalEnum::GLOBAL: {
703 :
704 0 : throw(AipsError("CTPatchedInterp::ctor: No non-Mueller/Jones yet."));
705 :
706 : nCTElem_=1;
707 : nMSElem_=1;
708 : break;
709 : }
710 0 : case VisCalEnum::MUELLER: {
711 0 : nCTElem_=nCTAnt_*(nCTAnt_+1)/2;
712 0 : nMSElem_=nMSAnt_*(nMSAnt_+1)/2;
713 0 : break;
714 : }
715 0 : case VisCalEnum::JONES: {
716 0 : nCTElem_=nCTAnt_;
717 0 : nMSElem_=nMSAnt_;
718 0 : break;
719 : }
720 : }
721 :
722 : // How many _Float_ parameters?
723 0 : isCmplx_=(ct_.keywordSet().asString("ParType")=="Complex"); // Complex input
724 0 : if (isCmplx_)
725 0 : nFPar_*=2; // interpolating 2X as many Float values
726 :
727 : // Set channel/freq info
728 0 : CTSpWindowColumns ctspw(ct_.spectralWindow());
729 0 : ctspw.numChan().getColumn(nChanIn_);
730 0 : freqIn_.resize(nCTSpw_);
731 0 : for (uInt iCTspw=0;iCTspw<ctspw.nrow();++iCTspw) {
732 0 : ctspw.chanFreq().get(iCTspw,freqIn_(iCTspw),true);
733 0 : if (relativeFreq_) {
734 0 : Vector<Double>& fIn(freqIn_(iCTspw));
735 0 : fIn-=mean(fIn); // assume mean freq is center, and apply offset
736 : // Flip LSB
737 0 : if (fIn.nelements()>1 && fIn(0)>fIn(1)) {
738 : //cout << " spw=" << iCTspw << " LSB!" << endl;
739 0 : fIn*=Double(-1.0);
740 : }
741 : }
742 : //cout << " freqIn_(" << iCTspw << ")=" << freqIn_(iCTspw) << endl;
743 : }
744 :
745 : // Manage 'byObs_' carefully
746 0 : if (byObs_) {
747 0 : Int nMSObs=ms.observation().nrow();
748 0 : Int nCTObs=ct_.observation().nrow();
749 :
750 : // Count _available_ obsids in caltable
751 0 : ROCTMainColumns ctmc(ct_);
752 0 : Vector<Int> obsid;
753 0 : ctmc.obsId().getColumn(obsid);
754 0 : Int nctobsavail=genSort(obsid,Sort::Ascending,
755 0 : (Sort::QuickSort | Sort::NoDuplicates));
756 :
757 :
758 0 : LogIO log;
759 0 : ostringstream msg;
760 :
761 0 : if (nctobsavail==1) {
762 0 : byObs_=false;
763 0 : msg << "Only one ObsId found in "
764 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
765 0 : << "; ignoring 'perobs' interpolation.";
766 0 : log << msg.str() << LogIO::WARN;
767 : }
768 : else {
769 :
770 : // Verify consistency between CT and MS
771 0 : if (nctobsavail==nCTObs &&
772 : nctobsavail==nMSObs) {
773 : // Everything ok
774 0 : nCTTimeSeg_=nCTObs;
775 0 : nMSTimeSeg_=nMSObs;
776 : }
777 : else {
778 : // only 1 obs, or available nobs doesn't match MS
779 0 : byObs_=false;
780 0 : msg << "Multiple ObsIds found in "
781 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
782 : << ", but they do not match the MS ObsIds;"
783 0 : << " turning off 'perobs'.";
784 0 : log << msg.str() << LogIO::WARN;
785 : }
786 : }
787 0 : }
788 :
789 : // Initialize caltable slices
790 0 : sliceTable();
791 :
792 : // Set spwmap
793 0 : setSpwMap(spwmap);
794 :
795 : // Set fldmap
796 0 : if (byField_)
797 0 : setFldMap(ms.field()); // on a trial basis
798 : else
799 0 : setDefFldMap();
800 :
801 : // Set defaultmaps
802 0 : setDefAntMap();
803 0 : setElemMap();
804 :
805 : // Resize working arrays
806 0 : result_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
807 0 : resFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
808 0 : timeResult_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
809 0 : timeResFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
810 0 : freqResult_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
811 0 : freqResFlag_.resize(nMSSpw_,nMSFld_,nMSTimeSeg_);
812 :
813 : // Figure out where we can duplicate field interpolators
814 0 : calcAltFld();
815 :
816 : // Setup mapped interpolators
817 : // TBD: defer this to later, so that spwmap, etc. can be revised
818 : // before committing to the interpolation engines
819 0 : makeInterpolators();
820 :
821 : //state();
822 :
823 0 : }
824 :
825 :
826 : // Destructor
827 0 : CTPatchedInterp::~CTPatchedInterp() {
828 :
829 : if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::~CTPatchedInterp()" << endl;
830 :
831 : {
832 0 : IPosition sh(tI_.shape());
833 0 : for (Int l=0;l<sh(3);++l)
834 0 : for (Int k=0;k<sh(2);++k)
835 0 : for (Int j=0;j<sh(1);++j)
836 0 : for (Int i=0;i<sh(0);++i) {
837 0 : IPosition ip(4,i,j,k,l);
838 0 : if (tIdel_(ip))
839 0 : delete tI_(ip);
840 0 : }
841 0 : tI_.resize();
842 0 : }
843 : {
844 0 : IPosition sh(ctSlices_.shape());
845 0 : for (Int l=0;l<sh(3);++l)
846 0 : for (Int k=0;k<sh(2);++k)
847 0 : for (Int j=0;j<sh(1);++j)
848 0 : for (Int i=0;i<sh(0);++i) {
849 0 : IPosition ip(4,i,j,k,l);
850 0 : if (ctSlices_(ip))
851 0 : delete ctSlices_(ip);
852 0 : }
853 0 : ctSlices_.resize();
854 0 : }
855 0 : }
856 :
857 0 : Bool CTPatchedInterp::interpolate(Int msobs, Int msscan, Int msfld, Int msspw, Double time, Double freq) {
858 :
859 : if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::interpolate(...)" << endl;
860 :
861 0 : Bool newcal(false);
862 0 : Int scanOrObsInt=thisTimeSeg(msobs,msscan);
863 0 : IPosition ip(4,0,msspw,msfld,scanOrObsInt);
864 :
865 : // Loop over _output_ elements
866 0 : for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
867 : // Call fully _patched_ time-interpolator, keeping track of 'newness'
868 : // fills timeResult_/timeResFlag_ implicitly
869 0 : ip(0)=iMSElem;
870 :
871 0 : if (!tI_(ip)) {
872 : //if (iMSElem==0) cout << "Flagging: " << ip << endl;
873 :
874 : // casa log post
875 0 : if (!tIMissingLogged_(ip)) { // only if these coords not logged before
876 0 : LogIO log;
877 0 : ostringstream msg;
878 0 : String scanOrObsStr=(byScan_ ? "scan" : "obs");
879 0 : msg << "MS " << scanOrObsStr << "=" << scanOrObsInt;
880 0 : if (byField_) msg << " in fld=" << msfld;
881 0 : msg << ",spw=" << msspw
882 0 : << ",ant=" << iMSElem
883 0 : << " is selected for processing, but has no available calibration in " << ctname_
884 0 : << " as mapped, and will be flagged.";
885 0 : log << msg.str() << LogIO::WARN;
886 0 : tIMissingLogged_(ip)=true;
887 0 : }
888 0 : newcal=true;
889 : }
890 : else {
891 0 : if (freq>0.0)
892 0 : newcal|=tI_(ip)->interpolate(time,freq);
893 : else
894 0 : newcal|=tI_(ip)->interpolate(time);
895 : }
896 : }
897 :
898 : // Whole result referred to time result:
899 0 : result_(msspw,msfld,scanOrObsInt).reference(timeResult_(msspw,msfld,scanOrObsInt));
900 0 : resFlag_(msspw,msfld,scanOrObsInt).reference(timeResFlag_(msspw,msfld,scanOrObsInt));
901 :
902 : // Detect if obs or fld changed, and cal is obs- or fld-dep
903 0 : Bool diffobsfld(false);
904 0 : diffobsfld|=(byField_ && msfld!=lastFld_(msspw)); // field-dep, and field changed
905 0 : diffobsfld|=(byObs_ && msobs!=lastObs_(msspw)); // obs-dep, and obs changed
906 0 : diffobsfld|=(byScan_ && msscan!=lastScan_(msspw)); // scan-dep, and scan changed
907 0 : newcal|=diffobsfld; // update newcal for return
908 :
909 : // Remember for next pass
910 0 : lastFld_(msspw)=msfld;
911 0 : lastObs_(msspw)=msobs;
912 0 : lastScan_(msspw)=msscan;
913 :
914 0 : return newcal;
915 0 : }
916 :
917 0 : Bool CTPatchedInterp::interpolate(Int msobs, Int msscan, Int msfld, Int msspw, Double time, const Vector<Double>& freq) {
918 :
919 : if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::interpolate(...,freq)" << endl;
920 :
921 : // obsid non-degenerate only if byObs_
922 :
923 : // The number of requested channels
924 0 : uInt nMSChan=freq.nelements();
925 :
926 : // Ensure freq result Array is properly sized
927 0 : if (freqResult_(msspw,msfld,thisTimeSeg(msobs,msscan)).nelements()==0) {
928 0 : Int thisAltFld=altFld_(msfld);
929 0 : if (freqResult_(msspw,thisAltFld,thisTimeSeg(msobs,msscan)).nelements()==0) {
930 0 : freqResult_(msspw,thisAltFld,thisTimeSeg(msobs,msscan)).resize(nFPar_,nMSChan,nMSElem_);
931 0 : freqResFlag_(msspw,thisAltFld,thisTimeSeg(msobs,msscan)).resize(nPar_,nMSChan,nMSElem_);
932 0 : freqResFlag_(msspw,thisAltFld,thisTimeSeg(msobs,msscan)).set(true);
933 : }
934 0 : if (thisAltFld!=msfld) {
935 0 : freqResult_(msspw,msfld,thisTimeSeg(msobs,msscan)).reference(freqResult_(msspw,thisAltFld,thisTimeSeg(msobs,msscan)));
936 0 : freqResFlag_(msspw,msfld,thisTimeSeg(msobs,msscan)).reference(freqResFlag_(msspw,thisAltFld,thisTimeSeg(msobs,msscan)));
937 : }
938 : }
939 :
940 :
941 :
942 :
943 : {
944 :
945 : // CAS-11744
946 : // In the following, we verify that there are sufficient soln channels to
947 : // do the requested freq-dep interpolation. If there are not, we
948 : // change the freq-dep interp mode to one that will work (linear or
949 : // nearest, for nchan=2 or 1, respectively), and issue a logger warning.
950 : // This is necessary specifically for the current _flag_ interpolation code.
951 : // The warning is thought warranted since freq-dep calibration (B, Tsys)
952 : // usually has adequate nchan, and if it does not, this may be unexpected
953 : // (e.g., use of too large a freq solint).
954 :
955 : // Check freq sampling adequate for specified freq interpolation
956 0 : Int nSolChan=timeResult_(msspw,msfld,thisTimeSeg(msobs,msscan)).shape()(1);
957 0 : LogIO log;
958 0 : ostringstream msg;
959 :
960 0 : switch (freqInterpMethod0_) {
961 0 : case LINEAR: {
962 0 : if (nSolChan<2 && freqInterpMethodVec_(msspw)!=NEAREST) {
963 0 : freqInterpMethodVec_(msspw)=NEAREST; // change to nearest for this msspw
964 0 : ostringstream spwmapstr;
965 0 : spwmapstr << " (mapped to CT spw=" << spwMap_(msspw) << ")";
966 0 : msg << "In caltable "
967 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
968 0 : << " (" << ct_.keywordSet().asString("VisCal") << ")"
969 0 : << ":" << endl
970 0 : << " Insufficient solution channel sampling (nchan=" << nSolChan << ") for frequency-dependent LINEAR interpolation " << endl
971 0 : << " of calibration for MS spw=" << msspw
972 0 : << ( msspw!=spwMap_(msspw) ? String(spwmapstr) : "" )
973 0 : << "; using NEAREST instead.";
974 0 : log << LogOrigin("CTPatchedInterp","interpolate") << msg.str() << LogIO::WARN;
975 0 : }
976 0 : break;
977 : }
978 0 : case CUBIC:
979 : case SPLINE: {
980 0 : if (nSolChan<2 && freqInterpMethodVec_(msspw)>NEAREST) {
981 0 : freqInterpMethodVec_(msspw)=NEAREST; // change to nearest for this msspw
982 0 : ostringstream spwmapstr;
983 0 : spwmapstr << " (mapped to CT spw=" << spwMap_(msspw) << ")";
984 0 : msg << "In caltable "
985 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
986 0 : << " (" << ct_.keywordSet().asString("VisCal") << ")"
987 0 : << ":" << endl
988 0 : << " Insufficient solution channel sampling (nchan=" << nSolChan << ") for frequency-dependent CUBIC/SPLINE interpolation " << endl
989 0 : << " of calibration for MS spw=" << msspw
990 0 : << ( msspw!=spwMap_(msspw) ? String(spwmapstr) : "" )
991 0 : << "; using NEAREST instead.";
992 0 : log << LogOrigin("CTPatchedInterp","interpolate") << msg.str() << LogIO::WARN;
993 0 : }
994 0 : else if (nSolChan<4 && freqInterpMethodVec_(msspw)>LINEAR) {
995 0 : freqInterpMethodVec_(msspw)=LINEAR; // change to nearest for this msspw
996 0 : ostringstream spwmapstr;
997 0 : spwmapstr << " (mapped to CT spw=" << spwMap_(msspw) << ")";
998 0 : msg << "In caltable "
999 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
1000 0 : << " (" << ct_.keywordSet().asString("VisCal") << ")"
1001 0 : << ":" << endl
1002 0 : << " Insufficient solution channel sampling (nchan=" << nSolChan << ") for frequency-dependent CUBIC/SPLINE interpolation " << endl
1003 0 : << " of calibration for MS spw=" << msspw
1004 0 : << ( msspw!=spwMap_(msspw) ? String(spwmapstr) : "" )
1005 0 : << "; using LINEAR instead.";
1006 0 : log << LogOrigin("CTPatchedInterp","interpolate") << msg.str() << LogIO::WARN;
1007 0 : }
1008 0 : break;
1009 : }
1010 0 : default: {
1011 0 : break;
1012 : }
1013 : }
1014 0 : }
1015 :
1016 : // Use the msspw-specific one!
1017 0 : freqInterpMethod_=freqInterpMethodVec_(msspw);
1018 :
1019 0 : Bool newcal(false);
1020 0 : Int scanOrObsInt=thisTimeSeg(msobs,msscan);
1021 0 : IPosition ip(4,0,msspw,msfld,scanOrObsInt);
1022 :
1023 : // Loop over _output_ antennas
1024 0 : for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1025 : // Call time interpolation calculation; resample in freq if new
1026 : // (fills timeResult_/timeResFlag_ implicitly)
1027 0 : ip(0)=iMSElem;
1028 0 : if (!tI_(ip)) {
1029 : //if (iMSElem==0) cout << "Flagging: " << ip << endl;
1030 :
1031 : // casa log post
1032 0 : if (!tIMissingLogged_(ip)) { // only if these coords not logged before
1033 0 : LogIO log;
1034 0 : ostringstream msg;
1035 0 : String scanOrObsStr=(byScan_ ? "scan" : "obs");
1036 0 : msg << "MS " << scanOrObsStr << "=" << scanOrObsInt;
1037 0 : if (byField_) msg << " in fld=" << msfld;
1038 0 : msg << ",spw=" << msspw
1039 0 : << ",ant=" << iMSElem
1040 0 : << " is selected for processing, but has no available calibration in " << ctname_
1041 0 : << " as mapped, and will be flagged.";
1042 0 : log << msg.str() << LogIO::WARN;
1043 0 : tIMissingLogged_(ip)=true;
1044 0 : }
1045 :
1046 :
1047 0 : newcal=true;
1048 : }
1049 : else {
1050 :
1051 0 : if (tI_(ip)->interpolate(time)) {
1052 : // Resample in frequency
1053 0 : Matrix<Float> fR(freqResult_(msspw,msfld,scanOrObsInt).xyPlane(iMSElem));
1054 0 : Matrix<Bool> fRflg(freqResFlag_(msspw,msfld,scanOrObsInt).xyPlane(iMSElem));
1055 0 : Matrix<Float> tR(timeResult_(msspw,msfld,scanOrObsInt).xyPlane(iMSElem));
1056 0 : Matrix<Bool> tRflg(timeResFlag_(msspw,msfld,scanOrObsInt).xyPlane(iMSElem));
1057 0 : resampleInFreq(fR,fRflg,freq,tR,tRflg,freqIn_(spwMap_(msspw)));
1058 :
1059 : // Calibration is new
1060 0 : newcal=true;
1061 0 : }
1062 : }
1063 : }
1064 :
1065 : // Whole result referred to freq result:
1066 0 : result_(msspw,msfld,scanOrObsInt).reference(freqResult_(msspw,msfld,scanOrObsInt));
1067 0 : resFlag_(msspw,msfld,scanOrObsInt).reference(freqResFlag_(msspw,msfld,scanOrObsInt));
1068 :
1069 : // Detect if obs or fld changed, and cal is obs- or fld-dep
1070 0 : Bool diffobsfld(false);
1071 0 : diffobsfld|=(byField_ && msfld!=lastFld_(msspw)); // field-dep, and field changed
1072 0 : diffobsfld|=(byObs_ && msobs!=lastObs_(msspw)); // obs-dep, and obs changed
1073 0 : diffobsfld|=(byScan_ && msscan!=lastScan_(msspw)); // scan-dep, and scan changed
1074 0 : newcal|=diffobsfld; // update newcal for return
1075 :
1076 : /*
1077 : if (newcal) {
1078 : Double t0(86400.0*floor(time/86400.0));
1079 : cout << boolalpha
1080 : << "fld="<<msfld
1081 : << " obs="<<scanOrObsInt
1082 : << " spw="<<msspw
1083 : << " time="<< time-t0
1084 : << " diffobsfld=" << diffobsfld
1085 : << " new=" << newcal
1086 : << " tI_(ip)=" << tI_(ip)
1087 : << " chan="<< nMSChan/2
1088 : << " result=" << result_(msspw,msfld,scanOrObsInt)(0,nMSChan/2,0)
1089 : << " addr=" << &result_(msspw,msfld,scanOrObsInt)(0,nMSChan/2,0)
1090 : << endl;
1091 : }
1092 : */
1093 : // Remember for next pass
1094 0 : lastFld_(msspw)=msfld;
1095 0 : lastObs_(msspw)=msobs;
1096 0 : lastScan_(msspw)=msscan;
1097 :
1098 0 : return newcal;
1099 : }
1100 :
1101 : // spwOK info for users
1102 0 : Bool CTPatchedInterp::spwOK(Int spw) const {
1103 :
1104 0 : if (spw<Int(spwMap_.nelements()))
1105 0 : return this->spwInOK(spwMap_(spw));
1106 :
1107 : // Something wrong...
1108 0 : return false;
1109 :
1110 : }
1111 0 : Bool CTPatchedInterp::spwInOK(Int spw) const {
1112 :
1113 0 : if (spw<Int(spwInOK_.nelements()))
1114 0 : return spwInOK_(spw);
1115 :
1116 : // Something wrong
1117 0 : return false;
1118 :
1119 : }
1120 :
1121 :
1122 : // Report state
1123 0 : void CTPatchedInterp::state() {
1124 :
1125 : if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::state()" << endl;
1126 :
1127 0 : cout << "-state--------" << endl;
1128 0 : cout << " ct_ = " << Path(ct_.tableName()).baseName().before(".tempMemCal") << endl;
1129 0 : cout << boolalpha;
1130 0 : cout << " isCmplx_ = " << isCmplx_ << endl;
1131 0 : cout << " nPar_ = " << nPar_ << endl;
1132 0 : cout << " nFPar_ = " << nFPar_ << endl;
1133 0 : cout << " nMSTimeSeg_ = " << nMSTimeSeg_ << endl;
1134 0 : cout << " nMSFld_ = " << nMSFld_ << endl;
1135 0 : cout << " nMSSpw_ = " << nMSSpw_ << endl;
1136 0 : cout << " nMSAnt_ = " << nMSAnt_ << endl;
1137 0 : cout << " nMSElem_ = " << nMSElem_ << endl;
1138 0 : cout << " nCTTimeSeg_ = " << nCTTimeSeg_ << endl;
1139 0 : cout << " nCTFld_ = " << nCTFld_ << endl;
1140 0 : cout << " nCTSpw_ = " << nCTSpw_ << endl;
1141 0 : cout << " nCTAnt_ = " << nCTAnt_ << endl;
1142 0 : cout << " nCTElem_ = " << nCTElem_ << endl;
1143 0 : cout << " fldMap_ = " << fldMap_ << endl;
1144 0 : cout << " spwMap_ = " << spwMap_ << endl;
1145 0 : cout << " antMap_ = " << antMap_ << endl;
1146 0 : cout << " byObs_ = " << byObs_ << endl;
1147 0 : cout << " byScan_ = " << byScan_ << endl;
1148 0 : cout << " byField_ = " << byField_ << endl;
1149 0 : cout << " altFld_ = " << altFld_ << endl;
1150 0 : cout << " timeType_ = " << timeType_ << endl;
1151 0 : cout << " freqTypeStr_ = " << freqTypeStr_ << endl;
1152 0 : cout << " freqInterpMethod0_ = " << freqInterpMethod0_ << endl;
1153 0 : cout << " freqInterpMethodVec_ = " << freqInterpMethodVec_ << endl;
1154 0 : cout << " spwInOK_ = " << spwInOK_ << endl;
1155 0 : }
1156 :
1157 0 : void CTPatchedInterp::sliceTable() {
1158 :
1159 : if (CTPATCHEDINTERPVERB) cout << " CTPatchedInterp::sliceTable()" << endl;
1160 :
1161 : // This method generates per-spw, per-antenna (and eventually per-field?)
1162 : // caltables.
1163 :
1164 : // Ensure time sort of input table
1165 : // TBD (here or inside loop?)
1166 :
1167 : // Indexed by the fields, spws, ants in the cal table (pre-mapped)
1168 0 : ctSlices_.resize(IPosition(4,nCTElem_,nCTSpw_,nCTFld_,nCTTimeSeg_));
1169 0 : ctSlices_.set(NULL);
1170 :
1171 : // Initialize spwInOK_
1172 0 : spwInOK_.resize(nCTSpw_);
1173 0 : spwInOK_.set(false);
1174 :
1175 : // Set up iterator
1176 : // TBD: handle baseline-based case!
1177 0 : Block<String> sortcol;
1178 0 : Int addobs( (byObs_ ? 1 : 0) ); // slicing by obs?
1179 0 : Int addscan( (byScan_ ? 1 : 0) ); // slicing by scan?
1180 0 : Int addfld( (byField_ ? 1 : 0) ); // slicing by field?
1181 :
1182 0 : switch(mtype_) {
1183 0 : case VisCalEnum::GLOBAL: {
1184 :
1185 0 : throw(AipsError("CTPatchedInterp::sliceTable: No non-Mueller/Jones yet."));
1186 :
1187 : sortcol.resize(1+addobs+addscan+addfld);
1188 : if (byObs_) sortcol[0]="OBSERVATION_ID"; // slicing by obs
1189 : if (byScan_) sortcol[0+addobs]="SCAN_NUMBER"; // slicing by scan
1190 : if (byField_) sortcol[0+addobs+addscan]="FIELD_ID"; // slicing by field
1191 : sortcol[0+addobs+addscan+addfld]="SPECTRAL_WINDOW_ID";
1192 : ROCTIter ctiter(ct_,sortcol);
1193 : while (!ctiter.pastEnd()) {
1194 : Int ispw=ctiter.thisSpw();
1195 : Int ifld = (byField_ ? ctiter.thisField() : 0); // use 0 if not slicing by field
1196 : Int iobs = (byObs_ ? ctiter.thisObs() : 0); // use 0 if not slicing by obs
1197 : iobs = (byScan_ ? ctiter.thisScan() : iobs); // XXX perscan overrides perobs?
1198 : IPosition ip(4,0,ispw,ifld,iobs);
1199 : ctSlices_(ip)= new NewCalTable(ctiter.table());
1200 : spwInOK_(ispw)=(spwInOK_(ispw) || ctSlices_(ip)->nrow()>0);
1201 : ctiter.next();
1202 : }
1203 : break;
1204 : }
1205 0 : case VisCalEnum::MUELLER: {
1206 0 : sortcol.resize(3+addobs+addscan+addfld);
1207 0 : if (byObs_) sortcol[0]="OBSERVATION_ID"; // slicing by obs
1208 0 : if (byScan_) sortcol[0+addobs]="SCAN_NUMBER"; // slicing by scan
1209 0 : if (byField_) sortcol[0+addobs+addscan]="FIELD_ID"; // slicing by field
1210 0 : sortcol[0+addobs+addscan+addfld]="SPECTRAL_WINDOW_ID";
1211 0 : sortcol[1+addobs+addscan+addfld]="ANTENNA1";
1212 0 : sortcol[2+addobs+addscan+addfld]="ANTENNA2";
1213 0 : ROCTIter ctiter(ct_,sortcol);
1214 0 : while (!ctiter.pastEnd()) {
1215 0 : Int ispw=ctiter.thisSpw();
1216 0 : Int iant1=ctiter.thisAntenna1();
1217 0 : Int iant2=ctiter.thisAntenna2();
1218 0 : Int ibln=blnidx(iant1,iant2,nCTAnt_);
1219 0 : Int ifld = (byField_ ? ctiter.thisField() : 0); // use 0 if not slicing by field
1220 0 : Int iobs = (byObs_ ? ctiter.thisObs() : 0); // use 0 if not slicing by obs
1221 0 : iobs = (byScan_ ? ctiter.thisScan() : iobs); // XXX perscan overrides perobs?
1222 0 : IPosition ip(4,ibln,ispw,ifld,iobs);
1223 0 : ctSlices_(ip)=new NewCalTable(ctiter.table());
1224 0 : spwInOK_(ispw)=(spwInOK_(ispw) || ctSlices_(ip)->nrow()>0);
1225 0 : ctiter.next();
1226 0 : }
1227 0 : break;
1228 0 : }
1229 0 : case VisCalEnum::JONES: {
1230 0 : sortcol.resize(2+addobs+addscan+addfld);
1231 0 : if (byObs_) sortcol[0]="OBSERVATION_ID"; // slicing by obs
1232 0 : if (byScan_) sortcol[0+addobs]="SCAN_NUMBER"; // slicing by scan
1233 0 : if (byField_) sortcol[0+addobs+addscan]="FIELD_ID"; // slicing by field
1234 0 : sortcol[0+addobs+addscan+addfld]="SPECTRAL_WINDOW_ID";
1235 0 : sortcol[1+addobs+addscan+addfld]="ANTENNA1";
1236 0 : ROCTIter ctiter(ct_,sortcol);
1237 0 : while (!ctiter.pastEnd()) {
1238 0 : Int ispw=ctiter.thisSpw();
1239 0 : Int iant=ctiter.thisAntenna1();
1240 0 : Int ifld = (byField_ ? ctiter.thisField() : 0); // use 0 if not slicing by field
1241 0 : Int iobs = (byObs_ ? ctiter.thisObs() : 0); // use 0 if not slicing by obs
1242 0 : iobs = (byScan_ ? ctiter.thisScan() : iobs); // XXX perscan overrides perobs?
1243 0 : IPosition ip(4,iant,ispw,ifld,iobs);
1244 0 : ctSlices_(ip)= new NewCalTable(ctiter.table());
1245 0 : spwInOK_(ispw)=(spwInOK_(ispw) || ctSlices_(ip)->nrow()>0);
1246 0 : ctiter.next();
1247 0 : }
1248 0 : break;
1249 0 : }
1250 : }
1251 :
1252 0 : }
1253 :
1254 : // Initialize by iterating over the supplied table
1255 0 : void CTPatchedInterp::makeInterpolators() {
1256 :
1257 : if (CTPATCHEDINTERPVERB) cout << " CTPatchedInterp::initialize()" << endl;
1258 :
1259 : // Save caltable name for log messages
1260 0 : ctname_=Path(ct_.antenna().tableName().before(".tempMem")).baseName();
1261 :
1262 : // Size/initialize interpolation engines
1263 0 : IPosition tIsize(4,nMSElem_,nMSSpw_,nMSFld_,nMSTimeSeg_);
1264 0 : tI_.resize(tIsize);
1265 0 : tI_.set(NULL);
1266 0 : tIdel_.resize(tIsize);
1267 0 : tIdel_.set(false);
1268 0 : tIMissingLogged_.resize(tIsize);
1269 0 : tIMissingLogged_.set(false);
1270 :
1271 0 : Bool reportBadSpw(false);
1272 0 : for (Int iMSTimeSeg=0;iMSTimeSeg<nMSTimeSeg_;++iMSTimeSeg) {
1273 0 : for (Int iMSFld=0;iMSFld<nMSFld_;++iMSFld) {
1274 :
1275 0 : if (altFld_(iMSFld)==iMSFld) { // i.e., this field has its own solutions
1276 :
1277 : // cout << "Making interpolators for " << iMSFld << " (mapped from " << fldMap_(iMSFld) << ")" << endl;
1278 :
1279 0 : std::set<uInt> spws;
1280 0 : if (msmc_)
1281 0 : spws = msmc_->msmd().getSpwsForField(iMSFld);
1282 :
1283 0 : for (Int iMSSpw=0;iMSSpw<nMSSpw_;++iMSSpw) {
1284 :
1285 : // Only if the required CT spw is available
1286 : // (spwmap applied in spwOK method)
1287 0 : if (this->spwOK(iMSSpw)) {
1288 :
1289 : // Size up the timeResult_ Cube (NB: channel shape matches Cal Table)
1290 0 : if (timeResult_(iMSSpw,iMSFld,iMSTimeSeg).nelements()==0) {
1291 0 : timeResult_(iMSSpw,iMSFld,iMSTimeSeg).resize(nFPar_,nChanIn_(spwMap_(iMSSpw)),nMSElem_);
1292 0 : timeResFlag_(iMSSpw,iMSFld,iMSTimeSeg).resize(nPar_,nChanIn_(spwMap_(iMSSpw)),nMSElem_);
1293 0 : timeResFlag_(iMSSpw,iMSFld,iMSTimeSeg).set(true);
1294 : }
1295 0 : for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1296 : // Realize the mapping
1297 0 : IPosition ictip(4,elemMap_(iMSElem),spwMap_(iMSSpw),fldMap_(iMSFld),iMSTimeSeg);
1298 0 : IPosition tIip(4,iMSElem,iMSSpw,iMSFld,iMSTimeSeg);
1299 0 : Matrix<Float> tR(timeResult_(iMSSpw,iMSFld,iMSTimeSeg).xyPlane(iMSElem));
1300 0 : Matrix<Bool> tRf(timeResFlag_(iMSSpw,iMSFld,iMSTimeSeg).xyPlane(iMSElem));
1301 :
1302 : // If the ct slice exists, set up an interpolator
1303 0 : if (ictip(0) >= 0 && ictip(0) < nCTElem_ &&
1304 0 : ictip(1) >= 0 && ictip(1) < nCTSpw_ &&
1305 0 : ictip(2) >= 0 && ictip(2) < nCTFld_ &&
1306 0 : ictip(3) >= 0 && ictip(3) < nCTTimeSeg_ &&
1307 0 : ctSlices_(ictip)) {
1308 0 : NewCalTable& ict(*ctSlices_(ictip));
1309 0 : if (!ict.isNull()) {
1310 0 : tI_(tIip)=(*cttifactoryptr_)(ict,timeType_,tR,tRf);
1311 0 : tIdel_(tIip)=true;
1312 : }
1313 : }
1314 : else {
1315 : // the required ct slice is empty, so arrange to flag it
1316 0 : tI_(tIip)=NULL;
1317 0 : tR.set(0.0);
1318 0 : tRf.set(true);
1319 : }
1320 0 : } // iMSElem
1321 : } // spwOK
1322 : else
1323 0 : reportBadSpw=true;
1324 : } // iMSSpw
1325 :
1326 0 : } // not re-using
1327 : else { // re-using
1328 : // Point to an existing interpolator group
1329 0 : Int thisAltFld=altFld_(iMSFld);
1330 :
1331 : // cout << "Reusing interpolators from " << thisAltFld << " for " << iMSFld << " (mapped to " << fldMap_(iMSFld) << ")" << endl;
1332 :
1333 0 : for (Int iMSSpw=0;iMSSpw<nMSSpw_;++iMSSpw) {
1334 0 : timeResult_(iMSSpw,iMSFld,iMSTimeSeg).reference(timeResult_(iMSSpw,thisAltFld,iMSTimeSeg));
1335 0 : timeResFlag_(iMSSpw,iMSFld,iMSTimeSeg).reference(timeResFlag_(iMSSpw,thisAltFld,iMSTimeSeg));
1336 0 : for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1337 0 : IPosition tIip0(4,iMSElem,iMSSpw,iMSFld,iMSTimeSeg),tIip1(4,iMSElem,iMSSpw,thisAltFld,iMSTimeSeg);
1338 0 : tI_(tIip0)=tI_(tIip1);
1339 : // if (!tI_(tIip0) && iMSElem==0) cout << "ouch---------------------" << "iMSTimeSeg="<<iMSTimeSeg<<" iMSFld="<<iMSFld<<" spw="<< iMSSpw << " ant="<<iMSElem<< endl;
1340 0 : }
1341 : }
1342 : }
1343 : } // iMSFld
1344 : } // iMSTimeSeg
1345 :
1346 :
1347 0 : if (reportBadSpw) {
1348 0 : cout << "The following MS spws have no corresponding cal spws in " << ctname_ << ": ";
1349 0 : for (Int iMSSpw=0;iMSSpw<nMSSpw_;++iMSSpw)
1350 : // (spwmap applied in spwOK method)
1351 0 : if (!this->spwOK(iMSSpw)) cout << iMSSpw << " ";
1352 0 : cout << endl;
1353 : }
1354 :
1355 0 : }
1356 :
1357 :
1358 :
1359 0 : void CTPatchedInterp::setFldMap(const MSField& msfld) {
1360 :
1361 0 : MSFieldColumns fcol(msfld);
1362 0 : setFldMap(fcol);
1363 :
1364 0 : }
1365 :
1366 0 : void CTPatchedInterp::setFldMap(const MSFieldColumns& fcol) {
1367 :
1368 : // Set the default fldmap
1369 0 : setDefFldMap();
1370 : // cout << "Nominal fldMap_ = " << fldMap_ << endl;
1371 :
1372 0 : ROCTColumns ctcol(ct_);
1373 :
1374 : // Discern _available_ fields in the CT
1375 0 : Vector<Int> ctFlds;
1376 0 : ctcol.fieldId().getColumn(ctFlds);
1377 0 : Int nAvFlds=genSort(ctFlds,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
1378 0 : ctFlds.resize(nAvFlds,true);
1379 :
1380 : //cout << "nAvFlds = " << nAvFlds << endl;
1381 : //cout << "ctFlds = " << ctFlds << endl;
1382 :
1383 : // If only one CT field, just use it
1384 :
1385 0 : Vector<Double> fseps(0);
1386 0 : if (nAvFlds==1)
1387 0 : fldMap_.set(ctFlds(0));
1388 : else {
1389 : // For each MS field, find the nearest available CT field
1390 0 : Int nMSFlds=fcol.nrow();
1391 0 : fseps.resize(nMSFlds);
1392 0 : fseps.set(0.0);
1393 0 : MDirection msdir,ctdir;
1394 0 : Vector<Double> sep(nAvFlds);
1395 0 : IPosition ipos(1,0); // get the first direction stored (no poly yet)
1396 0 : for (Int iMSFld=0;iMSFld<nMSFlds;++iMSFld) {
1397 0 : msdir=fcol.phaseDirMeas(iMSFld);
1398 : //cout << iMSFld << ":" << msdir.getValue() << endl;
1399 0 : sep.set(DBL_MAX);
1400 0 : for (Int iCTFld=0;iCTFld<nAvFlds;++iCTFld) {
1401 : // Get cal field direction, converted to ms field frame
1402 0 : ctdir=ctcol.field().phaseDirMeas(ctFlds(iCTFld));
1403 0 : MDirection::Convert(ctdir,msdir.getRef());
1404 : //cout << " c:" << ctFlds(iCTFld) << ":" << ctdir.getValue() << endl;
1405 0 : sep(iCTFld)=ctdir.getValue().separation(msdir.getValue());
1406 : }
1407 : // Sort separations
1408 0 : Vector<uInt> ord;
1409 0 : Int nsep=genSort(ord,sep,Sort::Ascending,(Sort::QuickSort | Sort::Ascending));
1410 :
1411 : /*
1412 : cout << " ord=" << ord << endl;
1413 : cout << " nsep=" << nsep << endl;
1414 : cout << " sep=" << sep << " " << sep*(180.0/C::pi)<< endl;
1415 : */
1416 :
1417 : // Trap case of duplication of nearest separation
1418 0 : if (nsep>1 && sep(ord(1))==sep(ord(0)))
1419 0 : throw(AipsError("Found more than one field at minimum distance, can't decide!"));
1420 :
1421 0 : fseps(iMSFld)=sep(ord(0));
1422 0 : fldMap_(iMSFld)=ctFlds(ord(0));
1423 0 : }
1424 0 : }
1425 0 : fseps*=(180.0/C::pi);
1426 0 : LogIO log;
1427 0 : ostringstream msg;
1428 0 : msg << "Calibration field mapping for "
1429 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
1430 0 : << " (via gainfield='nearest'): "
1431 0 : << fldMap_ << endl
1432 0 : << " Separations (deg): " << fseps;
1433 0 : log << msg.str() << LogIO::POST;
1434 :
1435 : // cout << ct_.tableName() << ": fldMap_ = " << fldMap_ << endl;
1436 0 : }
1437 :
1438 0 : void CTPatchedInterp::setFldMap(Vector<Int>& fldmap) {
1439 :
1440 : // Set the default spwmap first, then we'll ammend it
1441 0 : setDefFldMap();
1442 :
1443 0 : Int nfld=fldmap.nelements();
1444 :
1445 : // Must specify no more than needed, but at least one
1446 0 : AlwaysAssert(nfld>0,AipsError);
1447 0 : AlwaysAssert(nfld<=nMSFld_,AipsError);
1448 :
1449 : // Discern _available_ fields in the CT
1450 0 : ROCTColumns ctcol(ct_);
1451 0 : Vector<Int> ctFlds;
1452 0 : ctcol.fieldId().getColumn(ctFlds);
1453 0 : Int nAvFlds=genSort(ctFlds,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
1454 0 : ctFlds.resize(nAvFlds,true);
1455 :
1456 0 : for (Int i=0;i<nfld;++i) {
1457 0 : if (!anyEQ(ctFlds,fldmap(i)))
1458 0 : throw(AipsError("Specified fldmap contains an unavailable field: "+String(fldmap(i))));
1459 : else
1460 0 : fldMap_(i)=fldmap(i);
1461 : }
1462 0 : if (nfld<nMSFld_)
1463 : // Fill in the rest with last-specified
1464 0 : fldMap_(Slice(nfld,nMSFld_-nfld,1))=fldMap_(nfld-1);
1465 :
1466 :
1467 0 : LogIO log;
1468 0 : ostringstream msg;
1469 0 : msg << "Calibration field mapping for "
1470 0 : << Path(ct_.tableName()).baseName().before(".tempMemCal")
1471 0 : << " (via user specification): "
1472 0 : << fldMap_;
1473 0 : log << msg.str() << LogIO::POST;
1474 :
1475 0 : }
1476 :
1477 : // Calculate fldmap redundancy
1478 0 : void CTPatchedInterp::calcAltFld() {
1479 :
1480 0 : altFld_.resize(nMSFld_);
1481 0 : for (Int iMSFld=0;iMSFld<nMSFld_;++iMSFld) {
1482 0 : altFld_(iMSFld)=iMSFld; // nominally
1483 0 : for (Int ifld=0;ifld<iMSFld;++ifld)
1484 0 : if (fldMap_(ifld)==fldMap_(iMSFld))
1485 : //altFld_(iMSFld)=ifld;
1486 0 : altFld_(iMSFld)=altFld_(ifld);
1487 : }
1488 : /*
1489 : cout << "------------" << endl;
1490 : cout << "fldMap_ = " << fldMap_ << endl;
1491 : cout << "altFld_ = " << altFld_ << endl;
1492 : */
1493 0 : }
1494 :
1495 :
1496 :
1497 0 : void CTPatchedInterp::setSpwMap(Vector<Int>& spwmap) {
1498 :
1499 : // Set the default spwmap first, then we'll ammend it
1500 0 : setDefSpwMap();
1501 :
1502 0 : Int nspec=spwmap.nelements();
1503 :
1504 : // Do nothing, if nothing specified (and rely on default)
1505 0 : if (nspec==0) return;
1506 :
1507 : // Do nothing f a single -1 is specified
1508 0 : if (nspec==1 && spwmap(0)==-1) return;
1509 :
1510 :
1511 : // Alert user if too many spws specified
1512 0 : if (spwmap.nelements()>uInt(nMSSpw_))
1513 0 : throw(AipsError("Specified spwmap has more elements ("+String::toString(spwmap.nelements())+") than the number of spectral windows in the MS ("+String::toString(nMSSpw_)+")."));
1514 :
1515 : // Handle auto-fanout
1516 0 : if (spwmap(0)==-999) {
1517 : // Use first OK spw for all MS spws
1518 0 : Int gspw(0);
1519 0 : while (!spwInOK(gspw)) ++gspw;
1520 0 : spwMap_.set(gspw);
1521 : }
1522 : else {
1523 : // First trap out-of-range values
1524 0 : if (anyLT(spwmap,0))
1525 0 : throw(AipsError("Please specify positive indices in spwmap."));
1526 0 : if (anyGE(spwmap,nCTSpw_)) {
1527 0 : ostringstream o;
1528 0 : o << "Please specify spw indices <= maximum available ("
1529 0 : << (nCTSpw_-1) << " in " << Path(ct_.tableName()).baseName().before(".tempMemCal") << ")";
1530 0 : throw(AipsError(o.str()));
1531 0 : }
1532 :
1533 : // Now fill from spwmap
1534 0 : if (nspec==1)
1535 : // Use one value for all
1536 0 : spwMap_.set(spwmap(0));
1537 : else {
1538 : // set as many as are specified
1539 0 : IPosition blc(1,0);
1540 0 : IPosition trc(1,min(nspec-1,nMSSpw_-1));
1541 0 : spwMap_(blc,trc)=spwmap(blc,trc);
1542 0 : }
1543 : }
1544 :
1545 : //cout << "CTPatchedInterp::setSpwMap: Realized spwMap_ = " << spwMap_ << endl;
1546 :
1547 : }
1548 :
1549 :
1550 : // Resample in frequency
1551 0 : void CTPatchedInterp::resampleInFreq(Matrix<Float>& fres,Matrix<Bool>& fflg,const Vector<Double>& fout,
1552 : Matrix<Float>& tres,Matrix<Bool>& tflg,const Vector<Double>& fin) {
1553 :
1554 : if (CTPATCHEDINTERPVERB) cout << " CTPatchedInterp::resampleInFreq(...)" << endl;
1555 :
1556 : // if no good solutions coming in, return flagged
1557 0 : if (nfalse(tflg)==0) {
1558 0 : fflg.set(true);
1559 0 : return;
1560 : }
1561 :
1562 :
1563 0 : Int flparmod=nFPar_/nPar_; // for indexing the flag Matrices on the par axis
1564 :
1565 0 : Bool unWrapPhase=flparmod>1;
1566 :
1567 : // cout << "nFPar_,nPar_,flparmod = " << nFPar_ << "," << nPar_ << "," << flparmod << endl;
1568 :
1569 0 : fres=0.0;
1570 :
1571 0 : for (Int ifpar=0;ifpar<nFPar_;++ifpar) {
1572 :
1573 : // Slice by par (each Matrix row)
1574 0 : Vector<Float> fresi(fres.row(ifpar)), tresi(tres.row(ifpar));
1575 0 : Vector<Bool> fflgi(fflg.row(ifpar/flparmod)), tflgi(tflg.row(ifpar/flparmod));
1576 :
1577 : // Mask time result by flags
1578 0 : Vector<Double> mfin=fin(!tflgi).getCompressedArray();
1579 :
1580 0 : if (mfin.nelements()==0) {
1581 : // cout << ifpar << " All chans flagged!" << endl;
1582 : // Everything flagged this par
1583 : // Just flag, zero and go on to the next one
1584 0 : fflgi.set(true);
1585 0 : fresi.set(0.0);
1586 0 : continue;
1587 : }
1588 :
1589 0 : mfin/=1.0e9; // in GHz
1590 0 : Vector<Float> mtresi=tresi(!tflgi).getCompressedArray();
1591 :
1592 : // Trap case of same in/out frequencies
1593 0 : if (fout.nelements()==mfin.nelements() && allNear(fout,mfin,1.e-10)) {
1594 : // Just copy
1595 0 : fresi=mtresi;
1596 0 : fflgi.set(false); // none are flagged
1597 0 : continue;
1598 : }
1599 :
1600 0 : if (ifpar%2==1 && unWrapPhase) {
1601 0 : for (uInt i=1;i<mtresi.nelements();++i) {
1602 0 : while ( (mtresi(i)-mtresi(i-1))>C::pi ) mtresi(i)-=C::_2pi;
1603 0 : while ( (mtresi(i)-mtresi(i-1))<-C::pi ) mtresi(i)+=C::_2pi;
1604 : }
1605 : }
1606 :
1607 : // Set flags carefully
1608 0 : resampleFlagsInFreq(fflgi,fout,tflgi,fin);
1609 :
1610 :
1611 : // Always use nearest on edges
1612 : // TBD: trap cases where frequencies don't overlap at all
1613 : // (fout(hi)<mfin(0) || fout(lo)> mfin(ihi))..... already OK (lo>hi)?
1614 : // TBD: optimize the following by forming Slices in the
1615 : // while loops and doing Array assignment once afterwords
1616 :
1617 0 : Int nfreq=fout.nelements();
1618 0 : Int lo=0;
1619 0 : Int hi=fresi.nelements()-1;
1620 0 : Double inlo(mfin(0));
1621 0 : Int ihi=mtresi.nelements()-1;
1622 0 : Double inhi(mfin(ihi));
1623 :
1624 : // Handle 'nearest' extrapolation in sideband-dep way
1625 0 : Bool inUSB(inhi>inlo);
1626 0 : Bool outUSB(fout(hi)>fout(lo));
1627 0 : if (inUSB) {
1628 0 : if (outUSB) {
1629 0 : while (lo<nfreq && fout(lo)<=inlo) fresi(lo++)=mtresi(0);
1630 0 : while (hi>-1 && fout(hi)>=inhi) fresi(hi--)=mtresi(ihi);
1631 : }
1632 : else { // "outLSB"
1633 0 : while (lo<nfreq && fout(lo)>=inhi) fresi(lo++)=mtresi(ihi);
1634 0 : while (hi>-1 && fout(hi)<=inlo) fresi(hi--)=mtresi(0);
1635 : }
1636 : }
1637 : else { // "inLSB"
1638 0 : if (outUSB) {
1639 0 : while (lo<nfreq && fout(lo)<=inhi) fresi(lo++)=mtresi(ihi);
1640 0 : while (hi>-1 && fout(hi)>=inlo) fresi(hi--)=mtresi(0);
1641 : }
1642 : else { // "outLSB"
1643 0 : while (lo<nfreq && fout(lo)>=inlo) fresi(lo++)=mtresi(0);
1644 0 : while (hi>-1 && fout(hi)<=inhi) fresi(hi--)=mtresi(ihi);
1645 : }
1646 : }
1647 :
1648 : // cout << "lo, hi = " << lo << ","<<hi << endl;
1649 :
1650 0 : if (lo>hi) continue; // Frequencies didn't overlap, nearest was used
1651 :
1652 : // Use InterpolateArray1D to fill in the middle
1653 0 : IPosition blc(1,lo), trc(1,hi);
1654 0 : Vector<Float> slfresi(fresi(blc,trc));
1655 0 : Vector<Double> slfout(fout(blc,trc));
1656 :
1657 0 : InterpolateArray1D<Double,Float>::interpolate(slfresi,slfout,mfin,mtresi,freqInterpMethod_);
1658 :
1659 0 : }
1660 : }
1661 :
1662 0 : void CTPatchedInterp::resampleFlagsInFreq(Vector<Bool>& flgout,const Vector<Double>& fout,
1663 : Vector<Bool>& flgin,const Vector<Double>& fin) {
1664 :
1665 : // cout << "resampleFlagsInFreq" << endl;
1666 :
1667 0 : Vector<Double> finGHz=fin/1e9;
1668 :
1669 : // Handle chan-dep flags
1670 0 : if (freqTypeStr_.contains("flag")) {
1671 :
1672 : // Determine implied mode-dep flags indexed by channel registration
1673 0 : uInt nflg=flgin.nelements();
1674 0 : Vector<Bool> flreg(nflg,false);
1675 0 : switch (freqInterpMethod_) {
1676 0 : case NEAREST: {
1677 : // Just use input flags
1678 0 : flreg.reference(flgin);
1679 0 : break;
1680 : }
1681 0 : case LINEAR: {
1682 0 : for (uInt i=0;i<nflg-1;++i)
1683 0 : flreg[i]=(flgin[i] || flgin[i+1]);
1684 0 : flreg[nflg-1]=flreg[nflg-2];
1685 0 : break;
1686 : }
1687 0 : case CUBIC:
1688 : case SPLINE: {
1689 0 : for (uInt i=1;i<nflg-2;++i)
1690 0 : flreg[i]=(flgin[i-1] || flgin[i] || flgin[i+1] || flgin[i+2]);
1691 0 : flreg[0]=flreg[1];
1692 0 : flreg[nflg-2]=flreg[nflg-3];
1693 0 : flreg[nflg-1]=flreg[nflg-3];
1694 0 : break;
1695 : }
1696 : }
1697 :
1698 : // Now step through requested chans, setting flags
1699 0 : uInt ireg=0;
1700 0 : uInt nflgout=flgout.nelements();
1701 0 : for (uInt iflgout=0;iflgout<nflgout;++iflgout) {
1702 :
1703 : // Find nominal registration (the _index_ just left)
1704 0 : Bool exact(false);
1705 0 : ireg=binarySearch(exact,finGHz,fout(iflgout),nflg,0);
1706 :
1707 : // If registration is exact, assign verbatim
1708 : // NB: the calibration value calculation occurs agnostically w.r.t. flags,
1709 : // so the calculated value should also match
1710 : // TBD: Add "|| near(finGHz[ireg],fout(iflgout),1e-10) in case "exact" has
1711 : // precision issues?
1712 0 : if (exact) {
1713 0 : flgout[iflgout]=flgin[ireg];
1714 0 : continue;
1715 : }
1716 :
1717 : // Not exact, so carefully handle bracketing
1718 0 : if (ireg>0)
1719 0 : ireg-=1;
1720 0 : ireg=min(ireg,nflg-1);
1721 :
1722 : //while (finGHz(ireg)<=fout(iflgout) && ireg<nflg-1) {
1723 : // ireg+=1; // USB specific!
1724 : //}
1725 : //if (ireg>0 && finGHz(ireg)!=fout(iflgout)) --ireg; // registration is one sample prior
1726 :
1727 : // refine registration by interp type
1728 0 : switch (freqInterpMethod_) {
1729 0 : case NEAREST: {
1730 : // nearest might be forward sample
1731 0 : if ( ireg<(nflg-1) &&
1732 0 : abs(fout[iflgout]-finGHz[ireg])>abs(finGHz[ireg+1]-fout[iflgout]) )
1733 0 : ireg+=1;
1734 0 : break;
1735 : }
1736 0 : case LINEAR: {
1737 0 : if (ireg==(nflg-1)) // need one more sample to the right
1738 0 : ireg-=1;
1739 0 : break;
1740 : }
1741 0 : case CUBIC:
1742 : case SPLINE: {
1743 0 : if (ireg==0) ireg+=1; // need one more sample to the left
1744 0 : if (ireg>(nflg-3)) ireg=nflg-3; // need two more samples to the right
1745 0 : break;
1746 : }
1747 : }
1748 :
1749 : // Assign effective flag
1750 0 : flgout[iflgout]=flreg[ireg];
1751 :
1752 : /*
1753 : cout << iflgout << " "
1754 : << ireg << " "
1755 : << flreg[ireg]
1756 : << endl;
1757 : */
1758 :
1759 : }
1760 :
1761 0 : }
1762 : else
1763 : // We are interp/extrap-olating gaps absolutely
1764 0 : flgout.set(false);
1765 :
1766 0 : }
1767 :
1768 0 : void CTPatchedInterp::setElemMap() {
1769 :
1770 : // Ensure the antMap_ is set
1771 0 : if (antMap_.nelements()!=uInt(nMSAnt_))
1772 0 : setDefAntMap();
1773 :
1774 : // Handle cases
1775 0 : switch(mtype_) {
1776 0 : case VisCalEnum::GLOBAL: {
1777 :
1778 0 : throw(AipsError("CTPatchedInterp::sliceTable: No non-Mueller/Jones yet."));
1779 :
1780 : // There is only 1
1781 : AlwaysAssert(nMSElem_==1,AipsError);
1782 : elemMap_.resize(nMSElem_);
1783 : elemMap_.set(0);
1784 :
1785 : break;
1786 : }
1787 0 : case VisCalEnum::MUELLER: {
1788 0 : elemMap_.resize(nMSElem_);
1789 0 : conjTab_.resize(nMSElem_);
1790 0 : conjTab_.set(false);
1791 0 : Int iMSElem(0),a1in(0),a2in(0);
1792 0 : for (Int iMSAnt=0;iMSAnt<nMSAnt_;++iMSAnt) {
1793 0 : a1in=antMap_(iMSAnt);
1794 0 : for (Int jAntOut=iMSAnt;jAntOut<nMSAnt_;++jAntOut) {
1795 0 : a2in=antMap_(jAntOut);
1796 0 : if (a1in<=a2in)
1797 0 : elemMap_(iMSElem)=blnidx(a1in,a2in,nMSAnt_);
1798 : else {
1799 0 : elemMap_(iMSElem)=blnidx(a2in,a1in,nMSAnt_);
1800 0 : conjTab_(iMSElem)=true; // we must conjugate Complex params!
1801 : }
1802 0 : ++iMSElem;
1803 : } // jAntOut
1804 : } // iMSAnt
1805 0 : break;
1806 : }
1807 0 : case VisCalEnum::JONES: {
1808 : // Just reference the antMap_
1809 0 : elemMap_.reference(antMap_);
1810 0 : break;
1811 : }
1812 : } // switch
1813 0 : }
1814 :
1815 :
1816 0 : InterpolateArray1D<Double,Float>::InterpolationMethod CTPatchedInterp::ftype(String& strtype) {
1817 0 : if (strtype.contains("nearest"))
1818 0 : return InterpolateArray1D<Double,Float>::nearestNeighbour;
1819 0 : if (strtype.contains("linear"))
1820 0 : return InterpolateArray1D<Double,Float>::linear;
1821 0 : if (strtype.contains("cubic"))
1822 0 : return InterpolateArray1D<Double,Float>::cubic;
1823 0 : if (strtype.contains("spline"))
1824 0 : return InterpolateArray1D<Double,Float>::spline;
1825 :
1826 : // cout << "Using linear for freq interpolation as last resort." << endl;
1827 0 : return InterpolateArray1D<Double,Float>::linear;
1828 :
1829 :
1830 : }
1831 :
1832 :
1833 : } //# NAMESPACE CASA - END
|
