LCOV - code coverage report
Current view: top level - synthesis/CalTables - CLPatchPanel.cc (source / functions) Hit Total Coverage
Test: casacpp_coverage.info Lines: 555 972 57.1 %
Date: 2024-11-06 17:42:47 Functions: 40 58 69.0 %

          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      417562 : CalPatchKey::CalPatchKey(IPosition keyids) :
      54      417562 :   cpk_(keyids.asVector())
      55      417562 : {}
      56             : 
      57             : // Lexographical lessthan
      58     4458288 : 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    27028800 :   for (Int i=0;i<6;++i) {
      66    25401422 :     if (cpk_[i]>-1 && other.cpk_[i]>-1 && cpk_[i]!=other.cpk_[i])
      67             :       // both non-negative and not equal, so evaluate element <
      68     2830910 :       return cpk_[i]<other.cpk_[i];
      69             :   }
      70             :   // All apparently equal so "<" is false
      71     1627378 :   return false;
      72             : }
      73             : 
      74      414410 : MSCalPatchKey::MSCalPatchKey(Int obs,Int scan,Int fld,Int ent,Int spw,Int ant) :
      75      414410 :   CalPatchKey(IPosition(6,obs,scan,fld,ent,spw,ant)),
      76      414410 :   obs_(obs),scan_(scan),fld_(fld),ent_(ent),spw_(spw),ant_(ant)
      77      414410 : {}
      78             : 
      79             : // text output
      80           1 : String MSCalPatchKey::print() const {
      81           2 :   return "obs="+(obs_<0 ? "*" : String::toString(obs_))+" "
      82           4 :     "scan="+(scan_<0 ? "*" : String::toString(scan_))+" "
      83           4 :     "fld="+(fld_<0 ? "*" : String::toString(fld_))+" "
      84           4 :     "intent="+(ent_<0 ? "*" : String::toString(ent_))+" "
      85           4 :     "spw="+(spw_<0 ? "*" : String::toString(spw_))+" "
      86           3 :     "ant="+(ant_<0 ? "*" : String::toString(ant_));
      87             : }
      88             : 
      89             : 
      90        3152 : CTCalPatchKey::CTCalPatchKey(Int clsl,Int obs,Int scan,Int fld,Int spw,Int ant) :
      91        3152 :   CalPatchKey(IPosition(6,clsl,obs,scan,fld,spw,ant)),
      92        3152 :   clsl_(clsl),obs_(obs),scan_(scan),fld_(fld),spw_(spw),ant_(ant)
      93        3152 : {}
      94             : 
      95             : // text output
      96        1538 : String CTCalPatchKey::print() const {
      97        3076 :   return "cl="+(clsl_<0 ? "*" : String::toString(clsl_))+" "
      98        6152 :     "obs="+(obs_<0 ? "*" : String::toString(obs_))+" "
      99        6152 :     "scan="+(scan_<0 ? "*" : String::toString(scan_))+" "
     100        6152 :     "fld="+(fld_<0 ? "*" : String::toString(fld_))+" "
     101        6152 :     "spw="+(spw_<0 ? "*" : String::toString(spw_))+" "
     102        4614 :     "ant="+(ant_<0 ? "*" : String::toString(ant_));
     103             : }
     104             : 
     105             : 
     106         288 : CalMap::CalMap() :
     107         288 :   vcalmap_()
     108         288 : {}
     109             : 
     110          74 : CalMap::CalMap(const Vector<Int>& calmap) :
     111          74 :   vcalmap_(calmap)
     112             : {
     113             :   //  cout << "calmap addresses: " << calmap.data() << " " << vcalmap_.data()
     114             :   //       << " nrefs= " << calmap.nrefs() << " " << vcalmap_.nrefs() << endl;
     115          74 : }
     116             : 
     117        1618 : Int CalMap::operator()(Int msid) const {
     118             :   // TBD: reconsider algorithm (maybe just return msid if over-run?)
     119        1618 :   Int ncalmap=vcalmap_.nelements();
     120        3236 :   return (msid<ncalmap ? vcalmap_(msid) :
     121        3236 :           (ncalmap>0 ? vcalmap_(ncalmap-1) : msid) ); // Avoid going off end
     122             : }
     123             : 
     124         140 : Vector<Int> CalMap::ctids(const Vector<Int>& msids) const {
     125         140 :   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         144 :   if (ncalmap<1 ||
     131           4 :       (ncalmap==1 && vcalmap_[0]<0))
     132         113 :     return Vector<Int>(1,-1);
     133             : 
     134          27 :   Vector<Bool> calmask(ncalmap,false);
     135          51 :   if (msids.nelements()==1 &&
     136          24 :       msids[0]<0)
     137             :     // MS ids indefinite, so all are ok
     138          14 :     calmask.set(true);
     139             :   else {
     140             :     // just do the ones that are requested
     141          29 :     for (uInt i=0;i<msids.nelements();++i) {
     142          16 :       const uInt& thismsid=msids(i);
     143          16 :       calmask(thismsid<ncalmap?thismsid:ncalmap-1)=true;
     144             :     }
     145             :   }
     146          27 :   Vector<Int> reqids;
     147          27 :   reqids=vcalmap_(calmask).getCompressedArray();
     148          27 :   Int nsort=genSort(reqids,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
     149          27 :   reqids.resize(nsort,true);
     150          27 :   return reqids;
     151          27 : }
     152             : 
     153         154 : Vector<Int> CalMap::msids(Int ctid,const Vector<Int>& superset) const {
     154             : 
     155         154 :   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         158 :   if (ncalmap<1 ||
     160           4 :       (ncalmap==1 && vcalmap_[0]<0) )
     161         113 :     return Vector<Int>(1,-1);
     162             : 
     163             :   // mask vcalmap_ by specified ctid...
     164          41 :   Vector<Int> msidlist(ncalmap);
     165          41 :   indgen(msidlist);
     166          41 :   Vector<Bool> msidmask(vcalmap_==ctid);
     167             :   // ...and limit to superset, if specified
     168          41 :   if (superset.nelements()>0 && superset[0]>-1)
     169          44 :     for (Int i=0;i<Int(msidmask.nelements());++i)
     170          31 :       if (!anyEQ(superset,i))
     171          15 :         msidmask[i]=false;  // exclude if not in superset
     172             : 
     173             :   // Return the ms id list
     174          82 :   return msidlist(msidmask).getCompressedArray();
     175          41 : }
     176             : 
     177             : 
     178           0 : FieldCalMap::FieldCalMap() :
     179             :   CalMap(),
     180           0 :   fieldcalmap_("")
     181           0 : {}
     182             : 
     183          19 : FieldCalMap::FieldCalMap(const Vector<Int>& calmap) :
     184             :   CalMap(calmap),
     185          19 :   fieldcalmap_("")
     186          19 : {}
     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          22 : FieldCalMap::FieldCalMap(const String fieldcalmap, 
     204             :                          const MeasurementSet& ms, const NewCalTable& ct,
     205          22 :                          String& extfldsel) :
     206             :   CalMap(),
     207          22 :   fieldcalmap_(fieldcalmap)
     208             : {
     209             : 
     210          22 :   if (fieldcalmap_=="nearest") 
     211             :     // Calculate nearest map
     212          10 :     setNearestFieldMap(ms,ct);
     213             :   else 
     214             :     // Attempt field selection
     215          12 :     setSelectedFieldMap(fieldcalmap,ms,ct,extfldsel);
     216             : 
     217          22 : }
     218             : 
     219          10 : void FieldCalMap::setNearestFieldMap(const MeasurementSet& ms, const NewCalTable& ct) {
     220             :   // Access MS and CT columns
     221          10 :   MSFieldColumns msfc(ms.field());
     222          10 :   ROCTColumns ctc(ct);
     223          10 :   setNearestFieldMap(msfc,ctc);
     224          10 : }
     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          10 : void FieldCalMap::setNearestFieldMap(const MSFieldColumns& msfc, const ROCTColumns& ctc) {
     233             : 
     234             :   // Nominally, this many field need a map
     235          10 :   Int nMSFlds=msfc.nrow();
     236          10 :   vcalmap_.resize(nMSFlds);
     237          10 :   vcalmap_.set(-1);  // no map
     238             : 
     239             :   // Discern _available_ fields in the CT
     240          10 :   Vector<Int> ctFlds;
     241          10 :   ctc.fieldId().getColumn(ctFlds);
     242          10 :   Int nAvFlds=genSort(ctFlds,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
     243          10 :   ctFlds.resize(nAvFlds,true);
     244             : 
     245             :   // If only one CT field, just use it
     246          10 :   if (nAvFlds==1)
     247           0 :     vcalmap_.set(ctFlds(0));
     248             :   else {
     249             :     // For each MS field, find the nearest available CT field
     250          10 :     MDirection msdir,ctdir;
     251          10 :     Vector<Double> sep(nAvFlds);
     252          10 :     IPosition ipos(1,0);  // get the first direction stored (no poly yet)
     253          40 :     for (Int iMSFld=0;iMSFld<nMSFlds;++iMSFld) {
     254          30 :       msdir=msfc.phaseDirMeas(iMSFld); // MS fld dir
     255          30 :       sep.set(DBL_MAX);
     256          90 :       for (Int iCTFld=0;iCTFld<nAvFlds;++iCTFld) {
     257             :         // Get cal field direction, converted to ms field frame
     258          60 :         ctdir=ctc.field().phaseDirMeas(ctFlds(iCTFld));
     259          60 :         MDirection::Convert(ctdir,msdir.getRef());
     260          60 :         sep(iCTFld)=ctdir.getValue().separation(msdir.getValue());
     261             :       }
     262             :       // Sort separations
     263          30 :       Vector<uInt> ord;
     264          30 :       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          30 :       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          30 :       vcalmap_(iMSFld)=ctFlds(ord(0));
     276          30 :     }
     277          10 :   }
     278             :   //  cout << "FCM::setNearestFieldMap:***************  vcalmap_ = " << vcalmap_ << endl;
     279             : 
     280          10 : }
     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          12 : void FieldCalMap::setSelectedFieldMap(const String& fieldsel,
     308             :                                       const MeasurementSet& ms,
     309             :                                       const NewCalTable& ct,
     310             :                                       String& extfldsel) {
     311             : 
     312             :   
     313             :   try {
     314             : 
     315          12 :     CTInterface cti(ct);
     316          12 :     MSSelection mss;
     317          12 :     mss.setFieldExpr(fieldsel);
     318          12 :     TableExprNode ten=mss.toTableExprNode(&cti);
     319          12 :     Vector<Int> fieldlist=mss.getFieldList();
     320             : 
     321             :     // Nominally, don't apply selection externally
     322          12 :     extfldsel="";
     323             : 
     324          12 :     if (fieldlist.nelements()>1)
     325             :       // trigger external selection (multi-field case; see old gainfield)
     326           1 :       extfldsel=fieldsel;
     327          11 :     else if (fieldlist.nelements()==1) {
     328             :       // exactly 1 means fldmap=[f]*nfld
     329             : 
     330             :       // Nominally, this many fields need a map
     331          11 :       Int nMSFlds=ms.field().nrow();
     332          11 :       vcalmap_.resize(nMSFlds);
     333          11 :       vcalmap_.set(fieldlist[0]);
     334             :     }
     335             :     else
     336             :       // Field selection found no indices
     337           0 :       throw(AipsError(fieldsel+" matches no fields."));
     338             : 
     339          12 :   }
     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          12 :   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           1 : ScanCalMap::ScanCalMap(const String scancalmap, const MeasurementSet& ms) :
     375           1 :   CalMap()
     376             : {
     377           1 :   if (scancalmap=="self") {
     378           1 :     MSColumns msc(ms);
     379           1 :     Vector<Int> msScans;
     380           1 :     msc.scanNumber().getColumn(msScans);
     381           1 :     vcalmap_.resize(max(msScans) + 1);
     382           1 :     indgen(vcalmap_);
     383           1 :   } else {
     384           0 :     throw(AipsError("Scan mapping failure: unrecognized keyword '" + 
     385           0 :                     scancalmap + "'."));
     386             :   }
     387           1 : }
     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          53 : CalLibSlice::CalLibSlice(const Record& clslice,
     403             :                          const MeasurementSet& ms,
     404          53 :                          const NewCalTable& ct) :
     405          53 :   obs(),scan(),fld(),ent(),spw(),
     406          53 :   tinterp(),finterp(),
     407          53 :   obsmap(),scanmap(),fldmap(),spwmap(),antmap(),
     408          53 :   extfldsel("")
     409             : {
     410             : 
     411          53 :   if (clslice.isDefined("obs")) {
     412          53 :     obs=clslice.asString("obs");
     413             :   }
     414          53 :   if (clslice.isDefined("scan")) {
     415          14 :     scan=clslice.asString("scan");
     416             :   }
     417          53 :   if (clslice.isDefined("field")) {
     418          53 :     fld=clslice.asString("field");
     419             :   }
     420          53 :   if (clslice.isDefined("intent")) {
     421          53 :     ent=clslice.asString("intent");
     422             :   }
     423          53 :   if (clslice.isDefined("spw")) {
     424          53 :     spw=clslice.asString("spw");
     425             :   }
     426             : 
     427          53 :   if (clslice.isDefined("tinterp")) {
     428          53 :     tinterp=clslice.asString("tinterp");
     429          53 :     if (tinterp=="") {
     430           2 :       tinterp="linear";
     431             :     }
     432             :   }
     433          53 :   if (clslice.isDefined("finterp")) {
     434          53 :     finterp=clslice.asString("finterp");
     435             :   }
     436             : 
     437          53 :   if (clslice.isDefined("obsmap")) {
     438             :     //cout << "obsmap.dataType() = " << clslice.dataType("obsmap") << endl;
     439          53 :     if (clslice.dataType("obsmap")==TpArrayInt)
     440          14 :       obsmap=CalMap(Vector<Int>(clslice.asArrayInt("obsmap")));
     441          53 :     if (clslice.dataType("obsmap")==TpString)
     442           0 :       obsmap=ObsCalMap(clslice.asString("obsmap"),ms);
     443             :   }
     444          53 :   if (clslice.isDefined("scanmap")) {
     445             :     //cout << "scanmap.dataType() = " << clslice.dataType("scanmap") << endl;
     446          14 :     if (clslice.dataType("scanmap")==TpArrayInt)
     447          13 :       scanmap=CalMap(Vector<Int>(clslice.asArrayInt("scanmap")));
     448          14 :     if (clslice.dataType("scanmap")==TpString)
     449           1 :       scanmap=ScanCalMap(clslice.asString("scanmap"),ms);
     450             :   }
     451          53 :   if (clslice.isDefined("fldmap")) {
     452          53 :     if (clslice.dataType("fldmap")==TpArrayInt)
     453          19 :       fldmap=FieldCalMap(clslice.asArrayInt("fldmap"));
     454          53 :     if (clslice.dataType("fldmap")==TpString)
     455          22 :       fldmap=FieldCalMap(clslice.asString("fldmap"),ms,ct,extfldsel);
     456             :   }
     457          53 :   if (clslice.isDefined("spwmap")) {
     458             :     //    cout << "spwmap.dataType() = " << clslice.dataType("spwmap") << endl;
     459          53 :     if (clslice.dataType("spwmap")==TpArrayInt)
     460          14 :       spwmap=CalMap(clslice.asArrayInt("spwmap"));
     461             :   }
     462          53 :   if (clslice.isDefined("antmap")) {
     463             :     //    cout << "antmap.dataType() = " << clslice.dataType("antmap") << endl;
     464          53 :     if (clslice.dataType("antmap")==TpArrayInt)
     465          14 :       antmap=CalMap(clslice.asArrayInt("antmap"));
     466             :   }
     467             : 
     468             : 
     469          53 : }
     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          36 : Bool CalLibSlice::validateCLS(const Record& clslice) {
     493             : 
     494          36 :   String missing("");
     495             : 
     496          36 :   if (!clslice.isDefined("obs"))
     497           0 :     missing+="obs ";
     498          36 :   if (!clslice.isDefined("scan"))
     499          36 :     missing+="scan ";
     500          36 :   if (!clslice.isDefined("field"))
     501           0 :     missing+="field ";
     502          36 :   if (!clslice.isDefined("intent"))
     503           0 :     missing+="intent ";
     504          36 :   if (!clslice.isDefined("spw"))
     505           0 :     missing+="spw ";
     506             : 
     507          36 :   if (!clslice.isDefined("tinterp"))
     508           0 :     missing+="tinterp ";
     509          36 :   if (!clslice.isDefined("finterp"))
     510           0 :     missing+="finterp ";
     511             : 
     512          36 :   if (!clslice.isDefined("obsmap"))
     513           0 :     missing+="obsmap ";
     514          36 :   if (!clslice.isDefined("scanmap"))
     515          36 :     missing+="scanmap ";
     516          36 :   if (!clslice.isDefined("fldmap"))
     517           0 :     missing+="fldmap ";
     518          36 :   if (!clslice.isDefined("spwmap"))
     519           0 :     missing+="swmap ";
     520          36 :   if (!clslice.isDefined("antmap"))
     521           0 :     missing+="antmap";
     522             : 
     523             : 
     524          36 :   if (missing.length()>0)
     525          36 :     throw(AipsError(missing));
     526             : 
     527             :   // Everything is ok if we get here
     528           0 :   return true;
     529          36 : }
     530             : 
     531             : 
     532          53 : String CalLibSlice::state() {
     533             : 
     534          53 :   ostringstream o;
     535             : 
     536         106 :   o << "     MS: obs="+obs+" scan="+scan+" fld="+fld+" intent="+ent+" spw="+spw << endl
     537         106 :     << "     CT: tinterp="+tinterp << " finterp="+finterp << endl
     538         106 :     << "         obsmap=" << obsmap.vmap()
     539         106 :     << "         scanmap=" << scanmap.vmap()
     540          53 :     << "         fldmap=";
     541             : 
     542          53 :   if (extfldsel=="")
     543          52 :     o << fldmap.vmap();
     544             :   else
     545           1 :     o << extfldsel+" (by selection)";
     546             : 
     547          53 :   o << endl
     548          53 :     << "         spwmap=" << spwmap.vmap()
     549         106 :     << "         antmap=" << antmap.vmap()
     550          53 :     << endl;
     551             : 
     552         106 :   return String(o);
     553          53 : }
     554             : 
     555             : 
     556         349 : CLPPResult::CLPPResult() :
     557         349 :   result_(),resultFlag_()
     558         349 : {}
     559             : 
     560           0 : CLPPResult::CLPPResult(const IPosition& shape) :
     561           0 :   result_(shape,0.0),
     562           0 :   resultFlag_(shape,true)
     563           0 : {}
     564          76 : CLPPResult::CLPPResult(uInt nPar,uInt nFPar,uInt nChan,uInt nElem) :
     565          76 :   result_(nFPar,nChan,nElem,0.0),
     566          76 :   resultFlag_(nPar,nChan,nElem,true)
     567          76 : {}
     568             : 
     569         258 : CLPPResult& CLPPResult::operator=(const CLPPResult& other) {
     570             :   // Avoid Array deep copies!
     571         258 :   result_.reference(other.result_);
     572         258 :   resultFlag_.reference(other.resultFlag_);
     573         258 :   return *this;
     574             : }
     575          10 : void CLPPResult::resize(uInt nPar,uInt nFPar,uInt nChan,uInt nElem) {
     576          10 :   result_.resize(nFPar,nChan,nElem);
     577          10 :   result_.set(1.0);
     578          10 :   resultFlag_.resize(nPar,nChan,nElem);
     579          10 :   resultFlag_.set(true);
     580          10 : }
     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          43 : CLPatchPanel::CLPatchPanel(const String& ctname,
     878             :                            const MeasurementSet& ms,
     879             :                            const Record& callib,
     880             :                            VisCalEnum::MatrixType mtype,
     881             :                            Int nPar,
     882          43 :                            const CTTIFactoryPtr cttifactoryptr) :
     883          43 :   ct_(NewCalTable::createCT(ctname,Table::Old,Table::Memory)),
     884          43 :   ctasms_(), // null, in this context
     885          43 :   ms_(ms),
     886          43 :   mtype_(mtype),
     887          43 :   isCmplx_(false),
     888          43 :   nPar_(nPar),
     889          43 :   nFPar_(nPar),
     890             : 
     891          43 :   nChanIn_(),
     892          43 :   freqIn_(),
     893          43 :   refFreqIn_(),
     894             : 
     895          43 :   nMSObs_(ms_.observation().nrow()),
     896          43 :   nMSFld_(ms_.field().nrow()),
     897          43 :   nMSSpw_(ms_.spectralWindow().nrow()),
     898          43 :   nMSAnt_(ms_.antenna().nrow()),
     899          43 :   nMSElem_(0),  // set non-trivially in ctor body
     900          43 :   nCTObs_(ct_.observation().nrow()),
     901          43 :   nCTFld_(ct_.field().nrow()),
     902          43 :   nCTSpw_(ct_.spectralWindow().nrow()),
     903          43 :   nCTAnt_(ct_.antenna().nrow()),
     904          43 :   nCTElem_(0), // set non-trivially in ctor body
     905          43 :   lastresadd_(nMSSpw_,NULL),
     906         129 :   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          43 :   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          43 :   case VisCalEnum::JONES: {
     933          43 :     nCTElem_=nCTAnt_;
     934          43 :     nMSElem_=nMSAnt_;
     935          43 :     break;
     936             :   }
     937             :   }
     938             : 
     939             :   // How many _Float_ parameters?
     940          43 :   isCmplx_=(ct_.keywordSet().asString("ParType")=="Complex");
     941          43 :   if (isCmplx_) // Complex input
     942          41 :     nFPar_*=2;  // interpolating 2X as many Float values
     943             : 
     944             :   // Set CT channel/freq info
     945          86 :   CTSpWindowColumns ctspw(ct_.spectralWindow());
     946          43 :   ctspw.numChan().getColumn(nChanIn_);
     947          43 :   freqIn_.resize(nCTSpw_);
     948         101 :   for (uInt iCTspw=0;iCTspw<ctspw.nrow();++iCTspw)
     949          58 :     ctspw.chanFreq().get(iCTspw,freqIn_(iCTspw),true);
     950          43 :   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          43 :   uInt nkey=callib.nfields();
     961          43 :   Int icls=-1;
     962         182 :   for (uInt ikey=0;ikey<nkey;++ikey) {
     963             : 
     964             :     // CalTable name might be needed below (eg, for log messages)
     965         139 :     String ctname=Path(ct_.getPartNames()[0]).baseName().before(".tempMemCalTable");
     966             : 
     967             :     // Ignore non-Record members in the callib
     968         139 :     if (callib.type(ikey)!=TpRecord)
     969          86 :       continue;
     970             : 
     971          53 :     ++icls;
     972             : 
     973             :     // The current callib slice
     974         106 :     CalLibSlice cls(callib.asRecord(ikey),ms_,ct_);
     975             : 
     976             : 
     977             :     // Reference to the cal table that we'll use below
     978          53 :     NewCalTable ct0(ct_);
     979          53 :     if (cls.extfldsel!="") {
     980             :       // Select on the reference table
     981             :       try {
     982           1 :         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          53 :     MeasurementSet clsms(ms_);
     992             : 
     993             :     // Trap Null selection exceptions, as they are not needed downstream
     994             :     try {
     995          60 :       this->selectOnMS(clsms,ms_,cls.obs,cls.scan,cls.fld,cls.ent,cls.spw,"");
     996             :     }
     997           7 :     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          14 :       String msname=Path(ms_.getPartNames()[0]).baseName();
    1003           7 :       logsink_ << LogIO::WARN
    1004             :                << ".   The following callib entry matches no data" << endl
    1005          14 :                << ".   in the selected MS ("+msname+") and will be ignored:" << endl
    1006             :                << ".   " << icls << ":" << endl
    1007          14 :                << cls.state() << LogIO::POST;
    1008             : 
    1009             :       // Just jump to next slice (cal maps for this MS selection unneeded)
    1010           7 :       continue;
    1011           7 :     }
    1012             : 
    1013             :     // Report this relevant cal lib slice to the logger
    1014          46 :     logsink_ << LogIO::NORMAL << ".   " << icls << ":" << endl
    1015          46 :              << 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          46 :     Vector<Int> reqMSobs(1,-1); // assume all, indescriminately
    1028          46 :     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          46 :     Vector<Int> reqMSscan(1,-1); // assume all, indescriminately
    1035          46 :     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          46 :     Vector<Int> reqMSfld(1,-1); // assume all, indescriminately
    1042          46 :     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          13 :       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          46 :     Vector<Int> reqMSint(1,-1); // assume all
    1049          46 :     if (cls.ent.length()>0)
    1050           0 :       reqMSint.reference(this->getMSuniqueIds(clsms,"intent"));
    1051             :     //cout << "reqMSint = " << reqMSint << endl;
    1052          46 :     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          46 :     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          46 :     Vector<Int> reqMSant(nMSAnt_);
    1063          46 :     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          46 :     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          46 :     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          92 :     for (uInt iCTobs=0;iCTobs<reqCTobs.nelements();++iCTobs) {
    1079          46 :       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          46 :       Vector<Int> theseMSobs=cls.obsmap.msids(thisCTobs,reqMSobs);
    1084          46 :       if (theseMSobs.nelements()==1 && theseMSobs[0]<0)
    1085          46 :         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          46 :         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          46 :       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          46 :       NewCalTable scanselCT(obsselCT);
    1108          95 :       for (uInt iCTscan=0;iCTscan<reqCTscan.nelements();++iCTscan) {
    1109          49 :         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          49 :         Vector<Int> theseMSscan=cls.scanmap.msids(thisCTscan,reqMSscan);
    1114          49 :         if (theseMSscan.nelements()==1 && theseMSscan[0]<0)
    1115          45 :           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          49 :           if (thisCTscan!=-1) 
    1122           9 :             this->selectOnCT(scanselCT,obsselCT,"",String::toString(thisCTscan),"","","");
    1123             :         }
    1124           1 :         catch (...) {  //  MSSelectionNullSelection x ) {
    1125             : 
    1126             :           // Required CT scan does not exist in the caltable
    1127           1 :           recordBadMSIndices(reqMSobs,theseMSscan,reqMSfld,reqMSint,Vector<Int>(1,-1));  // all spws
    1128           1 :           continue;  // jump to next CT scan
    1129           1 :         }
    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          48 :         Vector<Int> reqCTfld=cls.fldmap.ctids(reqMSfld);
    1135             :         //cout << " reqCTfld = " << reqCTfld << endl;
    1136             : 
    1137             :         // For each required CT fld:
    1138          48 :         NewCalTable fldselCT(scanselCT);
    1139         107 :         for (uInt iCTfld=0;iCTfld<reqCTfld.nelements();++iCTfld) {
    1140          59 :           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          59 :           Vector<int> theseMSfld=cls.fldmap.msids(thisCTfld,reqMSfld);
    1145          59 :           if (theseMSfld.nelements()==1 && theseMSfld[0]<0)
    1146          22 :             theseMSfld.reference(reqMSfld);
    1147             : 
    1148             :           //cout << "  thisCTfld = " << thisCTfld << ": theseMSfld = " << theseMSfld << endl;
    1149             : 
    1150             :           // Apply thisCTfld selection to the CT
    1151             :           try {
    1152          59 :             if (thisCTfld!=-1)
    1153          37 :               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          59 :           NewCalTable spwselCT(fldselCT);
    1166         135 :           for (uInt iMSspw=0;iMSspw<reqMSspw.nelements();++iMSspw) {
    1167          76 :             Int& thisMSspw=reqMSspw(iMSspw);
    1168          76 :             Int thisCTspw=cls.spwmap(thisMSspw);
    1169          76 :             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          76 :               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          76 :             CTCalPatchKey iclTres(icls,thisCTobs,thisCTscan,thisCTfld,thisMSspw,-1);
    1194          76 :             clTres_[iclTres]=CLPPResult(nPar_,nFPar_,nChanIn_[thisCTspw],nMSElem_);
    1195             : 
    1196          76 :             NewCalTable antselCT(spwselCT);
    1197          76 :     Bool doLinkResults(true);  // initialize true, so it will happen if relevant code reached
    1198        1618 :             for (uInt iMSant=0;iMSant<reqMSant.nelements();++iMSant) {
    1199        1542 :               Int& thisMSant=reqMSant(iMSant);
    1200        1542 :               Int thisCTant=cls.antmap(thisMSant);
    1201        1542 :               if (thisCTant<0) thisCTant=thisMSant;
    1202             : 
    1203             :               // Apply thisCTant selection to CT
    1204             :               try {
    1205        1562 :                 this->selectOnCT(antselCT,spwselCT,"","","","",String::toString(thisCTant));
    1206             :               }
    1207           4 :               catch ( MSSelectionNullSelection x ) {
    1208             :                 // Log a warning about the missing antenna
    1209           4 :                 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           4 :                          << LogIO::POST;
    1212             :                 // Step to next antenna
    1213           4 :                 continue;
    1214           4 :               }
    1215             : 
    1216             :               // Make the Cal Interpolator (icls is the CL slice index):
    1217        1538 :               CTCalPatchKey ici0(icls,thisCTobs,thisCTscan,thisCTfld,thisCTspw,thisCTant);  // all CT indices
    1218        1538 :               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        1538 :               if (ci_.count(ici1)<1) {
    1222        1538 :                 ciname_[ici1]=ici0.print()+" rows="+String::toString(antselCT.nrow());
    1223        1538 :                 Array<Float> r(clTres_[iclTres].result(thisMSant));
    1224        1538 :                 Array<Bool> rf(clTres_[iclTres].resultFlag(thisMSant));
    1225        1538 :                 ci_[ici1]=(*cttifactoryptr_)(antselCT,cls.tinterp,r,rf);
    1226             :                 //if (iMSant==0) cout << "   Creating: CT("<<ici1.print() << ") --> CT(" << ici0.print() << ")  (all antennas)" << endl;
    1227        1538 :               }
    1228             :               else
    1229           0 :                 throw(AipsError("Attempted duplicate CTCalPatchKey!"));
    1230             : 
    1231             :               // Now distribute this CTTimeInterp1 instance to all relevant MS indices
    1232        3076 :               for (uInt iMSobs=0;iMSobs<theseMSobs.nelements();++iMSobs) {
    1233        1538 :                 Int& thisMSobs=theseMSobs(iMSobs);
    1234        3076 :                 for (uInt iMSscan=0;iMSscan<theseMSscan.nelements();++iMSscan) {
    1235        1538 :                   Int& thisMSscan=theseMSscan(iMSscan);
    1236        3460 :                   for (uInt iMSfld=0;iMSfld<theseMSfld.nelements();++iMSfld) {
    1237        1922 :                     Int& thisMSfld=theseMSfld(iMSfld);
    1238        3844 :                     for (uInt iMSint=0;iMSint<theseMSint.nelements();++iMSint) {
    1239        1922 :                       Int& thisMSint=theseMSint(iMSint);
    1240             : 
    1241        1922 :                       MSCalPatchKey ims(thisMSobs,thisMSscan,thisMSfld,thisMSint,thisMSspw,thisMSant);
    1242        1922 :                       if (msci_.count(ims)<1) {
    1243        1922 :                         msciname_[ims]=ciname_[ici1];
    1244        1922 :                         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        1922 :       if (doLinkResults) {
    1255          91 :                         MSCalPatchKey imsgroup(thisMSobs,thisMSscan,thisMSfld,thisMSint,thisMSspw,-1);
    1256          91 :                         msTres_[imsgroup]=clTres_[iclTres];
    1257          91 :                         msFres_[imsgroup]=CLPPResult(); // this will be resized on-demand
    1258          91 :                         ctspw_[imsgroup]=thisCTspw;
    1259          91 :                         finterp_[imsgroup]=cls.finterp;
    1260          91 :                       }
    1261        1922 :                     } // iMSint
    1262             :                   } // iMSfld
    1263             :                 } // iMSscan
    1264             :               } // iMSobs
    1265        1538 :         doLinkResults = False; // Don't do it again
    1266        1538 :             } // iMSant
    1267          76 :           } // iMSspw
    1268          59 :         } // iCTfld
    1269          49 :       } // iCTscan
    1270          46 :     } // iCTobs
    1271             : 
    1272             : 
    1273         160 :   } // icls
    1274             : 
    1275             : 
    1276             : 
    1277          43 : } // ctor
    1278             : 
    1279           1 : 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           2 :   for (uInt iobs=0;iobs<obs.nelements();++iobs) {
    1286           2 :     for (uInt iscan=0;iscan<scan.nelements();++iscan) {
    1287           2 :       for (uInt ifld=0;ifld<fld.nelements();++ifld) {
    1288           2 :         for (uInt ient=0;ient<ent.nelements();++ient) {
    1289           2 :           for (uInt ispw=0;ispw<spw.nelements();++ispw) {
    1290             : 
    1291           1 :             MSCalPatchKey ims(obs[iobs],scan[iscan],fld[ifld],ent[ient],spw[ispw],-1);  // All ants
    1292           1 :             if (badmsciname_.count(ims)<1) {
    1293           1 :               badmsciname_[ims]=ims.print();
    1294             :               //cout << "   Bad MS indices: " << ims.print() << endl;
    1295             :             }
    1296           1 :           }
    1297             :         }
    1298             :       }
    1299             :     }
    1300             :   }
    1301           1 :   return;
    1302             : }
    1303             : 
    1304             : 
    1305        1660 : 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        1660 :   String taql("");
    1311        1660 :   if (ant1.length()>0)
    1312        1542 :     taql="ANTENNA1=="+ant1;
    1313             : 
    1314             :   //  cout << "taql = >>" << taql << "<<" << endl;
    1315             : 
    1316             :   // Forward to generic method (sans intent)
    1317        1660 :   CTInterface cti(ctin);
    1318        1670 :   this->selectOnCTorMS(ctout,cti,obs,scan,fld,"",spw,"",taql);
    1319             : 
    1320        1665 : }
    1321             : 
    1322          53 : 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          53 :   MSInterface msi(msin);
    1329          60 :   this->selectOnCTorMS(msout,msi,obs,scan,fld,ent,spw,ant,"");
    1330             : 
    1331          53 : }
    1332        1713 : 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        1713 :   MSSelection mss;
    1339        1713 :   if (obs.length()>0)
    1340           0 :     mss.setObservationExpr(obs);
    1341        1713 :   if (scan.length()>0)
    1342           4 :     mss.setScanExpr(scan);
    1343        1713 :   if (fld.length()>0)
    1344          58 :     mss.setFieldExpr(fld);
    1345        1713 :   if (ent.length()>0)
    1346           0 :     mss.setStateExpr(ent);
    1347        1713 :   if (spw.length()>0)
    1348          76 :     mss.setSpwExpr(spw);
    1349        1713 :   if (ant.length()>0)
    1350           0 :     mss.setAntennaExpr(ant); // this will not behave as required for Jones caltables (its bl-based selection!)
    1351        1713 :   if (taql.length()>0)
    1352        1542 :     mss.setTaQLExpr(taql);
    1353             : 
    1354             : 
    1355        1713 :   TableExprNode ten=mss.toTableExprNode(&msst);
    1356        1725 :   getSelectedTable(ctout,*(msst.table()),ten,"");
    1357             : 
    1358        1725 : }
    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          59 : Vector<Int> CLPatchPanel::getMSuniqueIds(MeasurementSet& ms, String which) {
    1386             : 
    1387          59 :   MSColumns msc(ms);
    1388             : 
    1389             :   // Extract the requested column as a Vector
    1390          59 :   Vector<Int> colv;
    1391          59 :   if (which=="obs")
    1392           0 :     msc.observationId().getColumn(colv);
    1393          59 :   else if (which=="scan")
    1394           0 :     msc.scanNumber().getColumn(colv);
    1395          59 :   else if (which=="fld")
    1396          13 :     msc.fieldId().getColumn(colv);
    1397          46 :   else if (which=="intent")
    1398           0 :     msc.stateId().getColumn(colv);
    1399          46 :   else if (which=="spw")
    1400          46 :     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          59 :   uInt nuniq=genSort(colv,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
    1406          59 :   colv.resize(nuniq,true);
    1407             : 
    1408             :   // In case of spw, translate from ddid
    1409          59 :   if (which=="spw") {
    1410          46 :     Vector<Int> spwids;
    1411          46 :     msc.dataDescription().spectralWindowId().getColumn(spwids);
    1412         107 :     for (uInt i=0;i<colv.nelements();++i)
    1413          61 :       colv[i]=spwids[colv[i]];
    1414          46 :   }
    1415             : 
    1416             :   // return the value
    1417         118 :   return colv;
    1418             : 
    1419          59 : }
    1420             : 
    1421             : 
    1422             : 
    1423             : // Destructor
    1424          86 : CLPatchPanel::~CLPatchPanel() {
    1425             : 
    1426             :   if (CTPATCHPANELVERB) cout << "CLPatchPanel::~CLPatchPanel()" << endl;
    1427             : 
    1428             :   // Delete the atomic interpolators
    1429        1581 :   for (std::map<CTCalPatchKey,CTTimeInterp1*>::iterator it=ci_.begin(); it!=ci_.end(); ++it)
    1430        1538 :     delete it->second;
    1431             : 
    1432          86 : }
    1433             : 
    1434             : // Is specific calibration explicitly available for a obs,scan,fld,intent,spw,ant combination?
    1435       30582 : Bool CLPatchPanel::calAvailable(casacore::Int msobs, casacore::Int msscan, casacore::Int msfld,
    1436             :                                 casacore::Int msent, casacore::Int msspw, casacore::Int msant) {
    1437             : 
    1438       30582 :   const MSCalPatchKey key(msobs,msscan,msfld,msent,msspw,msant);
    1439             : 
    1440       30582 :   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       30582 :   return avail;
    1447             : 
    1448       30582 : }
    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       30042 : Bool CLPatchPanel::MSIndicesOK(casacore::Int msobs, casacore::Int msscan, casacore::Int msfld,
    1462             :                                casacore::Int msent, casacore::Int msspw, casacore::Int msant) {
    1463             : 
    1464       30042 :   const MSCalPatchKey key(msobs,msscan,msfld,msent,msspw,msant);
    1465       30042 :   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       30042 :   return !bad;
    1472             : 
    1473       30042 : }
    1474             : 
    1475       27942 : 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       27942 :   Bool newcal(false);
    1480             : 
    1481             :   // resultC and resFlag will be unchanged if newcal remains false
    1482       27942 :   Cube<Float> f; // temporary to reference Float interpolation result
    1483       27942 :   newcal=interpolate(f,resFlag,msobs,msscan,msfld,msent,msspw,time,freq);
    1484       27942 :   if (newcal)
    1485             :     // convert to complex and have resultC take over ownership
    1486       25870 :     resultC.reference(RIorAPArray(f).c());
    1487             : 
    1488       27942 :   return newcal;  // If T, calling scope should act nontrivially, if necessary
    1489             : 
    1490       27942 : }
    1491             : 
    1492         240 : 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         240 :   Bool newcal(false);
    1497             : 
    1498             :   // resultC and resFlag will be unchanged if newcal remains false
    1499         240 :   Cube<Float> f; // temporary to reference Float interpolation result
    1500         240 :   newcal=interpolate(f,resFlag,msobs,msscan,msfld,msent,msspw,time,freq);
    1501             : 
    1502         240 :   if (newcal)
    1503             :     // convert to complex and have resultC take over ownership
    1504          10 :     resultC.reference(RIorAPArray(f).c());
    1505             : 
    1506         240 :   return newcal;  // If T, calling scope should act nontrivially, if necessary
    1507             : 
    1508         240 : }
    1509             : 
    1510             : 
    1511       30342 : 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       30342 :   Bool newcal(false);
    1519             : 
    1520             :   // Suppled arrays reference the result (if available)
    1521       30342 :   MSCalPatchKey ires(msobs,msscan,msfld,msent,msspw,-1);
    1522             : 
    1523             :   // Trap lack of available calibration for requested obs,fld,intent,spw
    1524       30342 :   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       30342 :   if (lastresadd_(msspw)!=msTres_[ires].result_.data())
    1533          78 :     newcal=true;
    1534       30342 :   lastresadd_(msspw)=msTres_[ires].result_.data();
    1535             : 
    1536             :   // Loop over _output_ elements
    1537      344802 :   for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
    1538             :     // Call fully _patched_ time-interpolator, keeping track of 'newness'
    1539             :     //  fills ctTres_ implicitly
    1540      314460 :     MSCalPatchKey ims(msobs,msscan,msfld,msent,msspw,iMSElem);
    1541      314460 :     if (msci_.count(ims)>0) {
    1542      307152 :       if (freq>0.0)
    1543      280752 :         newcal|=msci_[ims]->interpolate(time,freq);
    1544             :       else
    1545       26400 :         newcal|=msci_[ims]->interpolate(time);
    1546             :     }
    1547      314460 :   }
    1548             : 
    1549       30342 :   if (newcal) {
    1550       27078 :     resultR.reference(msTres_[ires].result_);
    1551       27078 :     resFlag.reference(msTres_[ires].resultFlag_);
    1552             :   }
    1553             : 
    1554       30342 :   return newcal; // If true, calling scope should act
    1555       30342 : }
    1556             : 
    1557             : 
    1558         240 : 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         240 :   Bool newcal(false);
    1566             : 
    1567             :   // Suppled arrays reference the result (if available)
    1568         240 :   MSCalPatchKey ires(msobs,msscan,msfld,msent,msspw,-1);
    1569             : 
    1570             :   // Trap lack of available calibration for requested obs,scan,fld,intent,spw
    1571         240 :   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         240 :   if (lastresadd_(msspw)!=msTres_[ires].result_.data())
    1580          10 :     newcal=true;
    1581         240 :   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         240 :   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         240 :   uInt nMSChan=freq.nelements();    // The number of requested channels
    1588         240 :   if (msFres_.count(ires)>0) {
    1589         470 :     if (msFres_[ires].result_.nelements()==0 ||
    1590         470 :         msFres_[ires].result(0).ncolumn()!=nMSChan) {
    1591          10 :       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        6960 :   for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
    1599             :     // Call time interpolation calculation; resample in freq if new
    1600             :     //   (fills msTRes_ implicitly)
    1601             : 
    1602        6720 :     MSCalPatchKey ims(msobs,msscan,msfld,msent,msspw,iMSElem);
    1603        6720 :     if (msci_.count(ims)>0) {
    1604        6720 :       if (msci_[ims]->interpolate(time) || forceFinterp) {
    1605             : 
    1606             :         // Resample in frequency
    1607         280 :         Matrix<Float>   fR( msFres_[ires].result(iMSElem) );
    1608         280 :         Matrix<Bool> fRflg( msFres_[ires].resultFlag(iMSElem) );
    1609         280 :         Matrix<Float>   tR( msTres_[ires].result(iMSElem) );
    1610         280 :         Matrix<Bool> tRflg( msTres_[ires].resultFlag(iMSElem) );
    1611         280 :         resampleInFreq(fR,fRflg,freq,tR,tRflg,freqIn_(ctspw_[ires]),finterp_[ires]);
    1612             : 
    1613             :         // Calibration is new
    1614         280 :         newcal=true;
    1615         280 :       }
    1616             :     }
    1617        6720 :   }
    1618             : 
    1619         240 :   if (newcal) {
    1620             :     // Supplied arrays to reference the result
    1621          10 :     resultR.reference(msFres_[ires].result_);
    1622          10 :     resFlag.reference(msFres_[ires].resultFlag_);
    1623             :   }
    1624             : 
    1625         240 :   return newcal;
    1626         240 : }
    1627             : 
    1628          10 : Bool CLPatchPanel::getTresult(Cube<Float>& resultR, Cube<Bool>& resFlag,
    1629             :                               Int obs, Int scan, Int fld, Int ent, Int spw) {
    1630             : 
    1631          10 :   MSCalPatchKey mskey(obs,scan,fld,ent,spw,-1);
    1632             : 
    1633          10 :   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          10 :   resultR.reference(msTres_[mskey].result_);
    1640          10 :   resFlag.reference(msTres_[mskey].resultFlag_);
    1641             : 
    1642          10 :   return true;
    1643             : 
    1644          10 : }
    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         280 : 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         280 :   if (nfalse(tflg)==0) {
    1703          10 :     fflg.set(true);
    1704          10 :     return;
    1705             :   }
    1706             : 
    1707         270 :   Int flparmod=nFPar_/nPar_;    // for indexing the flag Matrices on the par axis
    1708             : 
    1709         270 :   Bool unWrapPhase=flparmod>1;
    1710             : 
    1711             :   //  cout << "nFPar_,nPar_,flparmod = " << nFPar_ << "," << nPar_ << "," << flparmod << endl;
    1712             : 
    1713         270 :   fres=0.0;
    1714             : 
    1715        1350 :   for (Int ifpar=0;ifpar<nFPar_;++ifpar) {
    1716             : 
    1717             :     // Slice by par (each Matrix row)
    1718        1080 :     Vector<Float> fresi(fres.row(ifpar)), tresi(tres.row(ifpar));
    1719        1080 :     Vector<Bool> fflgi(fflg.row(ifpar/flparmod)), tflgi(tflg.row(ifpar/flparmod));
    1720             : 
    1721             :     // Mask time result by flags
    1722        2160 :     Vector<Double> mfin=fin(!tflgi).getCompressedArray();
    1723             : 
    1724        1080 :     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        1080 :     mfin/=1.0e9; // in GHz
    1734        2160 :     Vector<Float> mtresi=tresi(!tflgi).getCompressedArray();
    1735             : 
    1736             :     // Trap case of same in/out frequencies
    1737        1080 :     if (fout.nelements()==mfin.nelements() && allNear(fout,mfin,1.e-10)) {
    1738             :       // Just copy
    1739        1080 :       fresi=mtresi;
    1740        1080 :       fflgi.set(false);  // none are flagged
    1741        1080 :       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        6480 :   }
    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

Generated by: LCOV version 1.16