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 33880 : CalPatchKey::CalPatchKey(IPosition keyids) :
54 33880 : cpk_(keyids.asVector())
55 33880 : {}
56 :
57 : // Lexographical lessthan
58 396418 : 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 2397950 : for (Int i=0;i<6;++i) {
66 2267598 : if (cpk_[i]>-1 && other.cpk_[i]>-1 && cpk_[i]!=other.cpk_[i])
67 : // both non-negative and not equal, so evaluate element <
68 266066 : return cpk_[i]<other.cpk_[i];
69 : }
70 : // All apparently equal so "<" is false
71 130352 : return false;
72 : }
73 :
74 33696 : MSCalPatchKey::MSCalPatchKey(Int obs,Int scan,Int fld,Int ent,Int spw,Int ant) :
75 33696 : CalPatchKey(IPosition(6,obs,scan,fld,ent,spw,ant)),
76 33696 : obs_(obs),scan_(scan),fld_(fld),ent_(ent),spw_(spw),ant_(ant)
77 33696 : {}
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 184 : CTCalPatchKey::CTCalPatchKey(Int clsl,Int obs,Int scan,Int fld,Int spw,Int ant) :
91 184 : CalPatchKey(IPosition(6,clsl,obs,scan,fld,spw,ant)),
92 184 : clsl_(clsl),obs_(obs),scan_(scan),fld_(fld),spw_(spw),ant_(ant)
93 184 : {}
94 :
95 : // text output
96 88 : String CTCalPatchKey::print() const {
97 176 : return "cl="+(clsl_<0 ? "*" : String::toString(clsl_))+" "
98 352 : "obs="+(obs_<0 ? "*" : String::toString(obs_))+" "
99 352 : "scan="+(scan_<0 ? "*" : String::toString(scan_))+" "
100 352 : "fld="+(fld_<0 ? "*" : String::toString(fld_))+" "
101 352 : "spw="+(spw_<0 ? "*" : String::toString(spw_))+" "
102 264 : "ant="+(ant_<0 ? "*" : String::toString(ant_));
103 : }
104 :
105 :
106 10 : CalMap::CalMap() :
107 10 : vcalmap_()
108 10 : {}
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 96 : Int CalMap::operator()(Int msid) const {
118 : // TBD: reconsider algorithm (maybe just return msid if over-run?)
119 96 : Int ncalmap=vcalmap_.nelements();
120 192 : return (msid<ncalmap ? vcalmap_(msid) :
121 192 : (ncalmap>0 ? vcalmap_(ncalmap-1) : msid) ); // Avoid going off end
122 : }
123 :
124 6 : Vector<Int> CalMap::ctids(const Vector<Int>& msids) const {
125 6 : 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 6 : if (ncalmap<1 ||
131 0 : (ncalmap==1 && vcalmap_[0]<0))
132 6 : 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 6 : Vector<Int> CalMap::msids(Int ctid,const Vector<Int>& superset) const {
154 :
155 6 : 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 6 : if (ncalmap<1 ||
160 0 : (ncalmap==1 && vcalmap_[0]<0) )
161 6 : 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/M_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 2 : CalLibSlice::CalLibSlice(const Record& clslice,
403 : const MeasurementSet& ms,
404 2 : const NewCalTable& ct) :
405 2 : obs(),scan(),fld(),ent(),spw(),
406 2 : tinterp(),finterp(),
407 2 : obsmap(),scanmap(),fldmap(),spwmap(),antmap(),
408 2 : extfldsel("")
409 : {
410 :
411 2 : if (clslice.isDefined("obs")) {
412 2 : obs=clslice.asString("obs");
413 : }
414 2 : if (clslice.isDefined("scan")) {
415 0 : scan=clslice.asString("scan");
416 : }
417 2 : if (clslice.isDefined("field")) {
418 2 : fld=clslice.asString("field");
419 : }
420 2 : if (clslice.isDefined("intent")) {
421 2 : ent=clslice.asString("intent");
422 : }
423 2 : if (clslice.isDefined("spw")) {
424 2 : spw=clslice.asString("spw");
425 : }
426 :
427 2 : if (clslice.isDefined("tinterp")) {
428 2 : tinterp=clslice.asString("tinterp");
429 2 : if (tinterp=="") {
430 0 : tinterp="linear";
431 : }
432 : }
433 2 : if (clslice.isDefined("finterp")) {
434 2 : finterp=clslice.asString("finterp");
435 : }
436 :
437 2 : if (clslice.isDefined("obsmap")) {
438 : //cout << "obsmap.dataType() = " << clslice.dataType("obsmap") << endl;
439 2 : if (clslice.dataType("obsmap")==TpArrayInt)
440 0 : obsmap=CalMap(Vector<Int>(clslice.asArrayInt("obsmap")));
441 2 : if (clslice.dataType("obsmap")==TpString)
442 0 : obsmap=ObsCalMap(clslice.asString("obsmap"),ms);
443 : }
444 2 : 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 2 : if (clslice.isDefined("fldmap")) {
452 2 : if (clslice.dataType("fldmap")==TpArrayInt)
453 0 : fldmap=FieldCalMap(clslice.asArrayInt("fldmap"));
454 2 : if (clslice.dataType("fldmap")==TpString)
455 0 : fldmap=FieldCalMap(clslice.asString("fldmap"),ms,ct,extfldsel);
456 : }
457 2 : if (clslice.isDefined("spwmap")) {
458 : // cout << "spwmap.dataType() = " << clslice.dataType("spwmap") << endl;
459 2 : if (clslice.dataType("spwmap")==TpArrayInt)
460 0 : spwmap=CalMap(clslice.asArrayInt("spwmap"));
461 : }
462 2 : if (clslice.isDefined("antmap")) {
463 : // cout << "antmap.dataType() = " << clslice.dataType("antmap") << endl;
464 2 : if (clslice.dataType("antmap")==TpArrayInt)
465 0 : antmap=CalMap(clslice.asArrayInt("antmap"));
466 : }
467 :
468 :
469 2 : }
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 2 : String CalLibSlice::state() {
533 :
534 2 : ostringstream o;
535 :
536 4 : o << " MS: obs="+obs+" scan="+scan+" fld="+fld+" intent="+ent+" spw="+spw << endl
537 4 : << " CT: tinterp="+tinterp << " finterp="+finterp << endl
538 4 : << " obsmap=" << obsmap.vmap()
539 4 : << " scanmap=" << scanmap.vmap()
540 2 : << " fldmap=";
541 :
542 2 : if (extfldsel=="")
543 2 : o << fldmap.vmap();
544 : else
545 0 : o << extfldsel+" (by selection)";
546 :
547 2 : o << endl
548 2 : << " spwmap=" << spwmap.vmap()
549 4 : << " antmap=" << antmap.vmap()
550 2 : << endl;
551 :
552 4 : return String(o);
553 2 : }
554 :
555 :
556 32 : CLPPResult::CLPPResult() :
557 32 : result_(),resultFlag_()
558 32 : {}
559 :
560 0 : CLPPResult::CLPPResult(const IPosition& shape) :
561 0 : result_(shape,0.0),
562 0 : resultFlag_(shape,true)
563 0 : {}
564 8 : CLPPResult::CLPPResult(uInt nPar,uInt nFPar,uInt nChan,uInt nElem) :
565 8 : result_(nFPar,nChan,nElem,0.0),
566 8 : resultFlag_(nPar,nChan,nElem,true)
567 8 : {}
568 :
569 24 : CLPPResult& CLPPResult::operator=(const CLPPResult& other) {
570 : // Avoid Array deep copies!
571 24 : result_.reference(other.result_);
572 24 : resultFlag_.reference(other.resultFlag_);
573 24 : 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 2 : CLPatchPanel::CLPatchPanel(const String& ctname,
878 : const MeasurementSet& ms,
879 : const Record& callib,
880 : VisCalEnum::MatrixType mtype,
881 : Int nPar,
882 2 : const CTTIFactoryPtr cttifactoryptr) :
883 2 : ct_(NewCalTable::createCT(ctname,Table::Old,Table::Memory)),
884 2 : ctasms_(), // null, in this context
885 2 : ms_(ms),
886 2 : mtype_(mtype),
887 2 : isCmplx_(false),
888 2 : nPar_(nPar),
889 2 : nFPar_(nPar),
890 :
891 2 : nChanIn_(),
892 2 : freqIn_(),
893 2 : refFreqIn_(),
894 :
895 2 : nMSObs_(ms_.observation().nrow()),
896 2 : nMSFld_(ms_.field().nrow()),
897 2 : nMSSpw_(ms_.spectralWindow().nrow()),
898 2 : nMSAnt_(ms_.antenna().nrow()),
899 2 : nMSElem_(0), // set non-trivially in ctor body
900 2 : nCTObs_(ct_.observation().nrow()),
901 2 : nCTFld_(ct_.field().nrow()),
902 2 : nCTSpw_(ct_.spectralWindow().nrow()),
903 2 : nCTAnt_(ct_.antenna().nrow()),
904 2 : nCTElem_(0), // set non-trivially in ctor body
905 2 : lastresadd_(nMSSpw_,NULL),
906 6 : 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 2 : 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 2 : case VisCalEnum::JONES: {
933 2 : nCTElem_=nCTAnt_;
934 2 : nMSElem_=nMSAnt_;
935 2 : break;
936 : }
937 : }
938 :
939 : // How many _Float_ parameters?
940 2 : isCmplx_=(ct_.keywordSet().asString("ParType")=="Complex");
941 2 : if (isCmplx_) // Complex input
942 0 : nFPar_*=2; // interpolating 2X as many Float values
943 :
944 : // Set CT channel/freq info
945 4 : CTSpWindowColumns ctspw(ct_.spectralWindow());
946 2 : ctspw.numChan().getColumn(nChanIn_);
947 2 : freqIn_.resize(nCTSpw_);
948 10 : for (uInt iCTspw=0;iCTspw<ctspw.nrow();++iCTspw)
949 8 : ctspw.chanFreq().get(iCTspw,freqIn_(iCTspw),true);
950 2 : 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 2 : uInt nkey=callib.nfields();
961 2 : Int icls=-1;
962 8 : for (uInt ikey=0;ikey<nkey;++ikey) {
963 :
964 : // CalTable name might be needed below (eg, for log messages)
965 6 : String ctname=Path(ct_.getPartNames()[0]).baseName().before(".tempMemCalTable");
966 :
967 : // Ignore non-Record members in the callib
968 6 : if (callib.type(ikey)!=TpRecord)
969 4 : continue;
970 :
971 2 : ++icls;
972 :
973 : // The current callib slice
974 4 : CalLibSlice cls(callib.asRecord(ikey),ms_,ct_);
975 :
976 :
977 : // Reference to the cal table that we'll use below
978 2 : NewCalTable ct0(ct_);
979 2 : 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 2 : MeasurementSet clsms(ms_);
992 :
993 : // Trap Null selection exceptions, as they are not needed downstream
994 : try {
995 2 : 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 2 : logsink_ << LogIO::NORMAL << ". " << icls << ":" << endl
1015 2 : << 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 2 : Vector<Int> reqMSobs(1,-1); // assume all, indescriminately
1028 2 : 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 2 : Vector<Int> reqMSscan(1,-1); // assume all, indescriminately
1035 2 : 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 2 : Vector<Int> reqMSfld(1,-1); // assume all, indescriminately
1042 2 : 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 2 : Vector<Int> reqMSint(1,-1); // assume all
1049 2 : if (cls.ent.length()>0)
1050 0 : reqMSint.reference(this->getMSuniqueIds(clsms,"intent"));
1051 : //cout << "reqMSint = " << reqMSint << endl;
1052 2 : 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 2 : 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 2 : Vector<Int> reqMSant(nMSAnt_);
1063 2 : 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 2 : 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 2 : 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 4 : for (uInt iCTobs=0;iCTobs<reqCTobs.nelements();++iCTobs) {
1079 2 : 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 2 : Vector<Int> theseMSobs=cls.obsmap.msids(thisCTobs,reqMSobs);
1084 2 : if (theseMSobs.nelements()==1 && theseMSobs[0]<0)
1085 2 : 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 2 : 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 2 : 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 2 : NewCalTable scanselCT(obsselCT);
1108 4 : for (uInt iCTscan=0;iCTscan<reqCTscan.nelements();++iCTscan) {
1109 2 : 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 2 : Vector<Int> theseMSscan=cls.scanmap.msids(thisCTscan,reqMSscan);
1114 2 : if (theseMSscan.nelements()==1 && theseMSscan[0]<0)
1115 2 : 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 2 : 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 2 : Vector<Int> reqCTfld=cls.fldmap.ctids(reqMSfld);
1135 : //cout << " reqCTfld = " << reqCTfld << endl;
1136 :
1137 : // For each required CT fld:
1138 2 : NewCalTable fldselCT(scanselCT);
1139 4 : for (uInt iCTfld=0;iCTfld<reqCTfld.nelements();++iCTfld) {
1140 2 : 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 2 : Vector<int> theseMSfld=cls.fldmap.msids(thisCTfld,reqMSfld);
1145 2 : if (theseMSfld.nelements()==1 && theseMSfld[0]<0)
1146 2 : theseMSfld.reference(reqMSfld);
1147 :
1148 : //cout << " thisCTfld = " << thisCTfld << ": theseMSfld = " << theseMSfld << endl;
1149 :
1150 : // Apply thisCTfld selection to the CT
1151 : try {
1152 2 : 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 2 : NewCalTable spwselCT(fldselCT);
1166 10 : for (uInt iMSspw=0;iMSspw<reqMSspw.nelements();++iMSspw) {
1167 8 : Int& thisMSspw=reqMSspw(iMSspw);
1168 8 : Int thisCTspw=cls.spwmap(thisMSspw);
1169 8 : 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 8 : 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 8 : CTCalPatchKey iclTres(icls,thisCTobs,thisCTscan,thisCTfld,thisMSspw,-1);
1194 8 : clTres_[iclTres]=CLPPResult(nPar_,nFPar_,nChanIn_[thisCTspw],nMSElem_);
1195 :
1196 8 : NewCalTable antselCT(spwselCT);
1197 8 : Bool doLinkResults(true); // initialize true, so it will happen if relevant code reached
1198 96 : for (uInt iMSant=0;iMSant<reqMSant.nelements();++iMSant) {
1199 88 : Int& thisMSant=reqMSant(iMSant);
1200 88 : Int thisCTant=cls.antmap(thisMSant);
1201 88 : if (thisCTant<0) thisCTant=thisMSant;
1202 :
1203 : // Apply thisCTant selection to CT
1204 : try {
1205 88 : 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 88 : CTCalPatchKey ici0(icls,thisCTobs,thisCTscan,thisCTfld,thisCTspw,thisCTant); // all CT indices
1218 88 : 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 88 : if (ci_.count(ici1)<1) {
1222 88 : ciname_[ici1]=ici0.print()+" rows="+String::toString(antselCT.nrow());
1223 88 : Array<Float> r(clTres_[iclTres].result(thisMSant));
1224 88 : Array<Bool> rf(clTres_[iclTres].resultFlag(thisMSant));
1225 88 : ci_[ici1]=(*cttifactoryptr_)(antselCT,cls.tinterp,r,rf);
1226 : //if (iMSant==0) cout << " Creating: CT("<<ici1.print() << ") --> CT(" << ici0.print() << ") (all antennas)" << endl;
1227 88 : }
1228 : else
1229 0 : throw(AipsError("Attempted duplicate CTCalPatchKey!"));
1230 :
1231 : // Now distribute this CTTimeInterp1 instance to all relevant MS indices
1232 176 : for (uInt iMSobs=0;iMSobs<theseMSobs.nelements();++iMSobs) {
1233 88 : Int& thisMSobs=theseMSobs(iMSobs);
1234 176 : for (uInt iMSscan=0;iMSscan<theseMSscan.nelements();++iMSscan) {
1235 88 : Int& thisMSscan=theseMSscan(iMSscan);
1236 176 : for (uInt iMSfld=0;iMSfld<theseMSfld.nelements();++iMSfld) {
1237 88 : Int& thisMSfld=theseMSfld(iMSfld);
1238 176 : for (uInt iMSint=0;iMSint<theseMSint.nelements();++iMSint) {
1239 88 : Int& thisMSint=theseMSint(iMSint);
1240 :
1241 88 : MSCalPatchKey ims(thisMSobs,thisMSscan,thisMSfld,thisMSint,thisMSspw,thisMSant);
1242 88 : if (msci_.count(ims)<1) {
1243 88 : msciname_[ims]=ciname_[ici1];
1244 88 : msci_[ims]=ci_[ici1];
1245 : }
1246 : else
1247 0 : throw(AipsError("Attempted duplicate MSCalPatchKey!"));
1248 :
1249 : // Link these obs,scan,fld,ant,spw to the correct results object
1250 : // (as a group over antennas; should move this out of ant loop, really)
1251 88 : if (doLinkResults) {
1252 8 : MSCalPatchKey imsgroup(thisMSobs,thisMSscan,thisMSfld,thisMSint,thisMSspw,-1);
1253 8 : msTres_[imsgroup]=clTres_[iclTres];
1254 8 : msFres_[imsgroup]=CLPPResult(); // this will be resized on-demand
1255 8 : ctspw_[imsgroup]=thisCTspw;
1256 8 : finterp_[imsgroup]=cls.finterp;
1257 8 : }
1258 88 : } // iMSint
1259 : } // iMSfld
1260 : } // iMSscan
1261 : } // iMSobs
1262 88 : doLinkResults = False; // Don't do it again
1263 88 : } // iMSant
1264 8 : } // iMSspw
1265 2 : } // iCTfld
1266 2 : } // iCTscan
1267 2 : } // iCTobs
1268 :
1269 :
1270 6 : } // icls
1271 :
1272 2 : } // ctor
1273 :
1274 0 : void CLPatchPanel::recordBadMSIndices(const Vector<Int>& obs, const Vector<Int>& scan,
1275 : const Vector<Int>& fld,
1276 : const Vector<Int>& ent, const Vector<Int>& spw) {
1277 :
1278 :
1279 : // Record bad _MS_ indices
1280 0 : for (uInt iobs=0;iobs<obs.nelements();++iobs) {
1281 0 : for (uInt iscan=0;iscan<scan.nelements();++iscan) {
1282 0 : for (uInt ifld=0;ifld<fld.nelements();++ifld) {
1283 0 : for (uInt ient=0;ient<ent.nelements();++ient) {
1284 0 : for (uInt ispw=0;ispw<spw.nelements();++ispw) {
1285 :
1286 0 : MSCalPatchKey ims(obs[iobs],scan[iscan],fld[ifld],ent[ient],spw[ispw],-1); // All ants
1287 0 : if (badmsciname_.count(ims)<1) {
1288 0 : badmsciname_[ims]=ims.print();
1289 : //cout << " Bad MS indices: " << ims.print() << endl;
1290 : }
1291 0 : }
1292 : }
1293 : }
1294 : }
1295 : }
1296 0 : return;
1297 : }
1298 :
1299 :
1300 96 : void CLPatchPanel::selectOnCT(NewCalTable& ctout,const NewCalTable& ctin,
1301 : const String& obs, const String& scan,
1302 : const String& fld,
1303 : const String& spw, const String& ant1) {
1304 :
1305 96 : String taql("");
1306 96 : if (ant1.length()>0)
1307 88 : taql="ANTENNA1=="+ant1;
1308 :
1309 : // cout << "taql = >>" << taql << "<<" << endl;
1310 :
1311 : // Forward to generic method (sans intent)
1312 96 : CTInterface cti(ctin);
1313 96 : this->selectOnCTorMS(ctout,cti,obs,scan,fld,"",spw,"",taql);
1314 :
1315 96 : }
1316 :
1317 2 : void CLPatchPanel::selectOnMS(MeasurementSet& msout,const MeasurementSet& msin,
1318 : const String& obs, const String& scan,
1319 : const String& fld, const String& ent,
1320 : const String& spw, const String& ant) {
1321 :
1322 : // Forward to generic method
1323 2 : MSInterface msi(msin);
1324 2 : this->selectOnCTorMS(msout,msi,obs,scan,fld,ent,spw,ant,"");
1325 :
1326 2 : }
1327 98 : void CLPatchPanel::selectOnCTorMS(Table& ctout,MSSelectableTable& msst,
1328 : const String& obs, const String& scan,
1329 : const String& fld, const String& ent,
1330 : const String& spw, const String& ant,
1331 : const String& taql) {
1332 :
1333 98 : MSSelection mss;
1334 98 : if (obs.length()>0)
1335 0 : mss.setObservationExpr(obs);
1336 98 : if (scan.length()>0)
1337 0 : mss.setScanExpr(scan);
1338 98 : if (fld.length()>0)
1339 0 : mss.setFieldExpr(fld);
1340 98 : if (ent.length()>0)
1341 0 : mss.setStateExpr(ent);
1342 98 : if (spw.length()>0)
1343 8 : mss.setSpwExpr(spw);
1344 98 : if (ant.length()>0)
1345 0 : mss.setAntennaExpr(ant); // this will not behave as required for Jones caltables (its bl-based selection!)
1346 98 : if (taql.length()>0)
1347 88 : mss.setTaQLExpr(taql);
1348 :
1349 :
1350 98 : TableExprNode ten=mss.toTableExprNode(&msst);
1351 98 : getSelectedTable(ctout,*(msst.table()),ten,"");
1352 :
1353 98 : }
1354 :
1355 0 : Vector<Int> CLPatchPanel::getCLuniqueIds(NewCalTable& ct, String vcol) {
1356 :
1357 0 : ROCTMainColumns ctmc(ct);
1358 :
1359 : // Extract the requested column as a Vector
1360 0 : Vector<Int> colv;
1361 0 : if (vcol=="obs")
1362 0 : ctmc.obsId().getColumn(colv);
1363 0 : else if (vcol=="scan")
1364 0 : ctmc.scanNo().getColumn(colv);
1365 0 : else if (vcol=="fld")
1366 0 : ctmc.fieldId().getColumn(colv);
1367 0 : else if (vcol=="spw")
1368 0 : ctmc.spwId().getColumn(colv);
1369 : else
1370 0 : throw(AipsError("Unsupported column in CLPatchPanel::getCLuniqueIds"));
1371 :
1372 : // Reduce to a unique list
1373 0 : uInt nuniq=genSort(colv,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
1374 0 : colv.resize(nuniq,true);
1375 :
1376 0 : return colv;
1377 :
1378 0 : }
1379 :
1380 2 : Vector<Int> CLPatchPanel::getMSuniqueIds(MeasurementSet& ms, String which) {
1381 :
1382 2 : MSColumns msc(ms);
1383 :
1384 : // Extract the requested column as a Vector
1385 2 : Vector<Int> colv;
1386 2 : if (which=="obs")
1387 0 : msc.observationId().getColumn(colv);
1388 2 : else if (which=="scan")
1389 0 : msc.scanNumber().getColumn(colv);
1390 2 : else if (which=="fld")
1391 0 : msc.fieldId().getColumn(colv);
1392 2 : else if (which=="intent")
1393 0 : msc.stateId().getColumn(colv);
1394 2 : else if (which=="spw")
1395 2 : msc.dataDescId().getColumn(colv); // these are actually ddids!
1396 : else
1397 0 : throw(AipsError("Unsupported column in CLPatchPanel::getCLuniqueIds"));
1398 :
1399 : // Reduce to a unique list
1400 2 : uInt nuniq=genSort(colv,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
1401 2 : colv.resize(nuniq,true);
1402 :
1403 : // In case of spw, translate from ddid
1404 2 : if (which=="spw") {
1405 2 : Vector<Int> spwids;
1406 2 : msc.dataDescription().spectralWindowId().getColumn(spwids);
1407 10 : for (uInt i=0;i<colv.nelements();++i)
1408 8 : colv[i]=spwids[colv[i]];
1409 2 : }
1410 :
1411 : // return the value
1412 4 : return colv;
1413 :
1414 2 : }
1415 :
1416 :
1417 :
1418 : // Destructor
1419 4 : CLPatchPanel::~CLPatchPanel() {
1420 :
1421 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::~CLPatchPanel()" << endl;
1422 :
1423 : // If we encountered data (via calAvailable) for which no calibration was arranged,
1424 : // report this with a warning
1425 2 : if (!passcal_.empty()) {
1426 0 : logsink_ << " " << endl << LogIO::POST;
1427 0 : logsink_ << LogIO::WARN << "The cal library instructions for caltable = " << endl;
1428 0 : logsink_ << ct_.keywordSet().asString("VisCal") << ": " << Path(ct_.tableName()).baseName().before(".tempMem") << endl;
1429 0 : logsink_ << "...excluded calibration instructions for the following joint data selection(s):" << LogIO::POST;
1430 0 : for (std::map<MSCalPatchKey,String>::iterator it=passcal_.begin(); it!=passcal_.end(); ++it)
1431 0 : logsink_ << LogIO::WARN << " " << it->second << LogIO::POST;
1432 0 : logsink_ << LogIO::WARN << "Therefore, these data were passed uncalibrated (and unflagged) by *this* caltable." << endl;
1433 0 : logsink_ << "Please check your results carefully to be sure this was intended." << LogIO::POST;
1434 : }
1435 :
1436 : // Delete the atomic interpolators
1437 90 : for (std::map<CTCalPatchKey,CTTimeInterp1*>::iterator it=ci_.begin(); it!=ci_.end(); ++it)
1438 88 : delete it->second;
1439 :
1440 4 : }
1441 :
1442 : // Is specific calibration explicitly available for a obs,scan,fld,intent,spw,ant combination?
1443 2400 : Bool CLPatchPanel::calAvailable(casacore::Int msobs, casacore::Int msscan, casacore::Int msfld,
1444 : casacore::Int msent, casacore::Int msspw, casacore::Int msant) {
1445 :
1446 2400 : const MSCalPatchKey key(msobs,msscan,msfld,msent,msspw,msant);
1447 :
1448 2400 : Bool avail=msTres_.count(key)>0;
1449 :
1450 : // Record lack of available calibration for this specific data selection,
1451 : // so we can report it at the end.
1452 : // Such data is passed unchanged by this caltable
1453 2400 : if (!avail) {
1454 0 : if (passcal_.count(key)<1) {
1455 : //cout << "passcal_.size() = " << passcal_.size() << "--->";
1456 0 : passcal_[key]=key.print();
1457 : //cout << passcal_.size() << endl;
1458 : }
1459 : }
1460 :
1461 2400 : return avail;
1462 :
1463 2400 : }
1464 :
1465 : // The specified MS indices are "OK" if not recorded as expected-but-absent as per
1466 : // the callibrary specification. Expected-but-absent occurs when a CL entry
1467 : // indicates (even implicitly) that the MS index combination is supported but
1468 : // the required CalTable indices (via *map params) are not actually available
1469 : // in the CalTable, e.g., a missing spw.
1470 : // Such data cannot be calibrated and this CL cannot be agnostic about it...
1471 : // Note that MSIndicesOK can return true for MS index combinations for which
1472 : // calibration is not actually available, if the CL entry does not purport
1473 : // to support calibrating them. In such cases, calAvailable() returns false, and
1474 : // this CL is agnostic w.r.t. such data and lets it pass thru unchanged.
1475 : // Returns TRUE when the specified MS indices are calibrateable or passable.
1476 2400 : Bool CLPatchPanel::MSIndicesOK(casacore::Int msobs, casacore::Int msscan, casacore::Int msfld,
1477 : casacore::Int msent, casacore::Int msspw, casacore::Int msant) {
1478 :
1479 2400 : const MSCalPatchKey key(msobs,msscan,msfld,msent,msspw,msant);
1480 2400 : Bool bad=badmsciname_.count(key)>0;
1481 : // if (bad) {
1482 : // cout << Path(ct_.tableName()).baseName().before(".tempMem") << " should but can't calibrate: " << key.print() << endl;
1483 : // }
1484 :
1485 : // Return TRUE if NOT bad
1486 2400 : return !bad;
1487 :
1488 2400 : }
1489 :
1490 0 : Bool CLPatchPanel::interpolate(Cube<Complex>& resultC, Cube<Bool>& resFlag,
1491 : Int msobs, Int msscan, Int msfld, Int msent, Int msspw,
1492 : Double time, Double freq) {
1493 :
1494 0 : Bool newcal(false);
1495 :
1496 : // resultC and resFlag will be unchanged if newcal remains false
1497 0 : Cube<Float> f; // temporary to reference Float interpolation result
1498 0 : newcal=interpolate(f,resFlag,msobs,msscan,msfld,msent,msspw,time,freq);
1499 0 : if (newcal)
1500 : // convert to complex and have resultC take over ownership
1501 0 : resultC.reference(RIorAPArray(f).c());
1502 :
1503 0 : return newcal; // If T, calling scope should act nontrivially, if necessary
1504 :
1505 0 : }
1506 :
1507 0 : Bool CLPatchPanel::interpolate(Cube<Complex>& resultC, Cube<Bool>& resFlag,
1508 : Int msobs, Int msscan, Int msfld, Int msent, Int msspw,
1509 : Double time, const Vector<Double>& freq) {
1510 :
1511 0 : Bool newcal(false);
1512 :
1513 : // resultC and resFlag will be unchanged if newcal remains false
1514 0 : Cube<Float> f; // temporary to reference Float interpolation result
1515 0 : newcal=interpolate(f,resFlag,msobs,msscan,msfld,msent,msspw,time,freq);
1516 :
1517 0 : if (newcal)
1518 : // convert to complex and have resultC take over ownership
1519 0 : resultC.reference(RIorAPArray(f).c());
1520 :
1521 0 : return newcal; // If T, calling scope should act nontrivially, if necessary
1522 :
1523 0 : }
1524 :
1525 :
1526 2400 : Bool CLPatchPanel::interpolate(Cube<Float>& resultR, Cube<Bool>& resFlag,
1527 : Int msobs, Int msscan, Int msfld, Int msent, Int msspw,
1528 : Double time, Double freq) {
1529 :
1530 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::interpolate(...)" << endl;
1531 :
1532 : // NB: resultR and resFlag will be unchanged if newcal remains false
1533 2400 : Bool newcal(false);
1534 :
1535 : // Suppled arrays reference the result (if available)
1536 2400 : MSCalPatchKey ires(msobs,msscan,msfld,msent,msspw,-1);
1537 :
1538 : // Trap lack of available calibration for requested obs,fld,intent,spw
1539 2400 : if (msTres_.count(ires)==0) {
1540 0 : throw(AipsError("No calibration arranged for "+ires.print()+
1541 0 : " in callib for caltable="+
1542 0 : Path(ct_.tableName()).baseName().before(".tempMemCalTable") ));
1543 : }
1544 :
1545 : // If result_ is at a new address (cf last time, this spw), treat as new
1546 : // TBD: do this is a less obscure way... (e.g., invent the CLCalGroup(nant))
1547 2400 : if (lastresadd_(msspw)!=msTres_[ires].result_.data())
1548 8 : newcal=true;
1549 2400 : lastresadd_(msspw)=msTres_[ires].result_.data();
1550 :
1551 : // Loop over _output_ elements
1552 28800 : for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1553 : // Call fully _patched_ time-interpolator, keeping track of 'newness'
1554 : // fills ctTres_ implicitly
1555 26400 : MSCalPatchKey ims(msobs,msscan,msfld,msent,msspw,iMSElem);
1556 26400 : if (msci_.count(ims)>0) {
1557 26400 : if (freq>0.0)
1558 0 : newcal|=msci_[ims]->interpolate(time,freq);
1559 : else
1560 26400 : newcal|=msci_[ims]->interpolate(time);
1561 : }
1562 26400 : }
1563 :
1564 2400 : if (newcal) {
1565 1208 : resultR.reference(msTres_[ires].result_);
1566 1208 : resFlag.reference(msTres_[ires].resultFlag_);
1567 : }
1568 :
1569 2400 : return newcal; // If true, calling scope should act
1570 2400 : }
1571 :
1572 :
1573 0 : Bool CLPatchPanel::interpolate(Cube<Float>& resultR, Cube<Bool>& resFlag,
1574 : Int msobs, Int msscan, Int msfld, Int msent, Int msspw,
1575 : Double time, const Vector<Double>& freq) {
1576 :
1577 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::interpolate(Cube<F>,...,Vector<D>)" << endl;
1578 :
1579 : // NB: resultR and resFlag will be unchanged if newcal remains false
1580 0 : Bool newcal(false);
1581 :
1582 : // Suppled arrays reference the result (if available)
1583 0 : MSCalPatchKey ires(msobs,msscan,msfld,msent,msspw,-1);
1584 :
1585 : // Trap lack of available calibration for requested obs,scan,fld,intent,spw
1586 0 : if (msTres_.count(ires)==0) {
1587 0 : throw(AipsError("No calibration arranged for "+ires.print()+
1588 0 : " in callib for caltable="+
1589 0 : Path(ct_.tableName()).baseName().before(".tempMemCalTable") ));
1590 : }
1591 :
1592 : // If result_ is at a new address (cf last time, this msspw), treat as new
1593 : // TBD: do this is a less obscure way...
1594 0 : if (lastresadd_(msspw)!=msTres_[ires].result_.data())
1595 0 : newcal=true;
1596 0 : lastresadd_(msspw)=msTres_[ires].result_.data();
1597 :
1598 : // Sometimes we need to force the freq interp, even if the time-interp isn't new
1599 0 : Bool forceFinterp=false || newcal;
1600 :
1601 : // The follow occurs unnecessarily in mosaics when newcal=false (i.e., when resampleInFreq won't be called below)
1602 0 : uInt nMSChan=freq.nelements(); // The number of requested channels
1603 0 : if (msFres_.count(ires)>0) {
1604 0 : if (msFres_[ires].result_.nelements()==0 ||
1605 0 : msFres_[ires].result(0).ncolumn()!=nMSChan) {
1606 0 : msFres_[ires].resize(nPar_,nFPar_,nMSChan,nMSElem_);
1607 : }
1608 : }
1609 : else
1610 0 : throw(AipsError("Error finding msFres_ in CLPatchPanel::interpolate(C<F>,...,V<D>)"));
1611 :
1612 : // Loop over _output_ antennas
1613 0 : for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1614 : // Call time interpolation calculation; resample in freq if new
1615 : // (fills msTRes_ implicitly)
1616 :
1617 0 : MSCalPatchKey ims(msobs,msscan,msfld,msent,msspw,iMSElem);
1618 0 : if (msci_.count(ims)>0) {
1619 0 : if (msci_[ims]->interpolate(time) || forceFinterp) {
1620 :
1621 : // Resample in frequency
1622 0 : Matrix<Float> fR( msFres_[ires].result(iMSElem) );
1623 0 : Matrix<Bool> fRflg( msFres_[ires].resultFlag(iMSElem) );
1624 0 : Matrix<Float> tR( msTres_[ires].result(iMSElem) );
1625 0 : Matrix<Bool> tRflg( msTres_[ires].resultFlag(iMSElem) );
1626 0 : resampleInFreq(fR,fRflg,freq,tR,tRflg,freqIn_(ctspw_[ires]),finterp_[ires]);
1627 :
1628 : // Calibration is new
1629 0 : newcal=true;
1630 0 : }
1631 : }
1632 0 : }
1633 :
1634 0 : if (newcal) {
1635 : // Supplied arrays to reference the result
1636 0 : resultR.reference(msFres_[ires].result_);
1637 0 : resFlag.reference(msFres_[ires].resultFlag_);
1638 : }
1639 :
1640 0 : return newcal;
1641 0 : }
1642 :
1643 0 : Bool CLPatchPanel::getTresult(Cube<Float>& resultR, Cube<Bool>& resFlag,
1644 : Int obs, Int scan, Int fld, Int ent, Int spw) {
1645 :
1646 0 : MSCalPatchKey mskey(obs,scan,fld,ent,spw,-1);
1647 :
1648 0 : if (msTres_.count(mskey)==0)
1649 0 : throw(AipsError("(CLPP::getTresult) No calibration arranged for "+mskey.print()+
1650 0 : " in callib for caltable="+
1651 0 : Path(ct_.tableName()).baseName().before(".tempMemCalTable") ));
1652 :
1653 : // Reference the requested Cube
1654 0 : resultR.reference(msTres_[mskey].result_);
1655 0 : resFlag.reference(msTres_[mskey].resultFlag_);
1656 :
1657 0 : return true;
1658 :
1659 0 : }
1660 :
1661 :
1662 :
1663 :
1664 0 : void CLPatchPanel::listmappings() {
1665 :
1666 0 : cout << "CalTable: " << ct_.tableName() << endl;
1667 0 : cout << "There are " << ci_.size() << " cal interpolation engines." << endl;
1668 0 : cout << "There are " << msci_.size() << " unique MS id combinations mapped to them:" << endl;
1669 0 : for (std::map<MSCalPatchKey,String>::iterator it=msciname_.begin(); it!=msciname_.end(); ++it)
1670 0 : cout << "MS (" << it->first.print() << ") --> CL (" << it->second << ")" << endl;
1671 :
1672 0 : cout << endl << "There are " << badmsciname_.size() << " expected but ABSENT MS id combinations (all ants):" << endl;
1673 0 : for (std::map<MSCalPatchKey,String>::iterator it=badmsciname_.begin(); it!=badmsciname_.end(); ++it)
1674 0 : cout << "MS (" << it->first.print() << ")" << endl;
1675 :
1676 0 : cout << endl;
1677 :
1678 0 : }
1679 :
1680 :
1681 :
1682 :
1683 : // Report state
1684 0 : void CLPatchPanel::state() {
1685 :
1686 : if (CTPATCHPANELVERB) cout << "CLPatchPanel::state()" << endl;
1687 :
1688 0 : cout << "-state--------" << endl;
1689 0 : cout << " ct_ = " << ct_.tableName() << endl;
1690 0 : cout << boolalpha;
1691 0 : cout << " isCmplx_ = " << isCmplx_ << endl;
1692 0 : cout << " nPar_ = " << nPar_ << endl;
1693 0 : cout << " nFPar_ = " << nFPar_ << endl;
1694 0 : cout << " nMSObs_ = " << nMSObs_ << endl;
1695 0 : cout << " nMSFld_ = " << nMSFld_ << endl;
1696 0 : cout << " nMSSpw_ = " << nMSSpw_ << endl;
1697 0 : cout << " nMSAnt_ = " << nMSAnt_ << endl;
1698 0 : cout << " nMSElem_ = " << nMSElem_ << endl;
1699 0 : cout << " nCTObs_ = " << nCTObs_ << endl;
1700 0 : cout << " nCTFld_ = " << nCTFld_ << endl;
1701 0 : cout << " nCTSpw_ = " << nCTSpw_ << endl;
1702 0 : cout << " nCTAnt_ = " << nCTAnt_ << endl;
1703 0 : cout << " nCTElem_ = " << nCTElem_ << endl;
1704 :
1705 : // cout << " timeType_ = " << timeType_ << endl;
1706 : // cout << " freqType_ = " << freqType_ << endl;
1707 0 : }
1708 :
1709 : // Resample in frequency
1710 0 : void CLPatchPanel::resampleInFreq(Matrix<Float>& fres,Matrix<Bool>& fflg,const Vector<Double>& fout,
1711 : Matrix<Float>& tres,Matrix<Bool>& tflg,const Vector<Double>& fin,
1712 : String finterp) {
1713 :
1714 : if (CTPATCHPANELVERB) cout << " CLPatchPanel::resampleInFreq(...)" << endl;
1715 :
1716 : // if no good solutions coming in, return flagged
1717 0 : if (nfalse(tflg)==0) {
1718 0 : fflg.set(true);
1719 0 : return;
1720 : }
1721 :
1722 0 : Int flparmod=nFPar_/nPar_; // for indexing the flag Matrices on the par axis
1723 :
1724 0 : Bool unWrapPhase=flparmod>1;
1725 :
1726 : // cout << "nFPar_,nPar_,flparmod = " << nFPar_ << "," << nPar_ << "," << flparmod << endl;
1727 :
1728 0 : fres=0.0;
1729 :
1730 0 : for (Int ifpar=0;ifpar<nFPar_;++ifpar) {
1731 :
1732 : // Slice by par (each Matrix row)
1733 0 : Vector<Float> fresi(fres.row(ifpar)), tresi(tres.row(ifpar));
1734 0 : Vector<Bool> fflgi(fflg.row(ifpar/flparmod)), tflgi(tflg.row(ifpar/flparmod));
1735 :
1736 : // Mask time result by flags
1737 0 : Vector<Double> mfin=fin(!tflgi).getCompressedArray();
1738 :
1739 0 : if (mfin.nelements()==0) {
1740 : // cout << ifpar << " All chans flagged!" << endl;
1741 : // Everything flagged this par
1742 : // Just flag, zero and go on to the next one
1743 0 : fflgi.set(true);
1744 0 : fresi.set(0.0);
1745 0 : continue;
1746 : }
1747 :
1748 0 : mfin/=1.0e9; // in GHz
1749 0 : Vector<Float> mtresi=tresi(!tflgi).getCompressedArray();
1750 :
1751 : // Trap case of same in/out frequencies
1752 0 : if (fout.nelements()==mfin.nelements() && allNear(fout,mfin,1.e-10)) {
1753 : // Just copy
1754 0 : fresi=mtresi;
1755 0 : fflgi.set(false); // none are flagged
1756 0 : continue;
1757 : }
1758 :
1759 0 : if (ifpar%2==1 && unWrapPhase) {
1760 0 : for (uInt i=1;i<mtresi.nelements();++i) {
1761 0 : while ( (mtresi(i)-mtresi(i-1))>M_PI ) mtresi(i)-=(2.0*M_PI);
1762 0 : while ( (mtresi(i)-mtresi(i-1))<-M_PI ) mtresi(i)+=(2.0*M_PI);
1763 : }
1764 : }
1765 :
1766 :
1767 : // TBD: ensure phases tracked on freq axis...
1768 :
1769 : // TBD: handle flags carefully!
1770 : // (i.e., flag gaps larger than user's "freach")
1771 : // For now,just unset them
1772 0 : fflgi.set(false);
1773 :
1774 :
1775 : // Set flags carefully
1776 0 : resampleFlagsInFreq(fflgi,fout,tflgi,fin,finterp);
1777 :
1778 :
1779 : // Always use nearest on edges
1780 : // TBD: trap cases where frequencies don't overlap at all
1781 : // (fout(hi)<mfin(0) || fout(lo)> mfin(ihi))..... already OK (lo>hi)?
1782 : // TBD: optimize the following by forming Slices in the
1783 : // while loops and doing Array assignment once afterwords
1784 :
1785 0 : Int nfreq=fout.nelements();
1786 0 : Int lo=0;
1787 0 : Int hi=fresi.nelements()-1;
1788 0 : Double inlo(mfin(0));
1789 0 : Int ihi=mtresi.nelements()-1;
1790 0 : Double inhi(mfin(ihi));
1791 :
1792 : // Handle 'nearest' extrapolation in sideband-dep way
1793 0 : Bool inUSB(inhi>inlo);
1794 0 : Bool outUSB(fout(hi)>fout(lo));
1795 0 : if (inUSB) {
1796 0 : if (outUSB) {
1797 0 : while (lo<nfreq && fout(lo)<=inlo) fresi(lo++)=mtresi(0);
1798 0 : while (hi>-1 && fout(hi)>=inhi) fresi(hi--)=mtresi(ihi);
1799 : }
1800 : else { // "outLSB"
1801 0 : while (lo<nfreq && fout(lo)>=inhi) fresi(lo++)=mtresi(ihi);
1802 0 : while (hi>-1 && fout(hi)<=inlo) fresi(hi--)=mtresi(0);
1803 : }
1804 : }
1805 : else { // "inLSB"
1806 0 : if (outUSB) {
1807 0 : while (lo<nfreq && fout(lo)<=inhi) fresi(lo++)=mtresi(ihi);
1808 0 : while (hi>-1 && fout(hi)>=inlo) fresi(hi--)=mtresi(0);
1809 : }
1810 : else { // "outLSB"
1811 0 : while (lo<nfreq && fout(lo)>=inlo) fresi(lo++)=mtresi(0);
1812 0 : while (hi>-1 && fout(hi)<=inhi) fresi(hi--)=mtresi(ihi);
1813 : }
1814 : }
1815 :
1816 : // cout << "lo, hi = " << lo << ","<<hi << endl;
1817 :
1818 0 : if (lo>hi) continue; // Frequencies didn't overlap, nearest was used
1819 :
1820 : // Use InterpolateArray1D to fill in the middle
1821 0 : IPosition blc(1,lo), trc(1,hi);
1822 0 : Vector<Float> slfresi(fresi(blc,trc));
1823 0 : Vector<Double> slfout(fout(blc,trc));
1824 :
1825 0 : InterpolateArray1D<Double,Float>::interpolate(slfresi,slfout,mfin,mtresi,this->ftype(finterp));
1826 :
1827 0 : }
1828 : }
1829 :
1830 :
1831 :
1832 0 : void CLPatchPanel::resampleFlagsInFreq(Vector<Bool>& flgout,const Vector<Double>& fout,
1833 : Vector<Bool>& flgin,const Vector<Double>& fin,
1834 : String finterp) {
1835 :
1836 : // cout << "resampleFlagsInFreq" << endl;
1837 :
1838 : #define NEAREST InterpolateArray1D<Double,Float>::nearestNeighbour
1839 : #define LINEAR InterpolateArray1D<Double,Float>::linear
1840 : #define CUBIC InterpolateArray1D<Double,Float>::cubic
1841 : #define SPLINE InterpolateArray1D<Double,Float>::spline
1842 :
1843 :
1844 : // Handle chan-dep flags
1845 0 : if (finterp.contains("flag")) {
1846 :
1847 0 : InterpolateArray1D<Double,Float>::InterpolationMethod interpMeth=this->ftype(finterp);
1848 :
1849 0 : Vector<Double> finGHz=fin/1e9;
1850 :
1851 : // Determine implied mode-dep flags indexed by channel registration
1852 0 : uInt nflg=flgin.nelements();
1853 0 : Vector<Bool> flreg(nflg,false);
1854 0 : switch (interpMeth) {
1855 0 : case NEAREST: {
1856 : // Just use input flags
1857 0 : flreg.reference(flgin);
1858 0 : break;
1859 : }
1860 0 : case LINEAR: {
1861 0 : for (uInt i=0;i<nflg-1;++i)
1862 0 : flreg[i]=(flgin[i] || flgin[i+1]);
1863 0 : flreg[nflg-1]=flreg[nflg-2];
1864 0 : break;
1865 : }
1866 0 : case CUBIC:
1867 : case SPLINE: {
1868 0 : for (uInt i=1;i<nflg-2;++i)
1869 0 : flreg[i]=(flgin[i-1] || flgin[i] || flgin[i+1] || flgin[i+2]);
1870 0 : flreg[0]=flreg[1];
1871 0 : flreg[nflg-2]=flreg[nflg-3];
1872 0 : flreg[nflg-1]=flreg[nflg-3];
1873 0 : break;
1874 : }
1875 : }
1876 :
1877 :
1878 :
1879 : // Now step through requested chans, setting flags
1880 0 : uInt ireg=0;
1881 0 : uInt nflgout=flgout.nelements();
1882 0 : for (uInt iflgout=0;iflgout<nflgout;++iflgout) {
1883 :
1884 : // Find nominal registration (the _index_ just left)
1885 0 : Bool exact(false);
1886 0 : ireg=binarySearch(exact,finGHz,fout(iflgout),nflg,0);
1887 :
1888 : // If registration is exact, assign verbatim
1889 : // NB: the calibration value calculation occurs agnostically w.r.t. flags,
1890 : // so the calculated value should also match
1891 : // TBD: Add "|| near(finGHz[ireg],fout(iflgout),1e-10) in case "exact" has
1892 : // precision issues?
1893 0 : if (exact) {
1894 0 : flgout[iflgout]=flgin[ireg];
1895 0 : continue;
1896 : }
1897 :
1898 : // Not exact, so carefully handle bracketing
1899 0 : if (ireg>0)
1900 0 : ireg-=1;
1901 0 : ireg=min(ireg,nflg-1);
1902 :
1903 : //while (finGHz(ireg)<=fout(iflgout) && ireg<nflg-1) {
1904 : // ireg+=1; // USB specific!
1905 : //}
1906 : //if (ireg>0 && finGHz(ireg)!=fout(iflgout)) --ireg; // registration is one sample prior
1907 :
1908 : // refine registration by interp type
1909 0 : switch (interpMeth) {
1910 0 : case NEAREST: {
1911 : // nearest might be forward sample
1912 0 : if ( ireg<(nflg-1) &&
1913 0 : abs(fout[iflgout]-finGHz[ireg])>abs(finGHz[ireg+1]-fout[iflgout]) )
1914 0 : ireg+=1;
1915 0 : break;
1916 : }
1917 0 : case LINEAR: {
1918 0 : if (ireg==(nflg-1)) // need one more sample to the right
1919 0 : ireg-=1;
1920 0 : break;
1921 : }
1922 0 : case CUBIC:
1923 : case SPLINE: {
1924 0 : if (ireg==0) ireg+=1; // need one more sample to the left
1925 0 : if (ireg>(nflg-3)) ireg=nflg-3; // need two more samples to the right
1926 0 : break;
1927 : }
1928 : }
1929 :
1930 : // Assign effective flag
1931 0 : flgout[iflgout]=flreg[ireg];
1932 :
1933 : /*
1934 : cout << iflgout << " "
1935 : << ireg << " "
1936 : << flreg[ireg]
1937 : << endl;
1938 : */
1939 :
1940 : }
1941 :
1942 0 : }
1943 : else
1944 : // We are interp/extrap-olating gaps absolutely
1945 0 : flgout.set(false);
1946 :
1947 0 : }
1948 :
1949 0 : InterpolateArray1D<Double,Float>::InterpolationMethod CLPatchPanel::ftype(String& strtype) {
1950 0 : if (strtype.contains("nearest"))
1951 0 : return InterpolateArray1D<Double,Float>::nearestNeighbour;
1952 0 : if (strtype.contains("linear"))
1953 0 : return InterpolateArray1D<Double,Float>::linear;
1954 0 : if (strtype.contains("cubic"))
1955 0 : return InterpolateArray1D<Double,Float>::cubic;
1956 0 : if (strtype.contains("spline"))
1957 0 : return InterpolateArray1D<Double,Float>::spline;
1958 :
1959 : // cout << "Using linear for freq interpolation as last resort." << endl;
1960 0 : return InterpolateArray1D<Double,Float>::linear;
1961 :
1962 :
1963 : }
1964 :
1965 :
1966 : } //# NAMESPACE CASA - END
|
