Commits

George Moellenbrock authored 8815a74e3b3
For CAS-13266:  Moved logging of missing calibration from caltable read / interpolation prep stages to interpolation calculation stage to avoid generating spurious log messages (e.g., cal missing for unselected data)
No tags

casatools/src/code/synthesis/CalTables/CTPatchedInterp.cc

Modified
53 53 CTPatchedInterp::CTPatchedInterp(NewCalTable& ct,
54 54 VisCalEnum::MatrixType mtype,
55 55 Int nPar,
56 56 const String& timetype,
57 57 const String& freqtype,
58 58 const String& fieldtype,
59 59 Vector<Int> spwmap,
60 60 Vector<Int> fldmap,
61 61 const CTTIFactoryPtr cttifactoryptr) :
62 62 ct_(ct),
63 + ctname_(),
63 64 msmc_(NULL),
64 65 mtype_(mtype),
65 66 isCmplx_(false),
66 67 nPar_(nPar),
67 68 nFPar_(nPar),
68 69 timeType_(timetype),
69 70 freqTypeStr_(freqtype),
70 71 relativeFreq_(freqtype.contains("rel")),
71 72 freqInterpMethod0_(ftype(freqTypeStr_)),
72 73 freqInterpMethod_(freqInterpMethod0_),
89 90 spwInOK_(),
90 91 fldMap_(),
91 92 spwMap_(),
92 93 antMap_(),
93 94 elemMap_(),
94 95 conjTab_(),
95 96 result_(),
96 97 resFlag_(),
97 98 tI_(),
98 99 tIdel_(),
100 + tIMissingLogged_(),
99 101 lastFld_(ct.spectralWindow().nrow(),-1),
100 102 lastObs_(ct.spectralWindow().nrow(),-1),
101 103 lastScan_(ct.spectralWindow().nrow(),-1),
102 104 cttifactoryptr_(cttifactoryptr)
103 105 {
104 106 if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::CTPatchedInterp(<no MS>)" << endl;
105 107
106 108 // cout << "freqInterpMethod_ = " << freqInterpMethod_ << endl;
107 109
108 110 freqInterpMethodVec_.resize(nMSSpw_);
276 278 VisCalEnum::MatrixType mtype,
277 279 Int nPar,
278 280 const String& timetype,
279 281 const String& freqtype,
280 282 const String& fieldtype,
281 283 const MeasurementSet& ms,
282 284 const MSMetaInfoForCal& msmc,
283 285 Vector<Int> spwmap,
284 286 const CTTIFactoryPtr cttifactoryptr) :
285 287 ct_(ct),
288 + ctname_(),
286 289 msmc_(&msmc),
287 290 mtype_(mtype),
288 291 isCmplx_(false),
289 292 nPar_(nPar),
290 293 nFPar_(nPar),
291 294 timeType_(timetype),
292 295 freqTypeStr_(freqtype),
293 296 relativeFreq_(freqtype.contains("rel")),
294 297 freqInterpMethod0_(ftype(freqTypeStr_)),
295 298 freqInterpMethod_(freqInterpMethod0_),
312 315 spwInOK_(),
313 316 fldMap_(),
314 317 spwMap_(),
315 318 antMap_(),
316 319 elemMap_(),
317 320 conjTab_(),
318 321 result_(),
319 322 resFlag_(),
320 323 tI_(),
321 324 tIdel_(),
325 + tIMissingLogged_(),
322 326 lastFld_(ms.spectralWindow().nrow(),-1),
323 327 lastObs_(ms.spectralWindow().nrow(),-1),
324 328 lastScan_(ms.spectralWindow().nrow(),-1),
325 329 cttifactoryptr_(cttifactoryptr)
326 330 {
327 331
328 332 if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::CTPatchedInterp(CT,MS)" << endl;
329 333
330 334 //freqInterpMethod_=ftype(freqTypeStr_);
331 335
450 454 << Path(ct_.tableName()).baseName().before(".tempMemCal")
451 455 << "; ignoring 'perscan' interpolation.";
452 456 log << msg.str() << LogIO::WARN;
453 457 } else {
454 458 std::set<Int> scanNumbers = msmc.msmd().getScanNumbers(0,0);
455 459 nCTTimeSeg_=maxScan+1;
456 460 nMSTimeSeg_=*scanNumbers.rbegin()+1;
457 461
458 462 }
459 463 }
460 -
464 +
461 465 // Initialize caltable slices
462 466 sliceTable();
463 467
464 468 // Set spwmap
465 469 setSpwMap(spwmap);
466 470
467 471 // Set fldmap
468 472 if (byField_)
469 473 setFldMap(ms.field()); // on a trial basis
470 474 else
497 501 CTPatchedInterp::CTPatchedInterp(NewCalTable& ct,
498 502 VisCalEnum::MatrixType mtype,
499 503 Int nPar,
500 504 const String& timetype,
501 505 const String& freqtype,
502 506 const String& fieldtype,
503 507 const MSColumns& mscol,
504 508 Vector<Int> spwmap,
505 509 const CTTIFactoryPtr cttifactoryptr) :
506 510 ct_(ct),
511 + ctname_(),
507 512 msmc_(NULL),
508 513 mtype_(mtype),
509 514 isCmplx_(false),
510 515 nPar_(nPar),
511 516 nFPar_(nPar),
512 517 timeType_(timetype),
513 518 freqTypeStr_(freqtype),
514 519 relativeFreq_(freqtype.contains("rel")),
515 520 freqInterpMethod0_(ftype(freqTypeStr_)),
516 521 freqInterpMethod_(freqInterpMethod0_),
531 536 spwInOK_(),
532 537 fldMap_(),
533 538 spwMap_(),
534 539 antMap_(),
535 540 elemMap_(),
536 541 conjTab_(),
537 542 result_(),
538 543 resFlag_(),
539 544 tI_(),
540 545 tIdel_(),
546 + tIMissingLogged_(),
541 547 lastFld_(mscol.spectralWindow().nrow(),-1),
542 548 lastObs_(mscol.spectralWindow().nrow(),-1),
543 549 lastScan_(mscol.spectralWindow().nrow(),-1),
544 550 cttifactoryptr_(cttifactoryptr)
545 551 {
546 552 if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::CTPatchedInterp(mscol)" << endl;
547 553
548 554 // cout << "freqInterpMethod_ = " << freqInterpMethod_ << endl;
549 555 freqInterpMethodVec_.resize(nMSSpw_);
550 556 freqInterpMethodVec_.set(freqInterpMethod0_);
633 639 CTPatchedInterp::CTPatchedInterp(NewCalTable& ct,
634 640 VisCalEnum::MatrixType mtype,
635 641 Int nPar,
636 642 const String& timetype,
637 643 const String& freqtype,
638 644 const String& fieldtype,
639 645 const MeasurementSet& ms,
640 646 Vector<Int> spwmap,
641 647 const CTTIFactoryPtr cttifactoryptr) :
642 648 ct_(ct),
649 + ctname_(),
643 650 msmc_(NULL),
644 651 mtype_(mtype),
645 652 isCmplx_(false),
646 653 nPar_(nPar),
647 654 nFPar_(nPar),
648 655 timeType_(timetype),
649 656 freqTypeStr_(freqtype),
650 657 relativeFreq_(freqtype.contains("rel")),
651 658 freqInterpMethod0_(ftype(freqTypeStr_)),
652 659 freqInterpMethod_(freqInterpMethod0_),
669 676 spwInOK_(),
670 677 fldMap_(),
671 678 spwMap_(),
672 679 antMap_(),
673 680 elemMap_(),
674 681 conjTab_(),
675 682 result_(),
676 683 resFlag_(),
677 684 tI_(),
678 685 tIdel_(),
686 + tIMissingLogged_(),
679 687 lastFld_(ms.spectralWindow().nrow(),-1),
680 688 lastObs_(ms.spectralWindow().nrow(),-1),
681 689 lastScan_(ms.spectralWindow().nrow(),-1),
682 690 cttifactoryptr_(cttifactoryptr)
683 691 {
684 692
685 693 if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::CTPatchedInterp(CT,MS)" << endl;
686 694
687 695 //freqInterpMethod_=ftype(freqTypeStr_);
688 696
845 853 }
846 854 ctSlices_.resize();
847 855 }
848 856 }
849 857
850 858 Bool CTPatchedInterp::interpolate(Int msobs, Int msscan, Int msfld, Int msspw, Double time, Double freq) {
851 859
852 860 if (CTPATCHEDINTERPVERB) cout << "CTPatchedInterp::interpolate(...)" << endl;
853 861
854 862 Bool newcal(false);
855 - IPosition ip(4,0,msspw,msfld,thisTimeSeg(msobs,msscan));
863 + Int scanOrObsInt=thisTimeSeg(msobs,msscan);
864 + IPosition ip(4,0,msspw,msfld,scanOrObsInt);
856 865
857 866 // Loop over _output_ elements
858 867 for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
859 868 // Call fully _patched_ time-interpolator, keeping track of 'newness'
860 869 // fills timeResult_/timeResFlag_ implicitly
861 870 ip(0)=iMSElem;
862 871
863 872 if (!tI_(ip)) {
864 - //cout << "Flagging: " << ip << endl;
873 + //if (iMSElem==0) cout << "Flagging: " << ip << endl;
874 +
875 + // casa log post
876 + if (!tIMissingLogged_(ip)) { // only if these coords not logged before
877 + LogIO log;
878 + ostringstream msg;
879 + String scanOrObsStr=(byScan_ ? "scan" : "obs");
880 + msg << "MS " << scanOrObsStr << "=" << scanOrObsInt;
881 + if (byField_) msg << " in fld=" << msfld;
882 + msg << ",spw=" << msspw
883 + << ",ant=" << iMSElem
884 + << " is selected for processing, but has no available calibration in " << ctname_
885 + << " as mapped, and will be flagged.";
886 + log << msg.str() << LogIO::WARN;
887 + tIMissingLogged_(ip)=true;
888 + }
865 889 newcal=true;
866 890 }
867 891 else {
868 892 if (freq>0.0)
869 893 newcal|=tI_(ip)->interpolate(time,freq);
870 894 else
871 895 newcal|=tI_(ip)->interpolate(time);
872 896 }
873 897 }
874 898
875 899 // Whole result referred to time result:
876 - result_(msspw,msfld,thisTimeSeg(msobs,msscan)).reference(timeResult_(msspw,msfld,thisTimeSeg(msobs,msscan)));
877 - resFlag_(msspw,msfld,thisTimeSeg(msobs,msscan)).reference(timeResFlag_(msspw,msfld,thisTimeSeg(msobs,msscan)));
900 + result_(msspw,msfld,scanOrObsInt).reference(timeResult_(msspw,msfld,scanOrObsInt));
901 + resFlag_(msspw,msfld,scanOrObsInt).reference(timeResFlag_(msspw,msfld,scanOrObsInt));
878 902
879 903 // Detect if obs or fld changed, and cal is obs- or fld-dep
880 904 Bool diffobsfld(false);
881 905 diffobsfld|=(byField_ && msfld!=lastFld_(msspw)); // field-dep, and field changed
882 906 diffobsfld|=(byObs_ && msobs!=lastObs_(msspw)); // obs-dep, and obs changed
883 907 diffobsfld|=(byScan_ && msscan!=lastScan_(msspw)); // scan-dep, and scan changed
884 908 newcal|=diffobsfld; // update newcal for return
885 909
886 910 // Remember for next pass
887 911 lastFld_(msspw)=msfld;
987 1011 default: {
988 1012 break;
989 1013 }
990 1014 }
991 1015 }
992 1016
993 1017 // Use the msspw-specific one!
994 1018 freqInterpMethod_=freqInterpMethodVec_(msspw);
995 1019
996 1020 Bool newcal(false);
997 - IPosition ip(4,0,msspw,msfld,thisTimeSeg(msobs,msscan));
1021 + Int scanOrObsInt=thisTimeSeg(msobs,msscan);
1022 + IPosition ip(4,0,msspw,msfld,scanOrObsInt);
1023 +
998 1024 // Loop over _output_ antennas
999 1025 for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1000 1026 // Call time interpolation calculation; resample in freq if new
1001 1027 // (fills timeResult_/timeResFlag_ implicitly)
1002 1028 ip(0)=iMSElem;
1003 1029 if (!tI_(ip)) {
1004 - // if (iMSElem==0) cout << "Flagging: " << ip << endl;
1030 + //if (iMSElem==0) cout << "Flagging: " << ip << endl;
1031 +
1032 + // casa log post
1033 + if (!tIMissingLogged_(ip)) { // only if these coords not logged before
1034 + LogIO log;
1035 + ostringstream msg;
1036 + String scanOrObsStr=(byScan_ ? "scan" : "obs");
1037 + msg << "MS " << scanOrObsStr << "=" << scanOrObsInt;
1038 + if (byField_) msg << " in fld=" << msfld;
1039 + msg << ",spw=" << msspw
1040 + << ",ant=" << iMSElem
1041 + << " is selected for processing, but has no available calibration in " << ctname_
1042 + << " as mapped, and will be flagged.";
1043 + log << msg.str() << LogIO::WARN;
1044 + tIMissingLogged_(ip)=true;
1045 + }
1046 +
1047 +
1005 1048 newcal=true;
1006 1049 }
1007 1050 else {
1008 -
1051 +
1009 1052 if (tI_(ip)->interpolate(time)) {
1010 1053 // Resample in frequency
1011 - Matrix<Float> fR(freqResult_(msspw,msfld,thisTimeSeg(msobs,msscan)).xyPlane(iMSElem));
1012 - Matrix<Bool> fRflg(freqResFlag_(msspw,msfld,thisTimeSeg(msobs,msscan)).xyPlane(iMSElem));
1013 - Matrix<Float> tR(timeResult_(msspw,msfld,thisTimeSeg(msobs,msscan)).xyPlane(iMSElem));
1014 - Matrix<Bool> tRflg(timeResFlag_(msspw,msfld,thisTimeSeg(msobs,msscan)).xyPlane(iMSElem));
1054 + Matrix<Float> fR(freqResult_(msspw,msfld,scanOrObsInt).xyPlane(iMSElem));
1055 + Matrix<Bool> fRflg(freqResFlag_(msspw,msfld,scanOrObsInt).xyPlane(iMSElem));
1056 + Matrix<Float> tR(timeResult_(msspw,msfld,scanOrObsInt).xyPlane(iMSElem));
1057 + Matrix<Bool> tRflg(timeResFlag_(msspw,msfld,scanOrObsInt).xyPlane(iMSElem));
1015 1058 resampleInFreq(fR,fRflg,freq,tR,tRflg,freqIn_(spwMap_(msspw)));
1016 1059
1017 1060 // Calibration is new
1018 1061 newcal=true;
1019 1062 }
1020 1063 }
1021 1064 }
1022 1065
1023 1066 // Whole result referred to freq result:
1024 - result_(msspw,msfld,thisTimeSeg(msobs,msscan)).reference(freqResult_(msspw,msfld,thisTimeSeg(msobs,msscan)));
1025 - resFlag_(msspw,msfld,thisTimeSeg(msobs,msscan)).reference(freqResFlag_(msspw,msfld,thisTimeSeg(msobs,msscan)));
1067 + result_(msspw,msfld,scanOrObsInt).reference(freqResult_(msspw,msfld,scanOrObsInt));
1068 + resFlag_(msspw,msfld,scanOrObsInt).reference(freqResFlag_(msspw,msfld,scanOrObsInt));
1026 1069
1027 1070 // Detect if obs or fld changed, and cal is obs- or fld-dep
1028 1071 Bool diffobsfld(false);
1029 1072 diffobsfld|=(byField_ && msfld!=lastFld_(msspw)); // field-dep, and field changed
1030 1073 diffobsfld|=(byObs_ && msobs!=lastObs_(msspw)); // obs-dep, and obs changed
1031 1074 diffobsfld|=(byScan_ && msscan!=lastScan_(msspw)); // scan-dep, and scan changed
1032 1075 newcal|=diffobsfld; // update newcal for return
1033 1076
1034 1077 /*
1035 1078 if (newcal) {
1036 1079 Double t0(86400.0*floor(time/86400.0));
1037 1080 cout << boolalpha
1038 1081 << "fld="<<msfld
1039 - << " obs="<<thisTimeSeg(msobs,msscan)
1082 + << " obs="<<scanOrObsInt
1040 1083 << " spw="<<msspw
1041 1084 << " time="<< time-t0
1042 1085 << " diffobsfld=" << diffobsfld
1043 1086 << " new=" << newcal
1044 1087 << " tI_(ip)=" << tI_(ip)
1045 1088 << " chan="<< nMSChan/2
1046 - << " result=" << result_(msspw,msfld,thisTimeSeg(msobs,msscan))(0,nMSChan/2,0)
1047 - << " addr=" << &result_(msspw,msfld,thisTimeSeg(msobs,msscan))(0,nMSChan/2,0)
1089 + << " result=" << result_(msspw,msfld,scanOrObsInt)(0,nMSChan/2,0)
1090 + << " addr=" << &result_(msspw,msfld,scanOrObsInt)(0,nMSChan/2,0)
1048 1091 << endl;
1049 1092 }
1050 1093 */
1051 1094 // Remember for next pass
1052 1095 lastFld_(msspw)=msfld;
1053 1096 lastObs_(msspw)=msobs;
1054 1097 lastScan_(msspw)=msscan;
1055 1098
1056 1099 return newcal;
1057 1100 }
1199 1242 Int iobs = (byObs_ ? ctiter.thisObs() : 0); // use 0 if not slicing by obs
1200 1243 iobs = (byScan_ ? ctiter.thisScan() : iobs); // XXX perscan overrides perobs?
1201 1244 IPosition ip(4,iant,ispw,ifld,iobs);
1202 1245 ctSlices_(ip)= new NewCalTable(ctiter.table());
1203 1246 spwInOK_(ispw)=(spwInOK_(ispw) || ctSlices_(ip)->nrow()>0);
1204 1247 ctiter.next();
1205 1248 }
1206 1249 break;
1207 1250 }
1208 1251 }
1209 -
1252 +
1210 1253 }
1211 1254
1212 1255 // Initialize by iterating over the supplied table
1213 1256 void CTPatchedInterp::makeInterpolators() {
1214 1257
1215 1258 if (CTPATCHEDINTERPVERB) cout << " CTPatchedInterp::initialize()" << endl;
1216 1259
1217 - // cal table name for messages
1218 - Path pathname(ct_.tableName());
1219 - String tabname=pathname.baseName().before(".tempMem");
1220 -
1260 + // Save caltable name for log messages
1261 + ctname_=Path(ct_.antenna().tableName().before(".tempMem")).baseName();
1262 +
1221 1263 // Size/initialize interpolation engines
1222 1264 IPosition tIsize(4,nMSElem_,nMSSpw_,nMSFld_,nMSTimeSeg_);
1223 1265 tI_.resize(tIsize);
1224 1266 tI_.set(NULL);
1225 1267 tIdel_.resize(tIsize);
1226 1268 tIdel_.set(false);
1227 -
1228 - Record summary;
1229 - std::set<Int> scans;
1230 - if (msmc_) {
1231 - summary = msmc_->msmd().getSummary();
1232 - scans = msmc_->msmd().getScanNumbers(0, 0);
1233 - }
1269 + tIMissingLogged_.resize(tIsize);
1270 + tIMissingLogged_.set(false);
1234 1271
1235 1272 Bool reportBadSpw(false);
1236 1273 for (Int iMSTimeSeg=0;iMSTimeSeg<nMSTimeSeg_;++iMSTimeSeg) {
1237 1274 for (Int iMSFld=0;iMSFld<nMSFld_;++iMSFld) {
1238 1275
1239 1276 if (altFld_(iMSFld)==iMSFld) { // i.e., this field has its own solutions
1240 1277
1241 1278 // cout << "Making interpolators for " << iMSFld << " (mapped from " << fldMap_(iMSFld) << ")" << endl;
1242 1279
1243 1280 std::set<uInt> spws;
1273 1310 if (!ict.isNull()) {
1274 1311 tI_(tIip)=(*cttifactoryptr_)(ict,timeType_,tR,tRf);
1275 1312 tIdel_(tIip)=true;
1276 1313 }
1277 1314 }
1278 1315 else {
1279 1316 // the required ct slice is empty, so arrange to flag it
1280 1317 tI_(tIip)=NULL;
1281 1318 tR.set(0.0);
1282 1319 tRf.set(true);
1283 -
1284 - // Get Obs Id from summary, from third set of keys get scan(s)
1285 - // Get antenna from scan and fields for scan
1286 - // Logic check -->
1287 - // Do you have obs/scan -> is your field in scans? Spw in field? antenna in scan?
1288 - // If all true print below, else don't print anything
1289 -
1290 - LogIO log;
1291 - ostringstream msg;
1292 -
1293 - if (spws.find(iMSSpw) != spws.end()) {
1294 - if (byScan_) {
1295 - if (scans.find(iMSTimeSeg) != scans.end()) {
1296 -
1297 - // casa log post
1298 - msg << "IF SELECTED, MS scan=" << iMSTimeSeg
1299 - << " in fld=" << iMSFld
1300 - << ",spw=" << iMSSpw
1301 - << ",ant=" << iMSElem
1302 - << " cannot be calibrated by " << tabname
1303 - << " as mapped, and will be flagged in this process.";
1304 - log << msg.str() << LogIO::WARN;
1305 - }
1306 - } else {
1307 - if (summary.isDefined("observationID=" + String::toString(iMSTimeSeg))) {
1308 -
1309 - // casa log post
1310 - msg << "MS obs=" << iMSTimeSeg
1311 - << ",fld=" << iMSFld
1312 - << ",spw=" << iMSSpw
1313 - << ",ant=" << iMSElem
1314 - << " cannot be calibrated by " << tabname
1315 - << " as mapped, and will be flagged in this process.";
1316 - log << msg.str() << LogIO::WARN;
1317 - }
1318 - }
1319 - }
1320 1320 }
1321 1321 } // iMSElem
1322 1322 } // spwOK
1323 1323 else
1324 1324 reportBadSpw=true;
1325 1325 } // iMSSpw
1326 1326
1327 1327 } // not re-using
1328 - else {
1328 + else { // re-using
1329 1329 // Point to an existing interpolator group
1330 1330 Int thisAltFld=altFld_(iMSFld);
1331 1331
1332 1332 // cout << "Reusing interpolators from " << thisAltFld << " for " << iMSFld << " (mapped to " << fldMap_(iMSFld) << ")" << endl;
1333 1333
1334 1334 for (Int iMSSpw=0;iMSSpw<nMSSpw_;++iMSSpw) {
1335 1335 timeResult_(iMSSpw,iMSFld,iMSTimeSeg).reference(timeResult_(iMSSpw,thisAltFld,iMSTimeSeg));
1336 1336 timeResFlag_(iMSSpw,iMSFld,iMSTimeSeg).reference(timeResFlag_(iMSSpw,thisAltFld,iMSTimeSeg));
1337 1337 for (Int iMSElem=0;iMSElem<nMSElem_;++iMSElem) {
1338 1338 IPosition tIip0(4,iMSElem,iMSSpw,iMSFld,iMSTimeSeg),tIip1(4,iMSElem,iMSSpw,thisAltFld,iMSTimeSeg);
1339 1339 tI_(tIip0)=tI_(tIip1);
1340 + // if (!tI_(tIip0) && iMSElem==0) cout << "ouch---------------------" << "iMSTimeSeg="<<iMSTimeSeg<<" iMSFld="<<iMSFld<<" spw="<< iMSSpw << " ant="<<iMSElem<< endl;
1340 1341 }
1341 1342 }
1342 1343 }
1343 1344 } // iMSFld
1344 1345 } // iMSTimeSeg
1345 1346
1346 1347
1347 1348 if (reportBadSpw) {
1348 - cout << "The following MS spws have no corresponding cal spws in " << tabname << ": ";
1349 + cout << "The following MS spws have no corresponding cal spws in " << ctname_ << ": ";
1349 1350 for (Int iMSSpw=0;iMSSpw<nMSSpw_;++iMSSpw)
1350 1351 // (spwmap applied in spwOK method)
1351 1352 if (!this->spwOK(iMSSpw)) cout << iMSSpw << " ";
1352 1353 cout << endl;
1353 1354 }
1354 1355
1355 1356 }
1356 1357
1357 1358
1358 1359

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut