Line data Source code
1 : //# CTPatchPanel.cc: Implementation of CTPatchPanel.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/CLPatchPanel.h>
29 : #include <synthesis/CalTables/CTInterface.h>
30 : #include <synthesis/CalTables/CTIter.h>
31 :
32 : #include <casacore/scimath/Mathematics/InterpolateArray1D.h>
33 : #include <casacore/casa/Utilities/GenSort.h>
34 : #include <casacore/casa/OS/Path.h>
35 :
36 : #include <casacore/ms/MSSel/MSSelectableTable.h>
37 : #include <casacore/ms/MSSel/MSSelection.h>
38 : #include <casacore/ms/MSSel/MSSelectionTools.h>
39 :
40 : #include <casacore/casa/aips.h>
41 :
42 : #define CTPATCHPANELVERB false
43 :
44 : //#include <casa/BasicSL/Constants.h>
45 : //#include <casa/OS/File.h>
46 : #include <casacore/casa/Logging/LogMessage.h>
47 : #include <casacore/casa/Logging/LogSink.h>
48 : #include <casacore/casa/Logging/LogIO.h>
49 :
50 : using namespace casacore;
51 : namespace casa { //# NAMESPACE CASA - BEGIN
52 :
53 0 : CalPatchKey::CalPatchKey(IPosition keyids) :
54 0 : cpk_(keyids.asVector())
55 0 : {}
56 :
57 : // Lexographical lessthan
58 0 : Bool CalPatchKey::operator<(const CalPatchKey& other) const {
59 :
60 : // This method does a lexigraphical less-than, wherein negative
61 : // elements in either operand behave as equal (reflexively not
62 : // less than)
63 :
64 : // Loop over elements in precendence order
65 0 : for (Int i=0;i<6;++i) {
66 0 : if (cpk_[i]>-1 && other.cpk_[i]>-1 && cpk_[i]!=other.cpk_[i])
67 : // both non-negative and not equal, so evaluate element <
68 0 : return cpk_[i]<other.cpk_[i];
69 : }
70 : // All apparently equal so "<" is false
71 0 : return false;
72 : }
73 :
74 0 : MSCalPatchKey::MSCalPatchKey(Int obs,Int scan,Int fld,Int ent,Int spw,Int ant) :
75 0 : CalPatchKey(IPosition(6,obs,scan,fld,ent,spw,ant)),
76 0 : obs_(obs),scan_(scan),fld_(fld),ent_(ent),spw_(spw),ant_(ant)
77 0 : {}
78 :
79 : // text output
80 0 : String MSCalPatchKey::print() const {
81 0 : return "obs="+(obs_<0 ? "*" : String::toString(obs_))+" "
82 0 : "scan="+(scan_<0 ? "*" : String::toString(scan_))+" "
83 0 : "fld="+(fld_<0 ? "*" : String::toString(fld_))+" "
84 0 : "intent="+(ent_<0 ? "*" : String::toString(ent_))+" "
85 0 : "spw="+(spw_<0 ? "*" : String::toString(spw_))+" "
86 0 : "ant="+(ant_<0 ? "*" : String::toString(ant_));
87 : }
88 :
89 :
90 0 : CTCalPatchKey::CTCalPatchKey(Int clsl,Int obs,Int scan,Int fld,Int spw,Int ant) :
91 0 : CalPatchKey(IPosition(6,clsl,obs,scan,fld,spw,ant)),
92 0 : clsl_(clsl),obs_(obs),scan_(scan),fld_(fld),spw_(spw),ant_(ant)
93 0 : {}
94 :
95 : // text output
96 0 : String CTCalPatchKey::print() const {
97 0 : return "cl="+(clsl_<0 ? "*" : String::toString(clsl_))+" "
98 0 : "obs="+(obs_<0 ? "*" : String::toString(obs_))+" "
99 0 : "scan="+(scan_<0 ? "*" : String::toString(scan_))+" "
100 0 : "fld="+(fld_<0 ? "*" : String::toString(fld_))+" "
101 0 : "spw="+(spw_<0 ? "*" : String::toString(spw_))+" "
102 0 : "ant="+(ant_<0 ? "*" : String::toString(ant_));
103 : }
104 :
105 :
106 0 : CalMap::CalMap() :
107 0 : vcalmap_()
108 0 : {}
109 :
110 0 : CalMap::CalMap(const Vector<Int>& calmap) :
111 0 : vcalmap_(calmap)
112 : {
113 : // cout << "calmap addresses: " << calmap.data() << " " << vcalmap_.data()
114 : // << " nrefs= " << calmap.nrefs() << " " << vcalmap_.nrefs() << endl;
115 0 : }
116 :
117 0 : Int CalMap::operator()(Int msid) const {
118 : // TBD: reconsider algorithm (maybe just return msid if over-run?)
119 0 : Int ncalmap=vcalmap_.nelements();
120 0 : return (msid<ncalmap ? vcalmap_(msid) :
121 0 : (ncalmap>0 ? vcalmap_(ncalmap-1) : msid) ); // Avoid going off end
122 : }
123 :
124 0 : Vector<Int> CalMap::ctids(const Vector<Int>& msids) const {
125 0 : uInt ncalmap=vcalmap_.nelements();
126 :
127 : // If vector map unspecified, the calling context must work it out
128 : // (for obs, scan, fld, it means no specific mapping, just use all avaiable;
129 : // for spw, ant, it means [probably] use the same id)
130 0 : if (ncalmap<1 ||
131 0 : (ncalmap==1 && vcalmap_[0]<0))
132 0 : return Vector<Int>(1,-1);
133 :
134 0 : Vector<Bool> calmask(ncalmap,false);
135 0 : if (msids.nelements()==1 &&
136 0 : msids[0]<0)
137 : // MS ids indefinite, so all are ok
138 0 : calmask.set(true);
139 : else {
140 : // just do the ones that are requested
141 0 : for (uInt i=0;i<msids.nelements();++i) {
142 0 : const uInt& thismsid=msids(i);
143 0 : calmask(thismsid<ncalmap?thismsid:ncalmap-1)=true;
144 : }
145 : }
146 0 : Vector<Int> reqids;
147 0 : reqids=vcalmap_(calmask).getCompressedArray();
148 0 : Int nsort=genSort(reqids,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
149 0 : reqids.resize(nsort,true);
150 0 : return reqids;
151 0 : }
152 :
153 0 : Vector<Int> CalMap::msids(Int ctid,const Vector<Int>& superset) const {
154 :
155 0 : uInt ncalmap=vcalmap_.nelements();
156 :
157 : // If vector map unspecified, return [-1] which signals to calling
158 : // context that it must figure it out on its own (all or identity)
159 0 : if (ncalmap<1 ||
160 0 : (ncalmap==1 && vcalmap_[0]<0) )
161 0 : return Vector<Int>(1,-1);
162 :
163 : // mask vcalmap_ by specified ctid...
164 0 : Vector<Int> msidlist(ncalmap);
165 0 : indgen(msidlist);
166 0 : Vector<Bool> msidmask(vcalmap_==ctid);
167 : // ...and limit to superset, if specified
168 0 : if (superset.nelements()>0 && superset[0]>-1)
169 0 : for (Int i=0;i<Int(msidmask.nelements());++i)
170 0 : if (!anyEQ(superset,i))
171 0 : msidmask[i]=false; // exclude if not in superset
172 :
173 : // Return the ms id list
174 0 : return msidlist(msidmask).getCompressedArray();
175 0 : }
176 :
177 :
178 0 : FieldCalMap::FieldCalMap() :
179 : CalMap(),
180 0 : fieldcalmap_("")
181 0 : {}
182 :
183 0 : FieldCalMap::FieldCalMap(const Vector<Int>& calmap) :
184 : CalMap(calmap),
185 0 : fieldcalmap_("")
186 0 : {}
187 :
188 0 : FieldCalMap::FieldCalMap(const String fieldcalmap,
189 0 : const MeasurementSet& ms, const NewCalTable& ct) :
190 : CalMap(),
191 0 : fieldcalmap_(fieldcalmap)
192 : {
193 :
194 0 : if (fieldcalmap_=="nearest")
195 : // Calculate nearest map
196 0 : setNearestFieldMap(ms,ct);
197 : else
198 : // Attempt field selection
199 0 : setSelectedFieldMap(fieldcalmap,ms,ct);
200 :
201 0 : }
202 :
203 0 : FieldCalMap::FieldCalMap(const String fieldcalmap,
204 : const MeasurementSet& ms, const NewCalTable& ct,
205 0 : String& extfldsel) :
206 : CalMap(),
207 0 : fieldcalmap_(fieldcalmap)
208 : {
209 :
210 0 : if (fieldcalmap_=="nearest")
211 : // Calculate nearest map
212 0 : setNearestFieldMap(ms,ct);
213 : else
214 : // Attempt field selection
215 0 : setSelectedFieldMap(fieldcalmap,ms,ct,extfldsel);
216 :
217 0 : }
218 :
219 0 : void FieldCalMap::setNearestFieldMap(const MeasurementSet& ms, const NewCalTable& ct) {
220 : // Access MS and CT columns
221 0 : MSFieldColumns msfc(ms.field());
222 0 : ROCTColumns ctc(ct);
223 0 : setNearestFieldMap(msfc,ctc);
224 0 : }
225 0 : void FieldCalMap::setNearestFieldMap(const NewCalTable& ctasms, const NewCalTable& ct) {
226 : // Access MS and CT columns
227 0 : ROCTFieldColumns msfc(ctasms.field());
228 0 : ROCTColumns ctc(ct);
229 0 : setNearestFieldMap(msfc,ctc);
230 0 : }
231 :
232 0 : void FieldCalMap::setNearestFieldMap(const MSFieldColumns& msfc, const ROCTColumns& ctc) {
233 :
234 : // Nominally, this many field need a map
235 0 : Int nMSFlds=msfc.nrow();
236 0 : vcalmap_.resize(nMSFlds);
237 0 : vcalmap_.set(-1); // no map
238 :
239 : // Discern _available_ fields in the CT
240 0 : Vector<Int> ctFlds;
241 0 : ctc.fieldId().getColumn(ctFlds);
242 0 : Int nAvFlds=genSort(ctFlds,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
243 0 : ctFlds.resize(nAvFlds,true);
244 :
245 : // If only one CT field, just use it
246 0 : if (nAvFlds==1)
247 0 : vcalmap_.set(ctFlds(0));
248 : else {
249 : // For each MS field, find the nearest available CT field
250 0 : MDirection msdir,ctdir;
251 0 : Vector<Double> sep(nAvFlds);
252 0 : IPosition ipos(1,0); // get the first direction stored (no poly yet)
253 0 : for (Int iMSFld=0;iMSFld<nMSFlds;++iMSFld) {
254 0 : msdir=msfc.phaseDirMeas(iMSFld); // MS fld dir
255 0 : sep.set(DBL_MAX);
256 0 : for (Int iCTFld=0;iCTFld<nAvFlds;++iCTFld) {
257 : // Get cal field direction, converted to ms field frame
258 0 : ctdir=ctc.field().phaseDirMeas(ctFlds(iCTFld));
259 0 : MDirection::Convert(ctdir,msdir.getRef());
260 0 : sep(iCTFld)=ctdir.getValue().separation(msdir.getValue());
261 : }
262 : // Sort separations
263 0 : Vector<uInt> ord;
264 0 : Int nsep=genSort(ord,sep,Sort::Ascending,Sort::QuickSort);
265 :
266 : //cout << iMSFld << ":" << endl;
267 : //cout << " ord=" << ord << endl;
268 : //cout << " nsep=" << nsep << endl;
269 : //cout << " sep=" << sep << " " << sep*(180.0/C::pi)<< endl;
270 :
271 : // Trap case of duplication of nearest separation
272 0 : if (nsep>1 && sep(ord(1))==sep(ord(0)))
273 0 : throw(AipsError("Found more than one field at minimum distance, can't decide!"));
274 :
275 0 : vcalmap_(iMSFld)=ctFlds(ord(0));
276 0 : }
277 0 : }
278 : // cout << "FCM::setNearestFieldMap:*************** vcalmap_ = " << vcalmap_ << endl;
279 :
280 0 : }
281 :
282 0 : void FieldCalMap::setSelectedFieldMap(const String& fieldsel,
283 : const MeasurementSet& ms,
284 : const NewCalTable& ct) {
285 :
286 : // Nominally, this many fields need a map
287 0 : Int nMSFlds=ms.field().nrow();
288 0 : vcalmap_.resize(nMSFlds);
289 0 : vcalmap_.set(-1); // no map
290 :
291 0 : CTInterface cti(ct);
292 0 : MSSelection mss;
293 0 : mss.setFieldExpr(fieldsel);
294 0 : TableExprNode ten=mss.toTableExprNode(&cti);
295 0 : Vector<Int> fieldlist=mss.getFieldList();
296 :
297 0 : if (fieldlist.nelements()>1)
298 0 : throw(AipsError("Field mapping by selection can support only one field."));
299 :
300 0 : if (fieldlist.nelements()>0)
301 0 : vcalmap_.set(fieldlist[0]);
302 :
303 : // cout << "FCM::setSelectedFieldMap:*************** vcalmap_ = " << vcalmap_ << endl;
304 :
305 0 : }
306 :
307 0 : void FieldCalMap::setSelectedFieldMap(const String& fieldsel,
308 : const MeasurementSet& ms,
309 : const NewCalTable& ct,
310 : String& extfldsel) {
311 :
312 :
313 : try {
314 :
315 0 : CTInterface cti(ct);
316 0 : MSSelection mss;
317 0 : mss.setFieldExpr(fieldsel);
318 0 : TableExprNode ten=mss.toTableExprNode(&cti);
319 0 : Vector<Int> fieldlist=mss.getFieldList();
320 :
321 : // Nominally, don't apply selection externally
322 0 : extfldsel="";
323 :
324 0 : if (fieldlist.nelements()>1)
325 : // trigger external selection (multi-field case; see old gainfield)
326 0 : extfldsel=fieldsel;
327 0 : else if (fieldlist.nelements()==1) {
328 : // exactly 1 means fldmap=[f]*nfld
329 :
330 : // Nominally, this many fields need a map
331 0 : Int nMSFlds=ms.field().nrow();
332 0 : vcalmap_.resize(nMSFlds);
333 0 : vcalmap_.set(fieldlist[0]);
334 : }
335 : else
336 : // Field selection found no indices
337 0 : throw(AipsError(fieldsel+" matches no fields."));
338 :
339 0 : }
340 0 : catch ( AipsError err ) {
341 :
342 0 : throw(AipsError("Field mapping by selection failure: "+err.getMesg()));
343 :
344 0 : }
345 :
346 :
347 : // cout << "FCM::setSelectedFieldMap:*************** vcalmap_ = " << vcalmap_ << endl;
348 :
349 0 : return;
350 : }
351 :
352 :
353 0 : ObsCalMap::ObsCalMap() :
354 0 : CalMap()
355 0 : {}
356 :
357 0 : ObsCalMap::ObsCalMap(const String obscalmap, const MeasurementSet& ms) :
358 0 : CalMap()
359 : {
360 0 : if (obscalmap=="self") {
361 0 : vcalmap_.resize(ms.observation().nrow());
362 0 : indgen(vcalmap_);
363 : } else {
364 0 : throw(AipsError("Observation mapping failure: unrecognized keyword '" +
365 0 : obscalmap + "'."));
366 : }
367 0 : }
368 :
369 :
370 0 : ScanCalMap::ScanCalMap() :
371 0 : CalMap()
372 0 : {}
373 :
374 0 : ScanCalMap::ScanCalMap(const String scancalmap, const MeasurementSet& ms) :
375 0 : CalMap()
376 : {
377 0 : if (scancalmap=="self") {
378 0 : MSColumns msc(ms);
379 0 : Vector<Int> msScans;
380 0 : msc.scanNumber().getColumn(msScans);
381 0 : vcalmap_.resize(max(msScans) + 1);
382 0 : indgen(vcalmap_);
383 0 : } else {
384 0 : throw(AipsError("Scan mapping failure: unrecognized keyword '" +
385 0 : scancalmap + "'."));
386 : }
387 0 : }
388 :
389 :
390 0 : CalLibSlice::CalLibSlice(String obs, String scan, String fld, String ent, String spw,
391 : String tinterp,String finterp,
392 : Vector<Int> obsmap, Vector<Int> scanmap,
393 : Vector<Int> fldmap,
394 0 : Vector<Int> spwmap, Vector<Int> antmap) :
395 0 : obs(obs),scan(scan),fld(fld),ent(ent),spw(spw),
396 0 : tinterp(tinterp),finterp(finterp),
397 0 : obsmap(obsmap),scanmap(scanmap),fldmap(fldmap),
398 0 : spwmap(spwmap),antmap(antmap),extfldsel("")
399 0 : {}
400 :
401 : // Construct from a Record
402 0 : CalLibSlice::CalLibSlice(const Record& clslice,
403 : const MeasurementSet& ms,
404 0 : const NewCalTable& ct) :
405 0 : obs(),scan(),fld(),ent(),spw(),
406 0 : tinterp(),finterp(),
407 0 : obsmap(),scanmap(),fldmap(),spwmap(),antmap(),
408 0 : extfldsel("")
409 : {
410 :
411 0 : if (clslice.isDefined("obs")) {
412 0 : obs=clslice.asString("obs");
413 : }
414 0 : if (clslice.isDefined("scan")) {
415 0 : scan=clslice.asString("scan");
416 : }
417 0 : if (clslice.isDefined("field")) {
418 0 : fld=clslice.asString("field");
419 : }
420 0 : if (clslice.isDefined("intent")) {
421 0 : ent=clslice.asString("intent");
422 : }
423 0 : if (clslice.isDefined("spw")) {
424 0 : spw=clslice.asString("spw");
425 : }
426 :
427 0 : if (clslice.isDefined("tinterp")) {
428 0 : tinterp=clslice.asString("tinterp");
429 0 : if (tinterp=="") {
430 0 : tinterp="linear";
431 : }
432 : }
433 0 : if (clslice.isDefined("finterp")) {
434 0 : finterp=clslice.asString("finterp");
435 : }
436 :
437 0 : if (clslice.isDefined("obsmap")) {
438 : //cout << "obsmap.dataType() = " << clslice.dataType("obsmap") << endl;
439 0 : if (clslice.dataType("obsmap")==TpArrayInt)
440 0 : obsmap=CalMap(Vector<Int>(clslice.asArrayInt("obsmap")));
441 0 : if (clslice.dataType("obsmap")==TpString)
442 0 : obsmap=ObsCalMap(clslice.asString("obsmap"),ms);
443 : }
444 0 : if (clslice.isDefined("scanmap")) {
445 : //cout << "scanmap.dataType() = " << clslice.dataType("scanmap") << endl;
446 0 : if (clslice.dataType("scanmap")==TpArrayInt)
447 0 : scanmap=CalMap(Vector<Int>(clslice.asArrayInt("scanmap")));
448 0 : if (clslice.dataType("scanmap")==TpString)
449 0 : scanmap=ScanCalMap(clslice.asString("scanmap"),ms);
450 : }
451 0 : if (clslice.isDefined("fldmap")) {
452 0 : if (clslice.dataType("fldmap")==TpArrayInt)
453 0 : fldmap=FieldCalMap(clslice.asArrayInt("fldmap"));
454 0 : if (clslice.dataType("fldmap")==TpString)
455 0 : fldmap=FieldCalMap(clslice.asString("fldmap"),ms,ct,extfldsel);
456 : }
457 0 : if (clslice.isDefined("spwmap")) {
458 : // cout << "spwmap.dataType() = " << clslice.dataType("spwmap") << endl;
459 0 : if (clslice.dataType("spwmap")==TpArrayInt)
460 0 : spwmap=CalMap(clslice.asArrayInt("spwmap"));
461 : }
462 0 : if (clslice.isDefined("antmap")) {
463 : // cout << "antmap.dataType() = " << clslice.dataType("antmap") << endl;
464 0 : if (clslice.dataType("antmap")==TpArrayInt)
465 0 : antmap=CalMap(clslice.asArrayInt("antmap"));
466 : }
467 :
468 :
469 0 : }
470 :
471 : // Return contents as a Record
472 0 : Record CalLibSlice::asRecord() {
473 0 : Record rec;
474 0 : rec.define("obs",obs);
475 0 : rec.define("scan",scan);
476 0 : rec.define("field",fld);
477 0 : rec.define("intent",ent);
478 0 : rec.define("spw",spw);
479 :
480 0 : rec.define("tinterp",tinterp);
481 0 : rec.define("finterp",finterp);
482 :
483 0 : rec.define("obsmap",obsmap.vmap());
484 0 : rec.define("scanmap",scanmap.vmap());
485 0 : rec.define("fldmap",fldmap.vmap());
486 0 : rec.define("spwmap",spwmap.vmap());
487 0 : rec.define("antmap",antmap.vmap());
488 :
489 0 : return rec;
490 0 : }
491 :
492 0 : Bool CalLibSlice::validateCLS(const Record& clslice) {
493 :
494 0 : String missing("");
495 :
496 0 : if (!clslice.isDefined("obs"))
497 0 : missing+="obs ";
498 0 : if (!clslice.isDefined("scan"))
499 0 : missing+="scan ";
500 0 : if (!clslice.isDefined("field"))
501 0 : missing+="field ";
502 0 : if (!clslice.isDefined("intent"))
503 0 : missing+="intent ";
504 0 : if (!clslice.isDefined("spw"))
505 0 : missing+="spw ";
506 :
507 0 : if (!clslice.isDefined("tinterp"))
508 0 : missing+="tinterp ";
509 0 : if (!clslice.isDefined("finterp"))
510 0 : missing+="finterp ";
511 :
512 0 : if (!clslice.isDefined("obsmap"))
513 0 : missing+="obsmap ";
514 0 : if (!clslice.isDefined("scanmap"))
515 0 : missing+="scanmap ";
516 0 : if (!clslice.isDefined("fldmap"))
517 0 : missing+="fldmap ";
518 0 : if (!clslice.isDefined("spwmap"))
519 0 : missing+="swmap ";
520 0 : if (!clslice.isDefined("antmap"))
521 0 : missing+="antmap";
522 :
523 :
524 0 : if (missing.length()>0)
525 0 : throw(AipsError(missing));
526 :
527 : // Everything is ok if we get here
528 0 : return true;
529 0 : }
530 :
531 :
532 0 : String CalLibSlice::state() {
533 :
534 0 : ostringstream o;
535 :
536 0 : o << " MS: obs="+obs+" scan="+scan+" fld="+fld+" intent="+ent+" spw="+spw << endl
537 0 : << " CT: tinterp="+tinterp << " finterp="+finterp << endl
538 0 : << " obsmap=" << obsmap.vmap()
539 0 : << " scanmap=" << scanmap.vmap()
540 0 : << " fldmap=";
541 :
542 0 : if (extfldsel=="")
543 0 : o << fldmap.vmap();
544 : else
545 0 : o << extfldsel+" (by selection)";
546 :
547 0 : o << endl
548 0 : << " spwmap=" << spwmap.vmap()
549 0 : << " antmap=" << antmap.vmap()
550 0 : << endl;
551 :
552 0 : return String(o);
553 0 : }
554 :
555 :
556 0 : CLPPResult::CLPPResult() :
557 0 : result_(),resultFlag_()
558 0 : {}
559 :
560 0 : CLPPResult::CLPPResult(const IPosition& shape) :
561 0 : result_(shape,0.0),
562 0 : resultFlag_(shape,true)
563 0 : {}
564 0 : CLPPResult::CLPPResult(uInt nPar,uInt nFPar,uInt nChan,uInt nElem) :
565 0 : result_(nFPar,nChan,nElem,0.0),
566 0 : resultFlag_(nPar,nChan,nElem,true)
567 0 : {}
568 :
569 0 : CLPPResult& CLPPResult::operator=(const CLPPResult& other) {
570 : // Avoid Array deep copies!
571 0 : result_.reference(other.result_);
572 0 : resultFlag_.reference(other.resultFlag_);
573 0 : return *this;
574 : }
575 0 : void CLPPResult::resize(uInt nPar,uInt nFPar,uInt nChan,uInt nElem) {
576 0 : result_.resize(nFPar,nChan,nElem);
577 0 : result_.set(1.0);
578 0 : resultFlag_.resize(nPar,nChan,nElem);
579 0 : resultFlag_.set(true);
580 0 : }
581 :
582 : // (caltable only, mainly for testing)
583 0 : CLPatchPanel::CLPatchPanel(const String& ctname,
584 : const Record& callib,
585 : VisCalEnum::MatrixType mtype,
586 : Int nPar,
587 0 : const CTTIFactoryPtr cttifactoryptr) :
588 0 : ct_(NewCalTable::createCT(ctname,Table::Old,Table::Memory)),
589 0 : ctasms_(NewCalTable::createCT(ctname,Table::Old,Table::Memory)),
590 0 : ms_(),
591 0 : mtype_(mtype),
592 0 : isCmplx_(false),
593 0 : nPar_(nPar),
594 0 : nFPar_(nPar),
595 :
596 0 : nChanIn_(),
597 0 : freqIn_(),
598 0 : refFreqIn_(),
599 :
600 0 : nMSObs_(ct_.observation().nrow()),
601 0 : nMSFld_(ct_.field().nrow()),
602 0 : nMSSpw_(ct_.spectralWindow().nrow()),
603 0 : nMSAnt_(ct_.antenna().nrow()),
604 0 : nMSElem_(0), // set non-trivially in ctor body
605 0 : nCTObs_(ct_.observation().nrow()),
606 0 : nCTFld_(ct_.field().nrow()),
607 0 : nCTSpw_(ct_.spectralWindow().nrow()),
608 0 : nCTAnt_(ct_.antenna().nrow()),
609 0 : nCTElem_(0), // set non-trivially in ctor body
610 0 : lastresadd_(nMSSpw_,NULL),
611 0 : cttifactoryptr_(cttifactoryptr)
612 :
613 : // elemMap_(),
614 : // conjTab_(),
615 :
616 : {
617 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::CLPatchPanel(<no MS>)" << endl;
618 :
619 : // ia1dmethod_=ftype(freqType_);
620 : // cout << "ia1dmethod_ = " << ia1dmethod_ << endl;
621 :
622 0 : switch(mtype_) {
623 0 : case VisCalEnum::GLOBAL: {
624 :
625 0 : throw(AipsError("CLPatchPanel::ctor: No non-Mueller/Jones yet."));
626 :
627 : nCTElem_=1;
628 : nMSElem_=1;
629 : break;
630 : }
631 0 : case VisCalEnum::MUELLER: {
632 0 : nCTElem_=nCTAnt_*(nCTAnt_+1)/2;
633 0 : nMSElem_=nMSAnt_*(nMSAnt_+1)/2;
634 0 : break;
635 : }
636 0 : case VisCalEnum::JONES: {
637 0 : nCTElem_=nCTAnt_;
638 0 : nMSElem_=nMSAnt_;
639 0 : break;
640 : }
641 : }
642 :
643 : // How many _Float_ parameters?
644 0 : isCmplx_=(ct_.keywordSet().asString("ParType")=="Complex");
645 0 : if (isCmplx_) // Complex input
646 0 : nFPar_*=2; // interpolating 2X as many Float values
647 :
648 : // Set CT channel/freq info
649 0 : CTSpWindowColumns ctspw(ct_.spectralWindow());
650 0 : ctspw.numChan().getColumn(nChanIn_);
651 0 : freqIn_.resize(nCTSpw_);
652 0 : for (uInt iCTspw=0;iCTspw<ctspw.nrow();++iCTspw)
653 0 : ctspw.chanFreq().get(iCTspw,freqIn_(iCTspw),true);
654 0 : ctspw.refFrequency().getColumn(refFreqIn_);
655 :
656 : // Initialize maps
657 : // (ctids, msids, result arrays)
658 : // TBD
659 :
660 : //void CLPatchPanel::unpackCLinstance(CalLibSlice cls) {
661 :
662 : // Loop over callib slices
663 : // (a callib slice is one MS selection and interp setup for the present caltable)
664 0 : uInt nkey=callib.nfields();
665 0 : Int icls=-1;
666 0 : for (uInt ikey=0;ikey<nkey;++ikey) {
667 :
668 : // skip
669 0 : if (callib.type(ikey)!=TpRecord)
670 0 : continue;
671 :
672 0 : ++icls;
673 :
674 : // The current callib slice
675 0 : CalLibSlice cls(callib.asRecord(ikey));
676 :
677 0 : logsink_ << LogIO::NORMAL << ". " << icls << ":" << endl
678 0 : << cls.state() << LogIO::POST;
679 :
680 : // Apply callib instance MS selection to the "MS" (in this case it is a CT)
681 0 : NewCalTable clsms(ctasms_);
682 0 : this->selectOnCT(clsms,ctasms_,cls.obs,cls.scan,cls.fld,cls.spw,"");
683 :
684 : // Extract the MS indices we must calibrate
685 0 : Vector<Int> reqMSobs(1,-1),reqMSscan(1,-1),reqMSfld(1,-1);
686 0 : if (cls.obs.length()>0)
687 0 : reqMSobs.reference(this->getCLuniqueObsIds(clsms));
688 0 : if (cls.scan.length()>0)
689 0 : reqMSscan.reference(this->getCLuniqueScanIds(clsms));
690 0 : if (cls.fld.length()>0)
691 0 : reqMSfld.reference(this->getCLuniqueFldIds(clsms));
692 0 : Vector<Int> reqMSspw(this->getCLuniqueSpwIds(clsms));
693 0 : Vector<Int> reqMSant(nMSAnt_);
694 0 : indgen(reqMSant);
695 :
696 : // cout << "reqMSobs = " << reqMSobs << endl;
697 : // cout << "reqMSscan = " << reqMSscan << endl;
698 : // cout << "reqMSfld = " << reqMSfld << endl;
699 : // cout << "reqMSspw = " << reqMSspw << endl;
700 : // cout << "reqMSant = " << reqMSant << endl;
701 :
702 : // The intent list from the callib instance, to be used to index the msci_
703 0 : Vector<Int> theseMSint(1,-1); // Details TBD
704 :
705 : // cout << "theseMSint = " << theseMSint << endl;
706 :
707 : // WE DO TIME-ISH (OBS,SCAN,FLD) AXES IN OUTER LOOPS
708 :
709 : // The net CT obs required for the MS obs according to the obsmap
710 : // (in principle, this may contain CT obs ids that aren't available)
711 0 : Vector<Int> reqCTobs=cls.obsmap.ctids(reqMSobs);
712 :
713 : // cout << "reqCTobs = " << reqCTobs << endl;
714 :
715 : // For each required CT obs (and thus the MS obs ids requiring it)
716 0 : NewCalTable obsselCT(ct_);
717 0 : for (uInt iCTobs=0;iCTobs<reqCTobs.nelements();++iCTobs) {
718 0 : Int& thisCTobs=reqCTobs(iCTobs);
719 :
720 : // Apply thisCTobs selection to the CT
721 : // (only if a meaningful obsid is specified)
722 0 : if (thisCTobs!=-1)
723 0 : this->selectOnCT(obsselCT,ct_,String::toString(thisCTobs),"","","","");
724 :
725 : // The MS obss to be calibrated by thisCTobs (limited by the req superset)
726 : // (could be [-1], which means all)
727 0 : Vector<Int> theseMSobs=cls.obsmap.msids(thisCTobs,reqMSobs);
728 0 : if (theseMSobs.nelements()==1 && theseMSobs[0]<0)
729 0 : theseMSobs.reference(reqMSobs);
730 : // cout << " thisCTobs = " << thisCTobs << ": theseMSobs = " << theseMSobs << endl;
731 :
732 : // Apply theseMSobs selection to the MS
733 :
734 : // The net CT scan required for the MS scan according to the scanmap
735 : // (in principle, this may contain CT scan ids that aren't available)
736 0 : Vector<Int> reqCTscan=cls.scanmap.ctids(reqMSscan);
737 :
738 : // cout << "reqCTscan = " << reqCTscan << endl;
739 :
740 : // For each required CT scan (and thus the MS scan ids requiring it)
741 0 : NewCalTable scanselCT(obsselCT);
742 0 : for (uInt iCTscan=0;iCTscan<reqCTscan.nelements();++iCTscan) {
743 0 : Int& thisCTscan=reqCTscan(iCTscan);
744 :
745 : // Apply thisCTscan selection to the CT
746 : // (only if a meaningful scan is specified)
747 0 : if (thisCTscan!=-1)
748 0 : this->selectOnCT(scanselCT,ct_,"",String::toString(thisCTscan),"","","");
749 :
750 : // The MS scans to be calibrated by thisCTscan (limited by the req superset)
751 : // (could be [-1], which means all)
752 0 : Vector<Int> theseMSscan=cls.scanmap.msids(thisCTscan,reqMSscan);
753 0 : if (theseMSscan.nelements()==1 && theseMSscan[0]<0)
754 0 : theseMSscan.reference(reqMSscan);
755 : // cout << " thisCTscan = " << thisCTscan << ": theseMSscan = " << theseMSscan << endl;
756 :
757 : // Apply theseMSscan selection to the MS
758 : // TBD: reqMSfld = ...
759 :
760 : // The net CT fld required for the MS fld according to the fldmap
761 : // (in principle, this may contain CT fld ids that aren't available)
762 : // NB: currently all [-1] or singles; "some" is TBD
763 0 : Vector<Int> reqCTfld=cls.fldmap.ctids(reqMSfld);
764 :
765 : // cout << " reqCTfld = " << reqCTfld << endl;
766 :
767 :
768 : // For each required CT fld:
769 0 : NewCalTable fldselCT(scanselCT);
770 0 : for (uInt iCTfld=0;iCTfld<reqCTfld.nelements();++iCTfld) {
771 0 : Int& thisCTfld=reqCTfld(iCTfld); // TBD: generalize to multiple fields?
772 :
773 : // Apply thisCTfld selection to the CT
774 0 : if (thisCTfld!=-1)
775 0 : this->selectOnCT(fldselCT,obsselCT,"","",String::toString(thisCTfld),"","");
776 :
777 :
778 : // The MS flds to be calibrated by thisCTfld
779 : // (could be [-1], which means all)
780 0 : Vector<int> theseMSfld=cls.fldmap.msids(thisCTfld,reqMSfld);
781 0 : if (theseMSfld.nelements()==1 && theseMSfld[0]<0)
782 0 : theseMSfld.reference(reqMSfld);
783 :
784 :
785 : // cout << " thisCTfld = " << thisCTfld << ": theseMSfld = " << theseMSfld << endl;
786 :
787 : // Apply theseMSfld selection to the MS
788 : // TBD: reqMSspw = ...
789 :
790 : // ...AND HARDWARE AXES (SPW,ANT) IN INNER LOOPS
791 :
792 : // For each required _MS_ spw:
793 0 : NewCalTable spwselCT(fldselCT);
794 0 : for (uInt iMSspw=0;iMSspw<reqMSspw.nelements();++iMSspw) {
795 0 : Int& thisMSspw=reqMSspw(iMSspw);
796 0 : Int thisCTspw=cls.spwmap(thisMSspw);
797 0 : if (thisCTspw<0) thisCTspw=thisMSspw; // MUST BE DEFINITE!
798 :
799 : // Apply thisCTspw selection to CT
800 0 : this->selectOnCT(spwselCT,fldselCT,"","","",String::toString(thisCTspw),"");
801 :
802 : // Create time-dep interp result container
803 0 : CTCalPatchKey iclTres(icls,thisCTobs,-1,thisCTfld,thisMSspw,-1);
804 0 : clTres_[iclTres]=CLPPResult(nPar_,nFPar_,nChanIn_[thisCTspw],nMSElem_);
805 :
806 0 : NewCalTable antselCT(spwselCT);
807 0 : for (uInt iMSant=0;iMSant<reqMSant.nelements();++iMSant) {
808 0 : Int& thisMSant=reqMSant(iMSant);
809 0 : Int thisCTant=cls.antmap(thisMSant);
810 0 : if (thisCTant<0) thisCTant=thisMSant;
811 :
812 : // Apply thisCTant selection to CT
813 0 : this->selectOnCT(antselCT,spwselCT,"","","","",String::toString(thisCTant));
814 :
815 : // (if null, warn and continue, or throw?)
816 :
817 : // Make the Cal Interpolator (icls is the CL slice index):
818 0 : CTCalPatchKey ici0(icls,thisCTobs,-1,thisCTfld,thisCTspw,thisCTant); // all CT indices
819 0 : CTCalPatchKey ici1(icls,thisCTobs,-1,thisCTfld,thisMSspw,thisMSant); // spw,ant are MS indices
820 : // (NB: Must use thisMSspw,thisMSant above to avoid duplication in resolved spwmap,antmap)
821 :
822 0 : if (ci_.count(ici1)<1) {
823 0 : ciname_[ici1]=ici0.print()+" rows="+String::toString(antselCT.nrow());
824 0 : Array<Float> r(clTres_[iclTres].result(thisMSant));
825 0 : Array<Bool> rf(clTres_[iclTres].resultFlag(thisMSant));
826 0 : ci_[ici1]=(*cttifactoryptr_)(antselCT,cls.tinterp,r,rf);
827 : // cout << "Creating: CT("<<ici1.print() << ") --> CT(" << ici0.print() << ")" << endl;
828 0 : }
829 : else
830 0 : throw(AipsError("Attempted duplicate CTCalPatchKey!"));
831 :
832 : // Now distribute this CTTimeInterp1 instance to all relevant MS indices
833 0 : for (uInt iMSobs=0;iMSobs<theseMSobs.nelements();++iMSobs) {
834 0 : Int& thisMSobs=theseMSobs(iMSobs);
835 0 : for (uInt iMSscan=0;iMSscan<theseMSscan.nelements();++iMSscan) {
836 0 : Int& thisMSscan=theseMSscan(iMSscan);
837 0 : for (uInt iMSfld=0;iMSfld<theseMSfld.nelements();++iMSfld) {
838 0 : Int& thisMSfld=theseMSfld(iMSfld);
839 0 : for (uInt iMSint=0;iMSint<theseMSint.nelements();++iMSint) {
840 0 : Int& thisMSint=theseMSint(iMSint);
841 :
842 0 : MSCalPatchKey ims(thisMSobs,thisMSscan,thisMSfld,thisMSint,thisMSspw,thisMSant);
843 0 : if (msci_.count(ims)<1) {
844 0 : msciname_[ims]=ciname_[ici1];
845 0 : msci_[ims]=ci_[ici1];
846 : }
847 : else
848 0 : throw(AipsError("Attempted duplicate MSCalPatchKey!"));
849 :
850 : // cout << " Patching: MS(" << ims.print() << ") --> CT(" << ici0.print() << ")" << endl;
851 :
852 : // Link these obs,fld,ant,spw to the correct results object
853 : // (as a group over antennas; should move this out of ant loop, really)
854 0 : if (iMSant==0) {
855 0 : MSCalPatchKey imsgroup(thisMSobs,thisMSscan,thisMSfld,thisMSint,thisMSspw,-1);
856 0 : msTres_[imsgroup]=clTres_[iclTres];
857 0 : msFres_[imsgroup]=CLPPResult(); // this will be resized on-demand
858 0 : ctspw_[imsgroup]=thisCTspw;
859 0 : finterp_[imsgroup]=cls.finterp;
860 0 : }
861 0 : } // iMSint
862 : } // iMSfld
863 : } // iMSscan
864 : } // iMSobs
865 :
866 0 : } // iMSant
867 0 : } // iMSspw
868 0 : } // iCTfld
869 0 : } // iCTscan
870 0 : } // iCTobs
871 :
872 0 : } // icls
873 :
874 0 : } // ctor
875 :
876 : // (MS sensitive)
877 0 : CLPatchPanel::CLPatchPanel(const String& ctname,
878 : const MeasurementSet& ms,
879 : const Record& callib,
880 : VisCalEnum::MatrixType mtype,
881 : Int nPar,
882 0 : const CTTIFactoryPtr cttifactoryptr) :
883 0 : ct_(NewCalTable::createCT(ctname,Table::Old,Table::Memory)),
884 0 : ctasms_(), // null, in this context
885 0 : ms_(ms),
886 0 : mtype_(mtype),
887 0 : isCmplx_(false),
888 0 : nPar_(nPar),
889 0 : nFPar_(nPar),
890 :
891 0 : nChanIn_(),
892 0 : freqIn_(),
893 0 : refFreqIn_(),
894 :
895 0 : nMSObs_(ms_.observation().nrow()),
896 0 : nMSFld_(ms_.field().nrow()),
897 0 : nMSSpw_(ms_.spectralWindow().nrow()),
898 0 : nMSAnt_(ms_.antenna().nrow()),
899 0 : nMSElem_(0), // set non-trivially in ctor body
900 0 : nCTObs_(ct_.observation().nrow()),
901 0 : nCTFld_(ct_.field().nrow()),
902 0 : nCTSpw_(ct_.spectralWindow().nrow()),
903 0 : nCTAnt_(ct_.antenna().nrow()),
904 0 : nCTElem_(0), // set non-trivially in ctor body
905 0 : lastresadd_(nMSSpw_,NULL),
906 0 : cttifactoryptr_(cttifactoryptr)
907 :
908 : // elemMap_(),
909 : // conjTab_(),
910 :
911 : {
912 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::CLPatchPanel(<w/MS>)" << endl;
913 :
914 : // ia1dmethod_=ftype(freqType_);
915 : // cout << "ia1dmethod_ = " << ia1dmethod_ << endl;
916 :
917 :
918 0 : switch(mtype_) {
919 0 : case VisCalEnum::GLOBAL: {
920 :
921 0 : throw(AipsError("CLPatchPanel::ctor: No non-Mueller/Jones yet."));
922 :
923 : nCTElem_=1;
924 : nMSElem_=1;
925 : break;
926 : }
927 0 : case VisCalEnum::MUELLER: {
928 0 : nCTElem_=nCTAnt_*(nCTAnt_+1)/2;
929 0 : nMSElem_=nMSAnt_*(nMSAnt_+1)/2;
930 0 : break;
931 : }
932 0 : case VisCalEnum::JONES: {
933 0 : nCTElem_=nCTAnt_;
934 0 : nMSElem_=nMSAnt_;
935 0 : break;
936 : }
937 : }
938 :
939 : // How many _Float_ parameters?
940 0 : isCmplx_=(ct_.keywordSet().asString("ParType")=="Complex");
941 0 : if (isCmplx_) // Complex input
942 0 : nFPar_*=2; // interpolating 2X as many Float values
943 :
944 : // Set CT channel/freq info
945 0 : CTSpWindowColumns ctspw(ct_.spectralWindow());
946 0 : ctspw.numChan().getColumn(nChanIn_);
947 0 : freqIn_.resize(nCTSpw_);
948 0 : for (uInt iCTspw=0;iCTspw<ctspw.nrow();++iCTspw)
949 0 : ctspw.chanFreq().get(iCTspw,freqIn_(iCTspw),true);
950 0 : ctspw.refFrequency().getColumn(refFreqIn_);
951 :
952 : // Initialize maps
953 : // (ctids, msids, result arrays)
954 : // TBD
955 :
956 : //void CLPatchPanel::unpackCLinstance(CalLibSlice cls) {
957 :
958 : // Loop over callib slices
959 : // (a callib slice is one MS selection and interp setup for the present caltable)
960 0 : uInt nkey=callib.nfields();
961 0 : Int icls=-1;
962 0 : for (uInt ikey=0;ikey<nkey;++ikey) {
963 :
964 : // CalTable name might be needed below (eg, for log messages)
965 0 : String ctname=Path(ct_.getPartNames()[0]).baseName().before(".tempMemCalTable");
966 :
967 : // Ignore non-Record members in the callib
968 0 : if (callib.type(ikey)!=TpRecord)
969 0 : continue;
970 :
971 0 : ++icls;
972 :
973 : // The current callib slice
974 0 : CalLibSlice cls(callib.asRecord(ikey),ms_,ct_);
975 :
976 :
977 : // Reference to the cal table that we'll use below
978 0 : NewCalTable ct0(ct_);
979 0 : if (cls.extfldsel!="") {
980 : // Select on the reference table
981 : try {
982 0 : this->selectOnCT(ct0,ct_,"","",cls.extfldsel,"","");
983 0 : } catch ( MSSelectionError err ) {
984 : // Selection failed somehow:
985 0 : throw(AipsError("Problem selecting for multi-field field mapping ('"+cls.extfldsel+"') in caltable="+ctname+": "+err.getMesg()));
986 0 : }
987 :
988 : }
989 :
990 : // Apply callib instance MS selection to the incoming (selected MS)
991 0 : MeasurementSet clsms(ms_);
992 :
993 : // Trap Null selection exceptions, as they are not needed downstream
994 : try {
995 0 : this->selectOnMS(clsms,ms_,cls.obs,cls.scan,cls.fld,cls.ent,cls.spw,"");
996 : }
997 0 : catch ( MSSelectionNullSelection x ) {
998 :
999 : // Warn in logger that this slice doesn't match anything
1000 : // in the selected MS
1001 : //String msname=Path(ms_.tableName()).baseName();
1002 0 : String msname=Path(ms_.getPartNames()[0]).baseName();
1003 0 : logsink_ << LogIO::WARN
1004 : << ". The following callib entry matches no data" << endl
1005 0 : << ". in the selected MS ("+msname+") and will be ignored:" << endl
1006 : << ". " << icls << ":" << endl
1007 0 : << cls.state() << LogIO::POST;
1008 :
1009 : // Just jump to next slice (cal maps for this MS selection unneeded)
1010 0 : continue;
1011 0 : }
1012 :
1013 : // Report this relevant cal lib slice to the logger
1014 0 : logsink_ << LogIO::NORMAL << ". " << icls << ":" << endl
1015 0 : << cls.state() << LogIO::POST;
1016 :
1017 : // What MS indices will be calibrated by this CL instance?
1018 :
1019 : // TBD: we will want to calculate these within the loops below
1020 : // so that per-obs and per-fld subsets will be correctly resolved
1021 : // (and, e.g., nearest field can be obs-dep, etc.)
1022 : // (gmoellen 2018/02/05: still true? Doesn't selection above
1023 : // reduce requirements to the union of obs/fld/intent accounted
1024 : // for by the current cl?)
1025 :
1026 : // OBS Ids in selected MS to be calibrated by this cl instance
1027 0 : Vector<Int> reqMSobs(1,-1); // assume all, indescriminately
1028 0 : if (cls.obs.length()>0)
1029 : // if CL is obs-specific, we must not be indescriminate
1030 0 : reqMSobs.reference(this->getMSuniqueIds(clsms,"obs"));
1031 : //cout << "reqMSobs = " << reqMSobs << endl;
1032 :
1033 : // SCAN Ids in selected MS to be calibrated by this cl instance
1034 0 : Vector<Int> reqMSscan(1,-1); // assume all, indescriminately
1035 0 : if (cls.scan.length()>0)
1036 : // if CL is scan-specific, we must not be indescriminate
1037 0 : reqMSscan.reference(this->getMSuniqueIds(clsms,"scan"));
1038 : //cout << "reqMSscan = " << reqMSscan << endl;
1039 :
1040 : // FIELD Ids in selected MS to be calibrated by this cl instance
1041 0 : Vector<Int> reqMSfld(1,-1); // assume all, indescriminately
1042 0 : if (cls.fld.length()>0) // if selected, maybe we only need a subset
1043 : // if CL is fld-specific, we must not be indescriminate
1044 0 : reqMSfld.reference(this->getMSuniqueIds(clsms,"fld"));
1045 : //cout << "reqMSfld = " << reqMSfld << endl;
1046 :
1047 : // INTENT Ids in selected MS to be calibrated by this cl instance
1048 0 : Vector<Int> reqMSint(1,-1); // assume all
1049 0 : if (cls.ent.length()>0)
1050 0 : reqMSint.reference(this->getMSuniqueIds(clsms,"intent"));
1051 : //cout << "reqMSint = " << reqMSint << endl;
1052 0 : Vector<Int> theseMSint(reqMSint);
1053 :
1054 : // SPW Ids in selected MS to be calibrated by this cl instance
1055 : // We are never indescriminate about spw
1056 0 : Vector<Int> reqMSspw(this->getMSuniqueIds(clsms,"spw"));
1057 : //cout << "reqMSspw = " << reqMSspw << endl;
1058 :
1059 : // ANTENNA in selected MS to be calibrated by this cl instance
1060 : // We are never indescriminate about ant
1061 : // (For now, we will do all MS ants)
1062 0 : Vector<Int> reqMSant(nMSAnt_);
1063 0 : indgen(reqMSant);
1064 : //cout << "reqMSant = " << reqMSant << endl;
1065 :
1066 : // SLICE CalTable by OBS, SCAN, FIELD, SPW, ANT, and map to
1067 : // the corresponding MS indicies
1068 :
1069 : // WE DO TIME-ISH (OBS,SCAN,FLD) AXES IN OUTER LOOPS
1070 :
1071 0 : NewCalTable obsselCT(ct0);
1072 :
1073 : // The net CT obs required for the MS obs according to the obsmap
1074 : // We will create separate interpolator groups for each
1075 0 : Vector<Int> reqCTobs=cls.obsmap.ctids(reqMSobs);
1076 : //cout << "reqCTobs = " << reqCTobs << endl;
1077 : // For each required CT obs (and thus the MS obs ids requiring it)
1078 0 : for (uInt iCTobs=0;iCTobs<reqCTobs.nelements();++iCTobs) {
1079 0 : Int& thisCTobs=reqCTobs(iCTobs);
1080 :
1081 : // The MS OBSs (subset of reqMSobs) to be calibrated by thisCTobs
1082 : // (could be [-1], which means all)
1083 0 : Vector<Int> theseMSobs=cls.obsmap.msids(thisCTobs,reqMSobs);
1084 0 : if (theseMSobs.nelements()==1 && theseMSobs[0]<0)
1085 0 : theseMSobs.reference(reqMSobs);
1086 : //cout << " thisCTobs = " << thisCTobs << ": theseMSobs = " << theseMSobs << endl;
1087 :
1088 : // Apply thisCTobs selection to the CT
1089 : // (only if a meaningful obsid is specified)
1090 : try {
1091 0 : if (thisCTobs!=-1)
1092 0 : this->selectOnCT(obsselCT,ct0,String::toString(thisCTobs),"","","","");
1093 : }
1094 0 : catch (...) { // MSSelectionNullSelection x ) {
1095 :
1096 : // Required CT obs does not exist in the caltable
1097 0 : recordBadMSIndices(theseMSobs,reqMSscan,reqMSfld,reqMSint,Vector<Int>(1,-1)); // all spws
1098 0 : continue; // jump to next CT obs
1099 0 : }
1100 :
1101 :
1102 : // The net CT scan required for the MS scan according to the scanmap
1103 : // We will create separate interpolator groups for each
1104 0 : Vector<Int> reqCTscan=cls.scanmap.ctids(reqMSscan);
1105 : //cout << "reqCTscan = " << reqCTscan << endl;
1106 : // For each required CT scan (and thus the MS scan ids requiring it)
1107 0 : NewCalTable scanselCT(obsselCT);
1108 0 : for (uInt iCTscan=0;iCTscan<reqCTscan.nelements();++iCTscan) {
1109 0 : Int& thisCTscan=reqCTscan(iCTscan);
1110 :
1111 : // The MS SCANs (subset of reqMSscan) to be calibrated by thisCTscan
1112 : // (could be [-1], which means all)
1113 0 : Vector<Int> theseMSscan=cls.scanmap.msids(thisCTscan,reqMSscan);
1114 0 : if (theseMSscan.nelements()==1 && theseMSscan[0]<0)
1115 0 : theseMSscan.reference(reqMSscan);
1116 : //cout << " thisCTscan = " << thisCTscan << ": theseMSscan = " << theseMSscan << endl;
1117 :
1118 : // Apply thisCTscan selection to the CT
1119 : // (only if a meaningful scanid is specified)
1120 : try {
1121 0 : if (thisCTscan!=-1)
1122 0 : this->selectOnCT(scanselCT,obsselCT,"",String::toString(thisCTscan),"","","");
1123 : }
1124 0 : catch (...) { // MSSelectionNullSelection x ) {
1125 :
1126 : // Required CT scan does not exist in the caltable
1127 0 : recordBadMSIndices(reqMSobs,theseMSscan,reqMSfld,reqMSint,Vector<Int>(1,-1)); // all spws
1128 0 : continue; // jump to next CT scan
1129 0 : }
1130 :
1131 :
1132 : // The net CT fld required for the MS fld according to the fldmap
1133 : // NB: currently all [-1] or singles; "some" is TBD
1134 0 : Vector<Int> reqCTfld=cls.fldmap.ctids(reqMSfld);
1135 : //cout << " reqCTfld = " << reqCTfld << endl;
1136 :
1137 : // For each required CT fld:
1138 0 : NewCalTable fldselCT(scanselCT);
1139 0 : for (uInt iCTfld=0;iCTfld<reqCTfld.nelements();++iCTfld) {
1140 0 : Int& thisCTfld=reqCTfld(iCTfld); // TBD: generalize to multiple fields?
1141 :
1142 : // The MS flds to be calibrated by thisCTfld
1143 : // (could be [-1], which means all)
1144 0 : Vector<int> theseMSfld=cls.fldmap.msids(thisCTfld,reqMSfld);
1145 0 : if (theseMSfld.nelements()==1 && theseMSfld[0]<0)
1146 0 : theseMSfld.reference(reqMSfld);
1147 :
1148 : //cout << " thisCTfld = " << thisCTfld << ": theseMSfld = " << theseMSfld << endl;
1149 :
1150 : // Apply thisCTfld selection to the CT
1151 : try {
1152 0 : if (thisCTfld!=-1)
1153 0 : this->selectOnCT(fldselCT,scanselCT,"","",String::toString(thisCTfld),"","");
1154 : }
1155 0 : catch (...) { // MSSelectionNullSelection x ) {
1156 :
1157 : // Required CT fld does not exist in the caltable (for current CT obs, scan)
1158 0 : recordBadMSIndices(theseMSobs,theseMSscan,theseMSfld,reqMSint,Vector<Int>(1,-1)); // all spws
1159 0 : continue; // jump to next fld
1160 0 : }
1161 :
1162 : // ...AND HARDWARE AXES (SPW,ANT) IN INNER LOOPS
1163 :
1164 : // For each required _MS_ spw:
1165 0 : NewCalTable spwselCT(fldselCT);
1166 0 : for (uInt iMSspw=0;iMSspw<reqMSspw.nelements();++iMSspw) {
1167 0 : Int& thisMSspw=reqMSspw(iMSspw);
1168 0 : Int thisCTspw=cls.spwmap(thisMSspw);
1169 0 : if (thisCTspw<0) thisCTspw=thisMSspw; // MUST BE DEFINITE!
1170 :
1171 : //cout << " thisCTspw=" << thisCTspw << "--> thisMSspw="<<thisMSspw<<endl;
1172 :
1173 : // Apply thisCTspw selection to CT
1174 : try {
1175 0 : this->selectOnCT(spwselCT,fldselCT,"","","",String::toString(thisCTspw),"");
1176 : }
1177 0 : catch (...) { // MSSelectionNullSelection x ) {
1178 :
1179 : // Required CT spw does not exist in the caltable (for current CT obs, scan, fld)
1180 0 : recordBadMSIndices(theseMSobs,theseMSscan,theseMSfld,reqMSint,Vector<Int>(1,thisMSspw)); // current spw
1181 0 : continue; // jump to next spw
1182 0 : }
1183 :
1184 : // If this selection fails (zero rows), and exception is thrown.
1185 : // What is the state of antselCT?
1186 : // Is it still the unselected-upon spwselCT?
1187 : // Or is an empty table?
1188 :
1189 :
1190 :
1191 : // Create time-dep interp result container
1192 : // Indexed by CTobs, CTscan, CTfld, MSspw (for all antennas)
1193 0 : CTCalPatchKey iclTres(icls,thisCTobs,thisCTscan,thisCTfld,thisMSspw,-1);
1194 0 : clTres_[iclTres]=CLPPResult(nPar_,nFPar_,nChanIn_[thisCTspw],nMSElem_);
1195 :
1196 0 : NewCalTable antselCT(spwselCT);
1197 0 : Bool doLinkResults(true); // initialize true, so it will happen if relevant code reached
1198 0 : for (uInt iMSant=0;iMSant<reqMSant.nelements();++iMSant) {
1199 0 : Int& thisMSant=reqMSant(iMSant);
1200 0 : Int thisCTant=cls.antmap(thisMSant);
1201 0 : if (thisCTant<0) thisCTant=thisMSant;
1202 :
1203 : // Apply thisCTant selection to CT
1204 : try {
1205 0 : this->selectOnCT(antselCT,spwselCT,"","","","",String::toString(thisCTant));
1206 : }
1207 0 : catch ( MSSelectionNullSelection x ) {
1208 : // Log a warning about the missing antenna
1209 0 : logsink_ << LogIO::WARN << " Found no calibration for MS ant Id=" << thisMSant << " (CT ant Id=" << thisCTant << ")"
1210 : << " in MS spw Id=" << thisMSspw << " (CT spw Id=" << thisCTspw << ") (" << ctname << ")"
1211 0 : << LogIO::POST;
1212 : // Step to next antenna
1213 0 : continue;
1214 0 : }
1215 :
1216 : // Make the Cal Interpolator (icls is the CL slice index):
1217 0 : CTCalPatchKey ici0(icls,thisCTobs,thisCTscan,thisCTfld,thisCTspw,thisCTant); // all CT indices
1218 0 : CTCalPatchKey ici1(icls,thisCTobs,thisCTscan,thisCTfld,thisMSspw,thisMSant); // spw,ant are MS indices
1219 : // (NB: Must use thisMSspw,thisMSant above to avoid duplication in resolved spwmap,antmap)
1220 :
1221 0 : if (ci_.count(ici1)<1) {
1222 0 : ciname_[ici1]=ici0.print()+" rows="+String::toString(antselCT.nrow());
1223 0 : Array<Float> r(clTres_[iclTres].result(thisMSant));
1224 0 : Array<Bool> rf(clTres_[iclTres].resultFlag(thisMSant));
1225 0 : ci_[ici1]=(*cttifactoryptr_)(antselCT,cls.tinterp,r,rf);
1226 : //if (iMSant==0) cout << " Creating: CT("<<ici1.print() << ") --> CT(" << ici0.print() << ") (all antennas)" << endl;
1227 0 : }
1228 : else
1229 0 : throw(AipsError("Attempted duplicate CTCalPatchKey!"));
1230 :
1231 : // Now distribute this CTTimeInterp1 instance to all relevant MS indices
1232 0 : for (uInt iMSobs=0;iMSobs<theseMSobs.nelements();++iMSobs) {
1233 0 : Int& thisMSobs=theseMSobs(iMSobs);
1234 0 : for (uInt iMSscan=0;iMSscan<theseMSscan.nelements();++iMSscan) {
1235 0 : Int& thisMSscan=theseMSscan(iMSscan);
1236 0 : for (uInt iMSfld=0;iMSfld<theseMSfld.nelements();++iMSfld) {
1237 0 : Int& thisMSfld=theseMSfld(iMSfld);
1238 0 : for (uInt iMSint=0;iMSint<theseMSint.nelements();++iMSint) {
1239 0 : Int& thisMSint=theseMSint(iMSint);
1240 :
1241 0 : MSCalPatchKey ims(thisMSobs,thisMSscan,thisMSfld,thisMSint,thisMSspw,thisMSant);
1242 0 : if (msci_.count(ims)<1) {
1243 0 : msciname_[ims]=ciname_[ici1];
1244 0 : msci_[ims]=ci_[ici1];
1245 : }
1246 : else
1247 0 : throw(AipsError("Attempted duplicate MSCalPatchKey!"));
1248 :
1249 : //if (doLinkResults)
1250 : // cout << " Patching: MS(" << ims.print() << ") --> CT(" << ici0.print() << ")" << endl;
1251 :
1252 : // Link these obs,scan,fld,ant,spw to the correct results object
1253 : // (as a group over antennas; should move this out of ant loop, really)
1254 0 : if (doLinkResults) {
1255 0 : MSCalPatchKey imsgroup(thisMSobs,thisMSscan,thisMSfld,thisMSint,thisMSspw,-1);
1256 0 : msTres_[imsgroup]=clTres_[iclTres];
1257 0 : msFres_[imsgroup]=CLPPResult(); // this will be resized on-demand
1258 0 : ctspw_[imsgroup]=thisCTspw;
1259 0 : finterp_[imsgroup]=cls.finterp;
1260 0 : }
1261 0 : } // iMSint
1262 : } // iMSfld
1263 : } // iMSscan
1264 : } // iMSobs
1265 0 : doLinkResults = False; // Don't do it again
1266 0 : } // iMSant
1267 0 : } // iMSspw
1268 0 : } // iCTfld
1269 0 : } // iCTscan
1270 0 : } // iCTobs
1271 :
1272 :
1273 0 : } // icls
1274 :
1275 :
1276 :
1277 0 : } // ctor
1278 :
1279 0 : void CLPatchPanel::recordBadMSIndices(const Vector<Int>& obs, const Vector<Int>& scan,
1280 : const Vector<Int>& fld,
1281 : const Vector<Int>& ent, const Vector<Int>& spw) {
1282 :
1283 :
1284 : // Record bad _MS_ indices
1285 0 : for (uInt iobs=0;iobs<obs.nelements();++iobs) {
1286 0 : for (uInt iscan=0;iscan<scan.nelements();++iscan) {
1287 0 : for (uInt ifld=0;ifld<fld.nelements();++ifld) {
1288 0 : for (uInt ient=0;ient<ent.nelements();++ient) {
1289 0 : for (uInt ispw=0;ispw<spw.nelements();++ispw) {
1290 :
1291 0 : MSCalPatchKey ims(obs[iobs],scan[iscan],fld[ifld],ent[ient],spw[ispw],-1); // All ants
1292 0 : if (badmsciname_.count(ims)<1) {
1293 0 : badmsciname_[ims]=ims.print();
1294 : //cout << " Bad MS indices: " << ims.print() << endl;
1295 : }
1296 0 : }
1297 : }
1298 : }
1299 : }
1300 : }
1301 0 : return;
1302 : }
1303 :
1304 :
1305 0 : void CLPatchPanel::selectOnCT(NewCalTable& ctout,const NewCalTable& ctin,
1306 : const String& obs, const String& scan,
1307 : const String& fld,
1308 : const String& spw, const String& ant1) {
1309 :
1310 0 : String taql("");
1311 0 : if (ant1.length()>0)
1312 0 : taql="ANTENNA1=="+ant1;
1313 :
1314 : // cout << "taql = >>" << taql << "<<" << endl;
1315 :
1316 : // Forward to generic method (sans intent)
1317 0 : CTInterface cti(ctin);
1318 0 : this->selectOnCTorMS(ctout,cti,obs,scan,fld,"",spw,"",taql);
1319 :
1320 0 : }
1321 :
1322 0 : void CLPatchPanel::selectOnMS(MeasurementSet& msout,const MeasurementSet& msin,
1323 : const String& obs, const String& scan,
1324 : const String& fld, const String& ent,
1325 : const String& spw, const String& ant) {
1326 :
1327 : // Forward to generic method
1328 0 : MSInterface msi(msin);
1329 0 : this->selectOnCTorMS(msout,msi,obs,scan,fld,ent,spw,ant,"");
1330 :
1331 0 : }
1332 0 : void CLPatchPanel::selectOnCTorMS(Table& ctout,MSSelectableTable& msst,
1333 : const String& obs, const String& scan,
1334 : const String& fld, const String& ent,
1335 : const String& spw, const String& ant,
1336 : const String& taql) {
1337 :
1338 0 : MSSelection mss;
1339 0 : if (obs.length()>0)
1340 0 : mss.setObservationExpr(obs);
1341 0 : if (scan.length()>0)
1342 0 : mss.setScanExpr(scan);
1343 0 : if (fld.length()>0)
1344 0 : mss.setFieldExpr(fld);
1345 0 : if (ent.length()>0)
1346 0 : mss.setStateExpr(ent);
1347 0 : if (spw.length()>0)
1348 0 : mss.setSpwExpr(spw);
1349 0 : if (ant.length()>0)
1350 0 : mss.setAntennaExpr(ant); // this will not behave as required for Jones caltables (its bl-based selection!)
1351 0 : if (taql.length()>0)
1352 0 : mss.setTaQLExpr(taql);
1353 :
1354 :
1355 0 : TableExprNode ten=mss.toTableExprNode(&msst);
1356 0 : getSelectedTable(ctout,*(msst.table()),ten,"");
1357 :
1358 0 : }
1359 :
1360 0 : Vector<Int> CLPatchPanel::getCLuniqueIds(NewCalTable& ct, String vcol) {
1361 :
1362 0 : ROCTMainColumns ctmc(ct);
1363 :
1364 : // Extract the requested column as a Vector
1365 0 : Vector<Int> colv;
1366 0 : if (vcol=="obs")
1367 0 : ctmc.obsId().getColumn(colv);
1368 0 : else if (vcol=="scan")
1369 0 : ctmc.scanNo().getColumn(colv);
1370 0 : else if (vcol=="fld")
1371 0 : ctmc.fieldId().getColumn(colv);
1372 0 : else if (vcol=="spw")
1373 0 : ctmc.spwId().getColumn(colv);
1374 : else
1375 0 : throw(AipsError("Unsupported column in CLPatchPanel::getCLuniqueIds"));
1376 :
1377 : // Reduce to a unique list
1378 0 : uInt nuniq=genSort(colv,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
1379 0 : colv.resize(nuniq,true);
1380 :
1381 0 : return colv;
1382 :
1383 0 : }
1384 :
1385 0 : Vector<Int> CLPatchPanel::getMSuniqueIds(MeasurementSet& ms, String which) {
1386 :
1387 0 : MSColumns msc(ms);
1388 :
1389 : // Extract the requested column as a Vector
1390 0 : Vector<Int> colv;
1391 0 : if (which=="obs")
1392 0 : msc.observationId().getColumn(colv);
1393 0 : else if (which=="scan")
1394 0 : msc.scanNumber().getColumn(colv);
1395 0 : else if (which=="fld")
1396 0 : msc.fieldId().getColumn(colv);
1397 0 : else if (which=="intent")
1398 0 : msc.stateId().getColumn(colv);
1399 0 : else if (which=="spw")
1400 0 : msc.dataDescId().getColumn(colv); // these are actually ddids!
1401 : else
1402 0 : throw(AipsError("Unsupported column in CLPatchPanel::getCLuniqueIds"));
1403 :
1404 : // Reduce to a unique list
1405 0 : uInt nuniq=genSort(colv,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
1406 0 : colv.resize(nuniq,true);
1407 :
1408 : // In case of spw, translate from ddid
1409 0 : if (which=="spw") {
1410 0 : Vector<Int> spwids;
1411 0 : msc.dataDescription().spectralWindowId().getColumn(spwids);
1412 0 : for (uInt i=0;i<colv.nelements();++i)
1413 0 : colv[i]=spwids[colv[i]];
1414 0 : }
1415 :
1416 : // return the value
1417 0 : return colv;
1418 :
1419 0 : }
1420 :
1421 :
1422 :
1423 : // Destructor
1424 0 : CLPatchPanel::~CLPatchPanel() {
1425 :
1426 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::~CLPatchPanel()" << endl;
1427 :
1428 : // Delete the atomic interpolators
1429 0 : for (std::map<CTCalPatchKey,CTTimeInterp1*>::iterator it=ci_.begin(); it!=ci_.end(); ++it)
1430 0 : delete it->second;
1431 :
1432 0 : }
1433 :
1434 : // Is specific calibration explicitly available for a obs,scan,fld,intent,spw,ant combination?
1435 0 : Bool CLPatchPanel::calAvailable(casacore::Int msobs, casacore::Int msscan, casacore::Int msfld,
1436 : casacore::Int msent, casacore::Int msspw, casacore::Int msant) {
1437 :
1438 0 : const MSCalPatchKey key(msobs,msscan,msfld,msent,msspw,msant);
1439 :
1440 0 : Bool avail=msTres_.count(key)>0;
1441 :
1442 : // if (!avail) {
1443 : // cout << Path(ct_.tableName()).baseName().before(".tempMem") << " stepped over: " << key.print() << endl;
1444 : // }
1445 :
1446 0 : return avail;
1447 :
1448 0 : }
1449 :
1450 : // The specified MS indices are "OK" if not recorded as expected-but-absent as per
1451 : // the callibrary specification. Expected-but-absent occurs when a CL entry
1452 : // indicates (even implicitly) that the MS index combination is supported but
1453 : // the required CalTable indices (via *map params) are not actually available
1454 : // in the CalTable, e.g., a missing spw.
1455 : // Such data cannot be calibrated and this CL cannot be agnostic about it...
1456 : // Note that MSIndicesOK can return true for MS index combinations for which
1457 : // calibration is not actually available, if the CL entry does not purport
1458 : // to support calibrating them. In such cases, calAvailable() returns false, and
1459 : // this CL is agnostic w.r.t. such data and lets it pass thru unchanged.
1460 : // Returns TRUE when the specified MS indices are calibrateable or passable.
1461 0 : Bool CLPatchPanel::MSIndicesOK(casacore::Int msobs, casacore::Int msscan, casacore::Int msfld,
1462 : casacore::Int msent, casacore::Int msspw, casacore::Int msant) {
1463 :
1464 0 : const MSCalPatchKey key(msobs,msscan,msfld,msent,msspw,msant);
1465 0 : Bool bad=badmsciname_.count(key)>0;
1466 : // if (bad) {
1467 : // cout << Path(ct_.tableName()).baseName().before(".tempMem") << " should but can't calibrate: " << key.print() << endl;
1468 : // }
1469 :
1470 : // Return TRUE if NOT bad
1471 0 : return !bad;
1472 :
1473 0 : }
1474 :
1475 0 : Bool CLPatchPanel::interpolate(Cube<Complex>& resultC, Cube<Bool>& resFlag,
1476 : Int msobs, Int msscan, Int msfld, Int msent, Int msspw,
1477 : Double time, Double freq) {
1478 :
1479 0 : Bool newcal(false);
1480 :
1481 : // resultC and resFlag will be unchanged if newcal remains false
1482 0 : Cube<Float> f; // temporary to reference Float interpolation result
1483 0 : newcal=interpolate(f,resFlag,msobs,msscan,msfld,msent,msspw,time,freq);
1484 0 : if (newcal)
1485 : // convert to complex and have resultC take over ownership
1486 0 : resultC.reference(RIorAPArray(f).c());
1487 :
1488 0 : return newcal; // If T, calling scope should act nontrivially, if necessary
1489 :
1490 0 : }
1491 :
1492 0 : Bool CLPatchPanel::interpolate(Cube<Complex>& resultC, Cube<Bool>& resFlag,
1493 : Int msobs, Int msscan, Int msfld, Int msent, Int msspw,
1494 : Double time, const Vector<Double>& freq) {
1495 :
1496 0 : Bool newcal(false);
1497 :
1498 : // resultC and resFlag will be unchanged if newcal remains false
1499 0 : Cube<Float> f; // temporary to reference Float interpolation result
1500 0 : newcal=interpolate(f,resFlag,msobs,msscan,msfld,msent,msspw,time,freq);
1501 :
1502 0 : if (newcal)
1503 : // convert to complex and have resultC take over ownership
1504 0 : resultC.reference(RIorAPArray(f).c());
1505 :
1506 0 : return newcal; // If T, calling scope should act nontrivially, if necessary
1507 :
1508 0 : }
1509 :
1510 :
1511 0 : Bool CLPatchPanel::interpolate(Cube<Float>& resultR, Cube<Bool>& resFlag,
1512 : Int msobs, Int msscan, Int msfld, Int msent, Int msspw,
1513 : Double time, Double freq) {
1514 :
1515 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::interpolate(...)" << endl;
1516 :
1517 : // NB: resultR and resFlag will be unchanged if newcal remains false
1518 0 : Bool newcal(false);
1519 :
1520 : // Suppled arrays reference the result (if available)
1521 0 : MSCalPatchKey ires(msobs,msscan,msfld,msent,msspw,-1);
1522 :
1523 : // Trap lack of available calibration for requested obs,fld,intent,spw
1524 0 : if (msTres_.count(ires)==0) {
1525 0 : throw(AipsError("No calibration arranged for "+ires.print()+
1526 0 : " in callib for caltable="+
1527 0 : Path(ct_.tableName()).baseName().before(".tempMemCalTable") ));
1528 : }
1529 :
1530 : // If result_ is at a new address (cf last time, this spw), treat as new
1531 : // TBD: do this is a less obscure way... (e.g., invent the CLCalGroup(nant))
1532 0 : if (lastresadd_(msspw)!=msTres_[ires].result_.data())
1533 0 : newcal=true;
1534 0 : lastresadd_(msspw)=msTres_[ires].result_.data();
1535 :
1536 : // Loop over _output_ elements
1537 0 : for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1538 : // Call fully _patched_ time-interpolator, keeping track of 'newness'
1539 : // fills ctTres_ implicitly
1540 0 : MSCalPatchKey ims(msobs,msscan,msfld,msent,msspw,iMSElem);
1541 0 : if (msci_.count(ims)>0) {
1542 0 : if (freq>0.0)
1543 0 : newcal|=msci_[ims]->interpolate(time,freq);
1544 : else
1545 0 : newcal|=msci_[ims]->interpolate(time);
1546 : }
1547 0 : }
1548 :
1549 0 : if (newcal) {
1550 0 : resultR.reference(msTres_[ires].result_);
1551 0 : resFlag.reference(msTres_[ires].resultFlag_);
1552 : }
1553 :
1554 0 : return newcal; // If true, calling scope should act
1555 0 : }
1556 :
1557 :
1558 0 : Bool CLPatchPanel::interpolate(Cube<Float>& resultR, Cube<Bool>& resFlag,
1559 : Int msobs, Int msscan, Int msfld, Int msent, Int msspw,
1560 : Double time, const Vector<Double>& freq) {
1561 :
1562 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::interpolate(Cube<F>,...,Vector<D>)" << endl;
1563 :
1564 : // NB: resultR and resFlag will be unchanged if newcal remains false
1565 0 : Bool newcal(false);
1566 :
1567 : // Suppled arrays reference the result (if available)
1568 0 : MSCalPatchKey ires(msobs,msscan,msfld,msent,msspw,-1);
1569 :
1570 : // Trap lack of available calibration for requested obs,scan,fld,intent,spw
1571 0 : if (msTres_.count(ires)==0) {
1572 0 : throw(AipsError("No calibration arranged for "+ires.print()+
1573 0 : " in callib for caltable="+
1574 0 : Path(ct_.tableName()).baseName().before(".tempMemCalTable") ));
1575 : }
1576 :
1577 : // If result_ is at a new address (cf last time, this msspw), treat as new
1578 : // TBD: do this is a less obscure way...
1579 0 : if (lastresadd_(msspw)!=msTres_[ires].result_.data())
1580 0 : newcal=true;
1581 0 : lastresadd_(msspw)=msTres_[ires].result_.data();
1582 :
1583 : // Sometimes we need to force the freq interp, even if the time-interp isn't new
1584 0 : Bool forceFinterp=false || newcal;
1585 :
1586 : // The follow occurs unnecessarily in mosaics when newcal=false (i.e., when resampleInFreq won't be called below)
1587 0 : uInt nMSChan=freq.nelements(); // The number of requested channels
1588 0 : if (msFres_.count(ires)>0) {
1589 0 : if (msFres_[ires].result_.nelements()==0 ||
1590 0 : msFres_[ires].result(0).ncolumn()!=nMSChan) {
1591 0 : msFres_[ires].resize(nPar_,nFPar_,nMSChan,nMSElem_);
1592 : }
1593 : }
1594 : else
1595 0 : throw(AipsError("Error finding msFres_ in CLPatchPanel::interpolate(C<F>,...,V<D>)"));
1596 :
1597 : // Loop over _output_ antennas
1598 0 : for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1599 : // Call time interpolation calculation; resample in freq if new
1600 : // (fills msTRes_ implicitly)
1601 :
1602 0 : MSCalPatchKey ims(msobs,msscan,msfld,msent,msspw,iMSElem);
1603 0 : if (msci_.count(ims)>0) {
1604 0 : if (msci_[ims]->interpolate(time) || forceFinterp) {
1605 :
1606 : // Resample in frequency
1607 0 : Matrix<Float> fR( msFres_[ires].result(iMSElem) );
1608 0 : Matrix<Bool> fRflg( msFres_[ires].resultFlag(iMSElem) );
1609 0 : Matrix<Float> tR( msTres_[ires].result(iMSElem) );
1610 0 : Matrix<Bool> tRflg( msTres_[ires].resultFlag(iMSElem) );
1611 0 : resampleInFreq(fR,fRflg,freq,tR,tRflg,freqIn_(ctspw_[ires]),finterp_[ires]);
1612 :
1613 : // Calibration is new
1614 0 : newcal=true;
1615 0 : }
1616 : }
1617 0 : }
1618 :
1619 0 : if (newcal) {
1620 : // Supplied arrays to reference the result
1621 0 : resultR.reference(msFres_[ires].result_);
1622 0 : resFlag.reference(msFres_[ires].resultFlag_);
1623 : }
1624 :
1625 0 : return newcal;
1626 0 : }
1627 :
1628 0 : Bool CLPatchPanel::getTresult(Cube<Float>& resultR, Cube<Bool>& resFlag,
1629 : Int obs, Int scan, Int fld, Int ent, Int spw) {
1630 :
1631 0 : MSCalPatchKey mskey(obs,scan,fld,ent,spw,-1);
1632 :
1633 0 : if (msTres_.count(mskey)==0)
1634 0 : throw(AipsError("No calibration arranged for "+mskey.print()+
1635 0 : " in callib for caltable="+
1636 0 : Path(ct_.tableName()).baseName().before(".tempMemCalTable") ));
1637 :
1638 : // Reference the requested Cube
1639 0 : resultR.reference(msTres_[mskey].result_);
1640 0 : resFlag.reference(msTres_[mskey].resultFlag_);
1641 :
1642 0 : return true;
1643 :
1644 0 : }
1645 :
1646 :
1647 :
1648 :
1649 0 : void CLPatchPanel::listmappings() {
1650 :
1651 0 : cout << "CalTable: " << ct_.tableName() << endl;
1652 0 : cout << "There are " << ci_.size() << " cal interpolation engines." << endl;
1653 0 : cout << "There are " << msci_.size() << " unique MS id combinations mapped to them:" << endl;
1654 0 : for (std::map<MSCalPatchKey,String>::iterator it=msciname_.begin(); it!=msciname_.end(); ++it)
1655 0 : cout << "MS (" << it->first.print() << ") --> CL (" << it->second << ")" << endl;
1656 :
1657 0 : cout << endl << "There are " << badmsciname_.size() << " expected but ABSENT MS id combinations (all ants):" << endl;
1658 0 : for (std::map<MSCalPatchKey,String>::iterator it=badmsciname_.begin(); it!=badmsciname_.end(); ++it)
1659 0 : cout << "MS (" << it->first.print() << ")" << endl;
1660 :
1661 0 : cout << endl;
1662 :
1663 0 : }
1664 :
1665 :
1666 :
1667 :
1668 : // Report state
1669 0 : void CLPatchPanel::state() {
1670 :
1671 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::state()" << endl;
1672 :
1673 0 : cout << "-state--------" << endl;
1674 0 : cout << " ct_ = " << ct_.tableName() << endl;
1675 0 : cout << boolalpha;
1676 0 : cout << " isCmplx_ = " << isCmplx_ << endl;
1677 0 : cout << " nPar_ = " << nPar_ << endl;
1678 0 : cout << " nFPar_ = " << nFPar_ << endl;
1679 0 : cout << " nMSObs_ = " << nMSObs_ << endl;
1680 0 : cout << " nMSFld_ = " << nMSFld_ << endl;
1681 0 : cout << " nMSSpw_ = " << nMSSpw_ << endl;
1682 0 : cout << " nMSAnt_ = " << nMSAnt_ << endl;
1683 0 : cout << " nMSElem_ = " << nMSElem_ << endl;
1684 0 : cout << " nCTObs_ = " << nCTObs_ << endl;
1685 0 : cout << " nCTFld_ = " << nCTFld_ << endl;
1686 0 : cout << " nCTSpw_ = " << nCTSpw_ << endl;
1687 0 : cout << " nCTAnt_ = " << nCTAnt_ << endl;
1688 0 : cout << " nCTElem_ = " << nCTElem_ << endl;
1689 :
1690 : // cout << " timeType_ = " << timeType_ << endl;
1691 : // cout << " freqType_ = " << freqType_ << endl;
1692 0 : }
1693 :
1694 : // Resample in frequency
1695 0 : void CLPatchPanel::resampleInFreq(Matrix<Float>& fres,Matrix<Bool>& fflg,const Vector<Double>& fout,
1696 : Matrix<Float>& tres,Matrix<Bool>& tflg,const Vector<Double>& fin,
1697 : String finterp) {
1698 :
1699 : if (CTPATCHPANELVERB) cout << " CLPatchPanel::resampleInFreq(...)" << endl;
1700 :
1701 : // if no good solutions coming in, return flagged
1702 0 : if (nfalse(tflg)==0) {
1703 0 : fflg.set(true);
1704 0 : return;
1705 : }
1706 :
1707 0 : Int flparmod=nFPar_/nPar_; // for indexing the flag Matrices on the par axis
1708 :
1709 0 : Bool unWrapPhase=flparmod>1;
1710 :
1711 : // cout << "nFPar_,nPar_,flparmod = " << nFPar_ << "," << nPar_ << "," << flparmod << endl;
1712 :
1713 0 : fres=0.0;
1714 :
1715 0 : for (Int ifpar=0;ifpar<nFPar_;++ifpar) {
1716 :
1717 : // Slice by par (each Matrix row)
1718 0 : Vector<Float> fresi(fres.row(ifpar)), tresi(tres.row(ifpar));
1719 0 : Vector<Bool> fflgi(fflg.row(ifpar/flparmod)), tflgi(tflg.row(ifpar/flparmod));
1720 :
1721 : // Mask time result by flags
1722 0 : Vector<Double> mfin=fin(!tflgi).getCompressedArray();
1723 :
1724 0 : if (mfin.nelements()==0) {
1725 : // cout << ifpar << " All chans flagged!" << endl;
1726 : // Everything flagged this par
1727 : // Just flag, zero and go on to the next one
1728 0 : fflgi.set(true);
1729 0 : fresi.set(0.0);
1730 0 : continue;
1731 : }
1732 :
1733 0 : mfin/=1.0e9; // in GHz
1734 0 : Vector<Float> mtresi=tresi(!tflgi).getCompressedArray();
1735 :
1736 : // Trap case of same in/out frequencies
1737 0 : if (fout.nelements()==mfin.nelements() && allNear(fout,mfin,1.e-10)) {
1738 : // Just copy
1739 0 : fresi=mtresi;
1740 0 : fflgi.set(false); // none are flagged
1741 0 : continue;
1742 : }
1743 :
1744 0 : if (ifpar%2==1 && unWrapPhase) {
1745 0 : for (uInt i=1;i<mtresi.nelements();++i) {
1746 0 : while ( (mtresi(i)-mtresi(i-1))>C::pi ) mtresi(i)-=C::_2pi;
1747 0 : while ( (mtresi(i)-mtresi(i-1))<-C::pi ) mtresi(i)+=C::_2pi;
1748 : }
1749 : }
1750 :
1751 :
1752 : // TBD: ensure phases tracked on freq axis...
1753 :
1754 : // TBD: handle flags carefully!
1755 : // (i.e., flag gaps larger than user's "freach")
1756 : // For now,just unset them
1757 0 : fflgi.set(false);
1758 :
1759 :
1760 : // Set flags carefully
1761 0 : resampleFlagsInFreq(fflgi,fout,tflgi,fin,finterp);
1762 :
1763 :
1764 : // Always use nearest on edges
1765 : // TBD: trap cases where frequencies don't overlap at all
1766 : // (fout(hi)<mfin(0) || fout(lo)> mfin(ihi))..... already OK (lo>hi)?
1767 : // TBD: optimize the following by forming Slices in the
1768 : // while loops and doing Array assignment once afterwords
1769 :
1770 0 : Int nfreq=fout.nelements();
1771 0 : Int lo=0;
1772 0 : Int hi=fresi.nelements()-1;
1773 0 : Double inlo(mfin(0));
1774 0 : Int ihi=mtresi.nelements()-1;
1775 0 : Double inhi(mfin(ihi));
1776 :
1777 : // Handle 'nearest' extrapolation in sideband-dep way
1778 0 : Bool inUSB(inhi>inlo);
1779 0 : Bool outUSB(fout(hi)>fout(lo));
1780 0 : if (inUSB) {
1781 0 : if (outUSB) {
1782 0 : while (lo<nfreq && fout(lo)<=inlo) fresi(lo++)=mtresi(0);
1783 0 : while (hi>-1 && fout(hi)>=inhi) fresi(hi--)=mtresi(ihi);
1784 : }
1785 : else { // "outLSB"
1786 0 : while (lo<nfreq && fout(lo)>=inhi) fresi(lo++)=mtresi(ihi);
1787 0 : while (hi>-1 && fout(hi)<=inlo) fresi(hi--)=mtresi(0);
1788 : }
1789 : }
1790 : else { // "inLSB"
1791 0 : if (outUSB) {
1792 0 : while (lo<nfreq && fout(lo)<=inhi) fresi(lo++)=mtresi(ihi);
1793 0 : while (hi>-1 && fout(hi)>=inlo) fresi(hi--)=mtresi(0);
1794 : }
1795 : else { // "outLSB"
1796 0 : while (lo<nfreq && fout(lo)>=inlo) fresi(lo++)=mtresi(0);
1797 0 : while (hi>-1 && fout(hi)<=inhi) fresi(hi--)=mtresi(ihi);
1798 : }
1799 : }
1800 :
1801 : // cout << "lo, hi = " << lo << ","<<hi << endl;
1802 :
1803 0 : if (lo>hi) continue; // Frequencies didn't overlap, nearest was used
1804 :
1805 : // Use InterpolateArray1D to fill in the middle
1806 0 : IPosition blc(1,lo), trc(1,hi);
1807 0 : Vector<Float> slfresi(fresi(blc,trc));
1808 0 : Vector<Double> slfout(fout(blc,trc));
1809 :
1810 0 : InterpolateArray1D<Double,Float>::interpolate(slfresi,slfout,mfin,mtresi,this->ftype(finterp));
1811 :
1812 0 : }
1813 : }
1814 :
1815 :
1816 :
1817 0 : void CLPatchPanel::resampleFlagsInFreq(Vector<Bool>& flgout,const Vector<Double>& fout,
1818 : Vector<Bool>& flgin,const Vector<Double>& fin,
1819 : String finterp) {
1820 :
1821 : // cout << "resampleFlagsInFreq" << endl;
1822 :
1823 : #define NEAREST InterpolateArray1D<Double,Float>::nearestNeighbour
1824 : #define LINEAR InterpolateArray1D<Double,Float>::linear
1825 : #define CUBIC InterpolateArray1D<Double,Float>::cubic
1826 : #define SPLINE InterpolateArray1D<Double,Float>::spline
1827 :
1828 :
1829 : // Handle chan-dep flags
1830 0 : if (finterp.contains("flag")) {
1831 :
1832 0 : InterpolateArray1D<Double,Float>::InterpolationMethod interpMeth=this->ftype(finterp);
1833 :
1834 0 : Vector<Double> finGHz=fin/1e9;
1835 :
1836 : // Determine implied mode-dep flags indexed by channel registration
1837 0 : uInt nflg=flgin.nelements();
1838 0 : Vector<Bool> flreg(nflg,false);
1839 0 : switch (interpMeth) {
1840 0 : case NEAREST: {
1841 : // Just use input flags
1842 0 : flreg.reference(flgin);
1843 0 : break;
1844 : }
1845 0 : case LINEAR: {
1846 0 : for (uInt i=0;i<nflg-1;++i)
1847 0 : flreg[i]=(flgin[i] || flgin[i+1]);
1848 0 : flreg[nflg-1]=flreg[nflg-2];
1849 0 : break;
1850 : }
1851 0 : case CUBIC:
1852 : case SPLINE: {
1853 0 : for (uInt i=1;i<nflg-2;++i)
1854 0 : flreg[i]=(flgin[i-1] || flgin[i] || flgin[i+1] || flgin[i+2]);
1855 0 : flreg[0]=flreg[1];
1856 0 : flreg[nflg-2]=flreg[nflg-3];
1857 0 : flreg[nflg-1]=flreg[nflg-3];
1858 0 : break;
1859 : }
1860 : }
1861 :
1862 :
1863 :
1864 : // Now step through requested chans, setting flags
1865 0 : uInt ireg=0;
1866 0 : uInt nflgout=flgout.nelements();
1867 0 : for (uInt iflgout=0;iflgout<nflgout;++iflgout) {
1868 :
1869 : // Find nominal registration (the _index_ just left)
1870 0 : Bool exact(false);
1871 0 : ireg=binarySearch(exact,finGHz,fout(iflgout),nflg,0);
1872 :
1873 : // If registration is exact, assign verbatim
1874 : // NB: the calibration value calculation occurs agnostically w.r.t. flags,
1875 : // so the calculated value should also match
1876 : // TBD: Add "|| near(finGHz[ireg],fout(iflgout),1e-10) in case "exact" has
1877 : // precision issues?
1878 0 : if (exact) {
1879 0 : flgout[iflgout]=flgin[ireg];
1880 0 : continue;
1881 : }
1882 :
1883 : // Not exact, so carefully handle bracketing
1884 0 : if (ireg>0)
1885 0 : ireg-=1;
1886 0 : ireg=min(ireg,nflg-1);
1887 :
1888 : //while (finGHz(ireg)<=fout(iflgout) && ireg<nflg-1) {
1889 : // ireg+=1; // USB specific!
1890 : //}
1891 : //if (ireg>0 && finGHz(ireg)!=fout(iflgout)) --ireg; // registration is one sample prior
1892 :
1893 : // refine registration by interp type
1894 0 : switch (interpMeth) {
1895 0 : case NEAREST: {
1896 : // nearest might be forward sample
1897 0 : if ( ireg<(nflg-1) &&
1898 0 : abs(fout[iflgout]-finGHz[ireg])>abs(finGHz[ireg+1]-fout[iflgout]) )
1899 0 : ireg+=1;
1900 0 : break;
1901 : }
1902 0 : case LINEAR: {
1903 0 : if (ireg==(nflg-1)) // need one more sample to the right
1904 0 : ireg-=1;
1905 0 : break;
1906 : }
1907 0 : case CUBIC:
1908 : case SPLINE: {
1909 0 : if (ireg==0) ireg+=1; // need one more sample to the left
1910 0 : if (ireg>(nflg-3)) ireg=nflg-3; // need two more samples to the right
1911 0 : break;
1912 : }
1913 : }
1914 :
1915 : // Assign effective flag
1916 0 : flgout[iflgout]=flreg[ireg];
1917 :
1918 : /*
1919 : cout << iflgout << " "
1920 : << ireg << " "
1921 : << flreg[ireg]
1922 : << endl;
1923 : */
1924 :
1925 : }
1926 :
1927 0 : }
1928 : else
1929 : // We are interp/extrap-olating gaps absolutely
1930 0 : flgout.set(false);
1931 :
1932 0 : }
1933 :
1934 0 : InterpolateArray1D<Double,Float>::InterpolationMethod CLPatchPanel::ftype(String& strtype) {
1935 0 : if (strtype.contains("nearest"))
1936 0 : return InterpolateArray1D<Double,Float>::nearestNeighbour;
1937 0 : if (strtype.contains("linear"))
1938 0 : return InterpolateArray1D<Double,Float>::linear;
1939 0 : if (strtype.contains("cubic"))
1940 0 : return InterpolateArray1D<Double,Float>::cubic;
1941 0 : if (strtype.contains("spline"))
1942 0 : return InterpolateArray1D<Double,Float>::spline;
1943 :
1944 : // cout << "Using linear for freq interpolation as last resort." << endl;
1945 0 : return InterpolateArray1D<Double,Float>::linear;
1946 :
1947 :
1948 : }
1949 :
1950 :
1951 : } //# NAMESPACE CASA - END
|
