//# CalCache.cc: Specialized PlotMSCache for filling CalTables //# Copyright (C) 2009 //# Associated Universities, Inc. Washington DC, USA. //# //# This library is free software; you can redistribute it and/or modify it //# under the terms of the GNU Library General Public License as published by //# the Free Software Foundation; either version 2 of the License, or (at your //# option) any later version. //# //# This library is distributed in the hope that it will be useful, but WITHOUT //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public //# License for more details. //# //# You should have received a copy of the GNU Library General Public License //# along with this library; if not, write to the Free Software Foundation, //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA. //# //# Correspondence concerning AIPS++ should be addressed as follows: //# Internet email: aips2-request@nrao.edu. //# Postal address: AIPS++ Project Office //# National Radio Astronomy Observatory //# 520 Edgemont Road //# Charlottesville, VA 22903-2475 USA //# #include <plotms/Data/CalCache.h> #include <plotms/Data/PlotMSAtm.h> #include <plotms/Data/PlotMSCTAverager.h> #include <plotms/Data/PlotMSIndexer.h> #include <plotms/PlotMS/PlotMS.h> #include <plotms/Threads/ThreadCommunication.h> #include <casa/OS/Timer.h> #include <casa/OS/HostInfo.h> #include <casa/OS/Memory.h> #include <casa/Quanta/MVTime.h> #include <casa/System/Aipsrc.h> #include <casa/Utilities/Sort.h> #include <casa/Arrays/ArrayMath.h> #include <graphics/GenericPlotter/Plotter.h> #include <ms/MSSel/MSSelectionTools.h> #include <synthesis/CalTables/CTColumns.h> #include <synthesis/MeasurementComponents/VisCalGlobals.h> #include <synthesis/MeasurementComponents/BPoly.h> #include <synthesis/MeasurementComponents/GSpline.h> #include <tables/Tables/Table.h> using namespace casacore; namespace casa { // Define external CLIC solvers #define cheb cheb_ extern "C" { void cheb(Int*, Double*, Double*, Int*); } CalCache::CalCache(PlotMSApp* parent): PlotMSCacheBase(parent), ci_p(nullptr), wci_p(nullptr), basis_("unknown"), parsAreComplex_(False), divZero_(False), msname_("") { } CalCache::~CalCache() {} String CalCache::polname(Int ipol) { if (polnRatio_) return "/"; if (basis_=="Linear") return ( (ipol%2==0) ? String("X") : String("Y") ); else if (basis_=="Circular") return ( (ipol%2==0) ? String("R") : String("L") ); else { // "unknown", or antenna positions if (calType_=="KAntPos Jones") { switch(ipol) { case 0: return "X"; case 1: return "Y"; case 2: return "Z"; default: return (String::toString(ipol)); } } else { return ( String::toString(ipol) ); } } } void CalCache::setFilename(String filename) { filename_ = filename; Table tab(filename); calType_= tab.tableInfo().subType(); if ((calType_=="T Jones") && (tab.keywordSet().isDefined("CAL_DESC"))) { throw AipsError(calType_ + " tables in the old cal table format are unsupported in plotms."); } if (tab.keywordSet().isDefined("MSName")) { casacore::String msname = tab.keywordSet().asString("MSName"); setMSName(msname); } } void CalCache::setMSName(String msname) { // set msname_ (with path) if valid Path filepath(filename_); String path(filepath.dirName()); if (!path.empty()) path += "/"; msname_ = path + msname; } //********************************* // protected method implementations void CalCache::loadIt(vector<PMS::Axis>& loadAxes, vector<PMS::DataColumn>& /*loadData*/, ThreadCommunication* thread) { // this also sets calType_: setFilename(filename_); logLoad("Plotting a " + calType_ + " calibration table."); // Trap unsupported cal types, averaging, transforms, poln ratio if ((calType_[0]=='X') && calType_.contains("Mueller")) { throw AipsError("Cal table type " + calType_ + " is unsupported in plotms. Please continue to use plotcal."); } // No averaging for BPOLY and GSPLINE if (((calType_=="BPOLY") || (calType_=="GSPLINE")) && averaging().anyAveraging()) { throw AipsError("Averaging not supported for cal table type " + calType_); } // No averaging with channel selection if (averaging().anyAveraging() && selection_.spw().contains(":")) { throw AipsError("Averaging not supported with channel selection for calibration tables"); } // Warn that transformations will be ignored if (transformations().anyTransform()) { logWarn("CalCache::loadIt", "Transformations ignored: not supported for calibration tables"); } // poln ratio polnRatio_ = false; if (selection_.corr()=="/") { if ((calType_=="BPOLY") || (calType_[0] == 'T') || ((calType_[0] == 'F') && !calType_.startsWith("Fringe"))) { throw(AipsError("Polarization ratio plots not supported for " + calType_ + " tables.")); } else { polnRatio_ = true; } } antnames_.resize(); stanames_.resize(); antstanames_.resize(); fldnames_.resize(); positions_.resize(); vector<PMS::DataColumn> loadData; loadData.assign(loadAxes.size(), PMS::DEFAULT_DATACOLUMN); if (calType_=="BPOLY") { loadBPoly(loadAxes, loadData, thread); } else if (calType_=="GSPLINE") { checkAxes(loadAxes); // check for invalid axis before proceeding loadGSpline(loadAxes, loadData, thread); } else { loadNewCalTable(loadAxes, loadData, thread); } } // ======================== NewCalTable ========================== void CalCache::loadNewCalTable(vector<PMS::Axis>& loadAxes, vector<PMS::DataColumn>& loadData, ThreadCommunication* thread) { // Load requested axes from NewCalTable // Get various names, properties TableLock lock(TableLock::AutoNoReadLocking); NewCalTable* ct = new NewCalTable(filename_, lock, Table::Old, Table::Plain); basis_ = ct->polBasis(); parsAreComplex_ = ct->isComplex(); ROCTColumns ctCol(*ct); antnames_ = ctCol.antenna().name().getColumn(); stanames_ = ctCol.antenna().station().getColumn(); positions_ = ctCol.antenna().position().getColumn(); nAnt_ = ctCol.antenna().nrow(); fldnames_ = ctCol.field().name().getColumn(); antstanames_ = antnames_ + String("@") + stanames_; // Apply selection to get selected cal table NewCalTable* selct = new NewCalTable(); selection_.apply(*ct, *selct); if (selct->nrow() == 0) { throw(AipsError("Selection resulted in zero rows")); } PlotMSAveraging pmsAveraging(averaging()); casacore::Bool readonly(True); setUpCalIter(*selct, pmsAveraging, readonly); ci_p->reset(); parshape_ = ci_p->flag().shape(); // Size cache arrays based on number of chunks if (pmsAveraging.anyAveraging()) { // Use PlotMSCTAverager casacore::Vector<int> nIterPerAve; // number of chunks per average countChunks(*ci_p, pmsAveraging, nIterPerAve, loadAxes, loadData, thread); if (!userCanceled_) { loadCalChunks(*ci_p, pmsAveraging, nIterPerAve, loadAxes, thread); } } else { countChunks(*ci_p, loadAxes, loadData, thread); loadCalChunks(*ci_p, loadAxes, thread); } // delete NCT and iter to release table locks if (ct != nullptr) { delete ct; ct = nullptr; } if (selct != nullptr) { delete selct; selct = nullptr; } if (ci_p != nullptr) { delete ci_p; ci_p = nullptr; } } void CalCache::setUpCalIter( NewCalTable& selct, PlotMSAveraging& pmsAveraging, Bool readonly) { // Set up cal table iterator for counting and loading chunks // Order of sort columns depends on averaging options Int nsortcol(3 + Int(!pmsAveraging.scan())), col(0); Block<String> columns(nsortcol); if (!pmsAveraging.scan()) { columns[col++]="SCAN_NUMBER"; } if (!pmsAveraging.field()) { columns[col++]="FIELD_ID"; } if (!pmsAveraging.spw()) { columns[col++]="SPECTRAL_WINDOW_ID"; } columns[col++]="TIME"; if (pmsAveraging.field()) { columns[col++]="FIELD_ID"; } if (pmsAveraging.spw()) { columns[col++]="SPECTRAL_WINDOW_ID"; } sortColumns_ = columns; if (readonly) { // Readonly version, for caching ci_p = new ROCTIter(selct, columns); wci_p = nullptr; } else { // Writable, e.g. for flagging wci_p = new CTIter(selct, columns); ci_p = wci_p; // const access } } void CalCache::countChunks(ROCTIter& ci, vector<PMS::Axis>& loadAxes, vector<PMS::DataColumn>& loadData, ThreadCommunication* thread) { // for NewCalTable (no averaging) // loadData not applicable but needed for setCache() if (thread) { thread->setStatus("Establishing cache size. Please wait..."); thread->setAllowedOperations(false,false,false); } // Iterate and count number of chunks. int chunk(0); ci.reset(); while (!ci.pastEnd()) { ++chunk; ci.next0(); } setCache(chunk, loadAxes, loadData); } void CalCache::countChunks(ROCTIter& ci, PlotMSAveraging& pmsAveraging, Vector<int>& nIterPerAve, vector<PMS::Axis>& loadAxes, vector<PMS::DataColumn>& loadData, ThreadCommunication* thread) { // for NewCalTable (with averaging) // loadData not applicable but needed for setCache() if (pmsAveraging.time() || pmsAveraging.baseline() || pmsAveraging.antenna() || pmsAveraging.spw() || pmsAveraging.scalarAve()) { // Set number of iterations per averaged chunk (nIterPerAve) // Allow user to cancel if (thread) { thread->setStatus("Establishing cache size. Please wait..."); thread->setAllowedOperations(false,false,true); } bool debug(false); Bool combScan(pmsAveraging.scan()); Bool combField(pmsAveraging.field()); Bool combSpw(pmsAveraging.spw()); // Keep track of time and averaging interval Double thistime(0.0), avetime0(-1.0); Double interval(0.0); if (pmsAveraging.time()) { interval = pmsAveraging.timeValue(); } // Keep track of other boundaries Int thisscan(-1),lastscan(-1); Int thisfield(-1), lastfield(-1); Int thisspw(-1),lastspw(-1); Int thisobsid(-1),lastobsid(-1); // Averaging stats Int chunk(0); Int maxAveNRows(0); nIterPerAve.resize(100); nIterPerAve = 0; Int nAveInterval(-1); ci.reset(); while (!ci.pastEnd()) { // If a thread is given, check if the user canceled. if (thread != nullptr) { if (thread->wasCanceled()) { dataLoaded_ = false; userCanceled_ = true; } else { // else users think it's hung... if ((chunk % 100) == 0) { thread->setProgress(chunk/100); } } } thistime = ci.thisTime(); thisscan = ci.thisScan(); thisfield = ci.thisField(); thisspw = ci.thisSpw(); thisobsid = ci.thisObs(); if ( ((thistime - avetime0) > interval) || // past avgtime interval ((thistime - avetime0) < 0.0) || // negative timestep (!combScan && (thisscan != lastscan)) || // new scan (!combField && (thisfield != lastfield)) || // new field (!combSpw && (thisspw != lastspw)) || // new spw (thisobsid != lastobsid) || // new obs id (nAveInterval == -1)) { // first interval // New averaging interval if (debug) { stringstream ss; ss << "--------------------------------\n"; ss << "New averaging interval\n"; ss << "time elapsed=" << ((thistime - avetime0) > interval) << " " << " neg step=" << ((thistime - avetime0) < 0.0) << " " << " scan=" << (!combScan && (thisscan != lastscan)) << " " << " spw=" << (!combSpw && (thisspw != lastspw)) << " " << " field=" << (!combField && (thisfield != lastfield)) << " " << " obsid=" << (thisobsid!=lastobsid) << " " << " first=" << (nAveInterval == -1) << "\n"; logInfo("count_chunks", ss.str()); } // If we have accumulated enough info, reset the ave'd row counter maxAveNRows = 0; nAveInterval++; if (debug) { stringstream ss; ss << "ave = " << nAveInterval << "\n"; logInfo("count_chunks", ss.str()); } // increase size of nIterPerAve array, if needed if (nIterPerAve.nelements() < uInt(nAveInterval + 1)) { nIterPerAve.resize(nIterPerAve.nelements()+100, true); } // initialize next ave interval nIterPerAve(nAveInterval) = 0; avetime0 = thistime; // first timestamp in this averaging interval } // Keep track of the maximum # of rows that might get averaged maxAveNRows = max(maxAveNRows, ci.nrow()); // Increment chunk-per-average count for current solution nIterPerAve(nAveInterval)++; if (debug) { stringstream ss; ss << "Completed chunk=" << chunk << "\n"; ss << "time=" << thistime << " "; ss << "scan" << thisscan << " "; ss << "fieldId=" << thisfield << " "; ss << "spw=" << thisspw << " "; ss << "obsId=" << thisobsid << "\n"; logInfo("count_chunks", ss.str()); } // Store last values for next iteration lastscan = thisscan; lastfield = thisfield; lastspw = thisspw; lastobsid = thisobsid; ci.next(); chunk++; } Int nAve(nAveInterval + 1); nIterPerAve.resize(nAve, True); setCache(nAve, loadAxes, loadData); // initialize cache size, nChunk_ if (debug) { stringstream ss; ss << "nIterPerAve = " << nIterPerAve << "\n"; ss << "Found " << nChunk_ << " chunks." << endl; logInfo("count_chunks", ss.str()); } } else { // In-row (channel) averaging does not change number of chunks countChunks(ci, loadAxes, loadData, thread); // Each chunk can be averaged separately nIterPerAve.resize(nChunk_); nIterPerAve = 1; } } void CalCache::loadCalChunks(ROCTIter& ci, const vector<PMS::Axis> loadAxes, ThreadCommunication* thread) { // for NewCalTable // permit cancel in progress meter: if(thread != nullptr) thread->setAllowedOperations(false,false,true); logLoad("Loading chunks......"); Int chunk(0), lastscan(0), thisscan(0), lastspw(-1), thisspw(0); chshapes_.resize(4,nChunk_); goodChunk_.resize(nChunk_); goodChunk_.set(False); double progress; // Reset iterator ci.reset(); while (!ci.pastEnd()) { // If a thread is given, check if the user canceled. if (thread != nullptr && thread->wasCanceled()) { dataLoaded_ = false; return; } // If a thread is given, update it. if (thread != nullptr && (nChunk_ <= (int)THREAD_SEGMENT || chunk % THREAD_SEGMENT == 0)) { thread->setStatus("Loading chunk " + String::toString(chunk) + " / " + String::toString(nChunk_) + "."); } // Discern npol/nchan shape IPosition pshape(ci.flag().shape()); // Use viscal to determine nPol String pol = selection_.corr(); String paramAxis = toVisCalAxis(PMS::AMP); size_t nPol; if (polnRatio_) { // pick one! nPol = getParSlice(paramAxis, "R").length(); } else { nPol = getParSlice(paramAxis, pol).length(); } // Apply channel selection to get nChan size_t nChan(pshape[1]); casacore::Matrix<casacore::Int> selectedChans = selection_.getSelectedChannels(); std::vector<casacore::Slice> chansel; if (selectedChans.empty()) { chansel.push_back(Slice()); } else { size_t nChanSelectedThisSpw(0); for (size_t i = 0; i < selectedChans.nrow(); ++i) { auto row = selectedChans.row(i); if (row(0) == ci.thisSpw()) { casacore::Slice chanSlice(row(1), row(2), row(3), false); chansel.push_back(chanSlice); nChanSelectedThisSpw += chanSlice.length(); } } if (nChanSelectedThisSpw) { nChan = nChanSelectedThisSpw; } } // Cache the data shapes chshapes_(0,chunk) = nPol; chshapes_(1,chunk) = nChan; chshapes_(2,chunk) = ci.nrow(); chshapes_(3,chunk) = nAnt_; goodChunk_(chunk) = True; for(unsigned int i = 0; i < loadAxes.size(); i++) { loadCalAxis(ci, chunk, loadAxes[i], pol, chansel); // print atm stats once per scan if (loadAxes[i]==PMS::ATM || loadAxes[i]==PMS::TSKY) { thisscan = ci.thisScan(); if (thisscan != lastscan) { printAtmStats(thisscan); lastscan = thisscan; } thisspw = ci.thisSpw(); if (thisspw != lastspw) { uInt vectorsize = ( loadAxes[i]==PMS::ATM ? (*atm_[chunk]).nelements() : (*tsky_[chunk]).nelements()); if (vectorsize==1) { logWarn("load_cache", "Setting " + PMS::axis(loadAxes[i]) + " for spw " + String::toString(thisspw) + " to zero because it has only one channel."); } lastspw = thisspw; } } } chunk++; ci.next(); // If a thread is given, update it. if ((thread != nullptr) && ((nChunk_ <= (int)THREAD_SEGMENT) || (chunk % THREAD_SEGMENT == 0))) { progress = ((double)chunk+1) / nChunk_; thread->setProgress((unsigned int)((progress * 100) + 0.5)); } } if (divZero_) { logWarn("CalCache::loadIt", "Caught divide-by-zero exception in ratio plots; result(s) set to 1.0 and flagged"); } } void CalCache::loadCalChunks(ROCTIter& ci, PlotMSAveraging& pmsAveraging, const casacore::Vector<int>& nIterPerAve, const std::vector<PMS::Axis> loadAxes, ThreadCommunication* thread) { // Load chunks using PlotMSCTAverager to accumulate and average chunks logLoad("Loading chunks with averaging....."); // Permit cancel in progress meter: if (thread != nullptr) { thread->setAllowedOperations(false,false,true); } // TBD: channel selection, for now select all std::vector<casacore::Slice> chansel; chansel.push_back(Slice()); // Access to header info and subtables when loading axes String partype = parsAreComplex_ ? "Complex" : "Float"; CTDesc caltabdesc(partype, msname_, calType_, basis_); chshapes_.resize(4, nChunk_); goodChunk_.resize(nChunk_); goodChunk_.set(false); String polsel(selection_.corr()); // for slicing axis data Int lastscan(0), thisscan(0); // print atm stats once per scan Int lastspw(-1), thisspw(0); // print atm warning once per spw double progress; ci.reset(); for (Int chunk = 0; chunk < nChunk_; ++chunk) { // Update progress with each chunk if ((thread != nullptr) && ((nChunk_ <= (int)THREAD_SEGMENT) || (chunk % THREAD_SEGMENT == 0))) { thread->setStatus("Loading chunk " + String::toString(chunk) + " / " + String::toString(nChunk_) + "."); progress = ((double)chunk + 1) / nChunk_; thread->setProgress((unsigned int)((progress * 100) + 0.5)); } // Set up CT Averager for each averaged chunk PlotMSCTAverager pmscta(pmsAveraging, nAnt_, parshape_(0)); if (calType_.contains("Mueller")) { pmscta.setBaselineBased(); } // Accumulate iterations into chunk Int iter(0); while (iter < nIterPerAve(chunk)) { pmscta.accumulate(ci); // Advance to next iteration unless finalize if ((iter + 1) < nIterPerAve(chunk)) { ci.next(); } ++iter; } // Finalize average pmscta.finalizeAverage(); // Get result as a memory NewCalTable with averaged main rows NewCalTable avgTable("avgcaltable.cal", caltabdesc, Table::Scratch, Table::Memory); pmscta.fillAvgCalTable(avgTable); if (avgTable.nrow() > 0) { // Attach iterator for accessor ROCTIter avgTableCti(avgTable, sortColumns_); avgTableCti.setCTColumns(ci.table()); avgTableCti.reset(); // Cache data shape IPosition avgShape(avgTableCti.flag().shape()); chshapes_(0, chunk) = avgShape(0); chshapes_(1, chunk) = avgShape(1); chshapes_(2, chunk) = avgShape(2); chshapes_(3, chunk) = nAnt_; goodChunk_(chunk) = true; // Load axes for (auto axis : loadAxes) { // Check for cancel before each axis if (thread && thread->wasCanceled()) { dataLoaded_ = false; userCanceled_ = true; goodChunk_(chunk) = false; return; } switch (axis) { case PMS::CHANNEL: { Int nchan(pmscta.nchan()); Vector<Int> chans(nchan); indgen(chans); *chan_[chunk] = chans; break; } case PMS::FREQUENCY: { Vector<Double> freqs = pmscta.avgfreq(); freqs /= 1.0e9; // GHz (*freq_[chunk]) = freqs; break; } case PMS::ATM: case PMS::TSKY: case PMS::IMAGESB: { // Use original caltable spectral window subtable for these axes Int spw = avgTableCti.thisSpw(); Int scan = avgTableCti.thisScan(); Vector<Double> freqs = pmscta.avgfreq(); casacore::Vector<casacore::Double> curve(1, 0.0); if (axis == PMS::ATM) { if (plotmsAtm_) { plotmsAtm_->calcAtmTskyCurve(curve, spw, scan, freqs); } *atm_[chunk] = curve; } else if (axis == PMS::TSKY) { if (plotmsAtm_) { plotmsAtm_->calcAtmTskyCurve(curve, spw, scan, freqs); } *tsky_[chunk] = curve; } else { if (plotmsAtm_) { plotmsAtm_->calcImageCurve(curve, spw, scan, freqs); } *imageSideband_[chunk] = curve; } break; } default: loadCalAxis(avgTableCti, chunk, axis, polsel, chansel); } if ((axis == PMS::ATM) || (axis == PMS::TSKY)) { // Print stats when scan changes thisscan = avgTableCti.thisScan(); if (thisscan != lastscan) { printAtmStats(thisscan); lastscan = thisscan; } // Print warning when one channel thisspw = avgTableCti.thisSpw(); if (thisspw != lastspw) { uInt nchan = (axis == PMS::ATM ? (*atm_[chunk]).nelements() : (*tsky_[chunk]).nelements()); if (nchan == 1) { logWarn("load_cache", "Setting " + PMS::axis(axis) + " for spw " + String::toString(thisspw) + " to zero because it has only one channel."); } lastspw = thisspw; } } } // end load axes } else { // No rows in result goodChunk_(chunk) = false; chshapes_.column(chunk) = 0; } // Advance to next chunk ci.next(); } // chunk loop } void CalCache::loadCalAxis(ROCTIter& cti, casacore::Int chunk, PMS::Axis axis, casacore::String& pol, std::vector<casacore::Slice>& chansel) { // for NewCalTable // Get polarization selection slice Slice parSlice1 = Slice(); Slice parSlice2 = Slice(); if (PMS::axisNeedsCalSlice(axis)) { String calAxis = toVisCalAxis(axis); if (polnRatio_) { parSlice1 = getParSlice(calAxis, "R"); parSlice2 = getParSlice(calAxis, "L"); } else { parSlice1 = getParSlice(calAxis, pol); } } switch(axis) { case PMS::SCAN: // assumes scan unique scan_(chunk) = cti.thisScan(); break; case PMS::FIELD: field_(chunk) = cti.thisField(); break; case PMS::TIME: // assumes time unique time_(chunk) = cti.thisTime(); break; case PMS::TIME_INTERVAL: // assumes timeInterval unique in cti chunk timeIntr_(chunk) = cti.thisInterval(); break; case PMS::SPW: spw_(chunk) = cti.thisSpw(); break; case PMS::CHANNEL: { casacore::Vector<casacore::Int> channels; cti.chan(channels); // Apply channel slices casacore::Vector<casacore::Int> selectedChans; for (auto& chan_slice : chansel) { casacore::Vector<casacore::Int> chanSlice = channels(chan_slice); ConcatArrays<casacore::Int>(selectedChans, chanSlice); } *chan_[chunk] = selectedChans; break; } case PMS::FREQUENCY: { // TBD: Convert freq to desired frame casacore::Vector<casacore::Double> freqs; cti.freq(freqs); // Apply channel slices casacore::Vector<casacore::Double> selectedFreqs; for (auto& chan_slice : chansel) { casacore::Vector<casacore::Double> freqSlice = freqs(chan_slice); ConcatArrays<casacore::Double>(selectedFreqs, freqSlice); } *freq_[chunk] = selectedFreqs / 1.0e9; // in GHz break; } /* case PMS::VELOCITY: { // Convert freq in the vb to velocity vbu_.toVelocity(*vel_[chunk], vb, transformations_.frame(), MVFrequency(transformations_.restFreqHz()), transformations_.veldef()); (*vel_[chunk]) /= 1.0e3; // in km/s break; } */ case PMS::CORR: { corr_[chunk]->resize(chshapes_(0,chunk)); if (pol=="" || pol=="RL" || pol=="XY") { indgen(*corr_[chunk]); } else if (pol== "R" || pol=="X") { corr_[chunk]->resize(1); corr_[chunk]->set(0); } else if (pol== "L" || pol=="Y") { corr_[chunk]->resize(1); corr_[chunk]->set(1); } else if (pol=="/") { corr_[chunk]->resize(1); corr_[chunk]->set(-1); // ??? } break; } case PMS::ANTENNA1: *antenna1_[chunk] = cti.antenna1(); break; case PMS::ANTENNA2: *antenna2_[chunk] = cti.antenna2(); break; case PMS::BASELINE: { baseline_[chunk]->resize(cti.nrow()); Vector<Int> bl(*baseline_[chunk]); if (averaging().baseline()) { bl.set(0); } else { Vector<Int> a1(cti.antenna1()); Vector<Int> a2(cti.antenna2()); if (allEQ(a2, -1)) { bl = a1; } else { for (Int irow = 0; irow < cti.nrow(); ++irow) { // Same hash as in MSCache: if (a1(irow) < 0) a1(irow) = chshapes_(3, 0); if (a2(irow) < 0) a2(irow) = chshapes_(3, 0); bl(irow) = (chshapes_(3,0)+1)*a1(irow) - (a1(irow)*(a1(irow) + 1))/2 + a2(irow); } } } break; } case PMS::U: { *u_[chunk] = cti.uvw().row(0); break; } case PMS::V: { *v_[chunk] = cti.uvw().row(1); break; } case PMS::W: { *w_[chunk] = cti.uvw().row(2); break; } case PMS::UVDIST: { Array<Double> u(cti.uvw().row(0)); Array<Double> v(cti.uvw().row(1)); *uvdist_[chunk] = sqrt(u*u+v*v); break; } case PMS::UVDIST_L: { Array<Double> u(cti.uvw().row(0)); Array<Double> v(cti.uvw().row(1)); Vector<Double> uvdistM = sqrt(u*u + v*v); uvdistM /= C::c; casacore::Vector<casacore::Double> freqs; cti.freq(freqs); auto nrow = cti.nrow(); uvdistL_[chunk]->resize(cti.nchan(), nrow); Vector<Double> uvrow; for (int irow = 0; irow < nrow; ++irow) { uvrow.reference(uvdistL_[chunk]->column(irow)); uvrow.set(uvdistM(irow)); uvrow *= freqs; } break; } case PMS::UWAVE: { Vector<Double> uM(cti.uvw().row(0)); uM /= C::c; casacore::Vector<casacore::Double> freqs; cti.freq(freqs); auto nrow = cti.nrow(); uwave_[chunk]->resize(cti.nchan(), nrow); Vector<Double> urow; for (int irow = 0; irow < nrow; ++irow) { urow.reference(uwave_[chunk]->column(irow)); urow.set(uM(irow)); urow *= freqs; } break; } case PMS::VWAVE: { Vector<Double> vM(cti.uvw().row(1)); vM /= C::c; casacore::Vector<casacore::Double> freqs; cti.freq(freqs); auto nrow = cti.nrow(); vwave_[chunk]->resize(cti.nchan(), nrow); Vector<Double> vrow; for (int irow = 0; irow < nrow; ++irow) { vrow.reference(vwave_[chunk]->column(irow)); vrow.set(vM(irow)); vrow *= freqs; } break; } case PMS::WWAVE: { Vector<Double> wM(cti.uvw().row(2)); wM /= C::c; casacore::Vector<casacore::Double> freqs; cti.freq(freqs); auto nrow = cti.nrow(); wwave_[chunk]->resize(cti.nchan(), nrow); Vector<Double> wrow; for (int irow = 0; irow < nrow; ++irow) { wrow.reference(wwave_[chunk]->column(irow)); wrow.set(wM(irow)); wrow *= freqs; } break; } case PMS::ANTPOS: { if (!calType_.startsWith("KAntPos")) { throw(AipsError( "ANTPOS has no meaning for this table")); } Cube<Float> fArray = cti.fparam(); Array<Float> selectedAntPos; for (auto& chan_slice : chansel) { Array<Float> antposSlice = fArray(parSlice1, chan_slice, Slice()); ConcatArrays<casacore::Float>(selectedAntPos, antposSlice); } *antpos_[chunk] = selectedAntPos; break; } case PMS::GAMP: case PMS::AMP: { Cube<Complex> cArray; Cube<Float> fArray; bool isComplex(parsAreComplex()); if (isComplex) { cArray = cti.cparam(); } else { fArray = cti.fparam(); } // Take each channel slice and concat into selectedAmp Array<Float> selectedAmp; for (auto& chan_slice : chansel) { Array<Float> ampSlice; if (polnRatio_) { if (isComplex) { ampSlice = amplitude(cArray(parSlice1, chan_slice, Slice()) / cArray(parSlice2, chan_slice, Slice())); checkRatioArray(ampSlice, chunk); } else { if (calType_ == "Fringe Jones") { // subtract ampSlice = fArray(parSlice1, chan_slice, Slice()) - fArray(parSlice2, chan_slice, Slice()); } else { ampSlice = fArray(parSlice1, chan_slice, Slice()) / fArray(parSlice2, chan_slice, Slice()); checkRatioArray(ampSlice, chunk); } } } else { if (isComplex) { ampSlice = amplitude( cArray(parSlice1, chan_slice, Slice())); } else { ampSlice = fArray(parSlice1, chan_slice, Slice()); } } ConcatArrays<casacore::Float>(selectedAmp, ampSlice); } if (calType_[0] == 'F') { // F Jones TEC table selectedAmp /= Float(1e+16); } *amp_[chunk] = selectedAmp; break; } case PMS::GPHASE: case PMS::PHASE: { if (parsAreComplex()) { Cube<Complex> cArray = cti.cparam(); Array<Float> selectedPhase; for (auto& chan_slice : chansel) { Array<Float> phaseSlice; if (polnRatio_) { phaseSlice = phase(cArray(parSlice1, chan_slice, Slice()) / cArray(parSlice2, chan_slice, Slice())); checkRatioArray(phaseSlice, chunk); } else { phaseSlice = phase(cArray(parSlice1, chan_slice, Slice())); } ConcatArrays<casacore::Float>(selectedPhase, phaseSlice); } *pha_[chunk] = selectedPhase * Float(180.0/C::pi); } else if (calType_ == "Fringe Jones") { Cube<Float> fArray = cti.fparam(); Array<Float> selectedPhase; for (auto& chan_slice : chansel) { Array<Float> phaseSlice; if (polnRatio_) { phaseSlice = fArray(parSlice1, chan_slice, Slice()) - fArray(parSlice2, chan_slice, Slice()); } else { phaseSlice = fArray(parSlice1, chan_slice, Slice()); } ConcatArrays<casacore::Float>(selectedPhase, phaseSlice); } *pha_[chunk] = selectedPhase * Float(180.0/C::pi); } else { throw(AipsError("phase has no meaning for this table")); } break; } case PMS::GREAL: case PMS::REAL: { if (calType_ == "Fringe Jones") { // do not use float for this axis throw(AipsError("real has no meaning for this table")); } Cube<Complex> cArray; Cube<Float> fArray; bool isComplex(parsAreComplex()); if (isComplex) { cArray = cti.cparam(); } else { fArray = cti.fparam(); } Array<Float> selectedReal; for (auto& chan_slice : chansel) { Array<Float> realSlice; if (polnRatio_) { if (isComplex) { realSlice = real(cArray(parSlice1, chan_slice, Slice()) / cArray(parSlice2, chan_slice, Slice())); } else { // use float for single dish cal tables realSlice = fArray(parSlice1, chan_slice, Slice()) / fArray(parSlice2, chan_slice, Slice()); } checkRatioArray(realSlice, chunk); } else { if (isComplex) { realSlice = real(cArray(parSlice1, chan_slice, Slice())); } else { // use float for single dish cal tables realSlice = fArray(parSlice1, chan_slice, Slice()); } } ConcatArrays<casacore::Float>(selectedReal, realSlice); } *real_[chunk] = selectedReal; break; } case PMS::GIMAG: case PMS::IMAG: { if (parsAreComplex()) { Cube<Complex> cArray = cti.cparam(); Array<Float> selectedImag; for (auto& chan_slice : chansel) { Array<Float> imagSlice; if (polnRatio_) { imagSlice = imag(cArray(parSlice1, chan_slice, Slice()) / cArray(parSlice2, chan_slice, Slice())); checkRatioArray(imagSlice, chunk); } else { imagSlice = imag(cArray(parSlice1, chan_slice, Slice())); } ConcatArrays<casacore::Float>(selectedImag, imagSlice); } *imag_[chunk] = selectedImag; } else { throw(AipsError("imag has no meaning for this table")); } break; } case PMS::DELAY:{ if (!parsAreComplex()) { Cube<Float> fArray = cti.fparam(); Array<Float> selectedDelay; for (auto& chan_slice : chansel) { Array<Float> delaySlice; if (polnRatio_) { delaySlice = fArray(parSlice1, chan_slice, Slice()) - fArray(parSlice2, chan_slice, Slice()); } else { delaySlice = fArray(parSlice1, chan_slice, Slice()); } ConcatArrays<casacore::Float>(selectedDelay, delaySlice); } *par_[chunk] = selectedDelay; } else { throw(AipsError("delay has no meaning for this table")); } break; } case PMS::DELAY_RATE: { if (calType_.startsWith("Fringe") && !parsAreComplex()) { Cube<Float> fArray = cti.fparam(); Array<Float> selectedRate; for (auto& chan_slice : chansel) { Array<Float> rateSlice; if (polnRatio_) { rateSlice = fArray(parSlice1, chan_slice, Slice()) - fArray(parSlice2, chan_slice, Slice()); } else { rateSlice = fArray(parSlice1, chan_slice, Slice()); } ConcatArrays<casacore::Float>(selectedRate, rateSlice); } *par_[chunk] = selectedRate / 1.0e-12; } else { throw(AipsError("delay rate has no meaning for this table")); } break; } case PMS::DISP_DELAY: { if (calType_.startsWith("Fringe") && !parsAreComplex()) { Cube<Float> fArray = cti.fparam(); Array<Float> selectedDelay; for (auto& chan_slice : chansel) { Array<Float> delaySlice; if (polnRatio_) { delaySlice = fArray(parSlice1, chan_slice, Slice()) - fArray(parSlice2, chan_slice, Slice()); checkRatioArray(delaySlice, chunk); } else { delaySlice = fArray(parSlice1, chan_slice, Slice()); } ConcatArrays<casacore::Float>(selectedDelay, delaySlice); } // Divisor from PlotCal.cc *par_[chunk] = selectedDelay / 1.334537; } else { throw(AipsError("dispersive delay has no meaning for this table")); } break; } case PMS::OPAC: { if (!parsAreComplex() && calType_.contains("Opac")) { Cube<Float> fArray = cti.fparam(); Array<Float> selectedOpac; for (auto& chan_slice : chansel) { Array<Float> opacSlice = fArray(parSlice1, chan_slice, Slice()); ConcatArrays<casacore::Float>(selectedOpac, opacSlice); } *par_[chunk] = selectedOpac; } else { throw(AipsError( "opacity has no meaning for this table")); } break; } case PMS::SWP: { // "SPGAIN" in plotcal if ( !parsAreComplex() && calType_.contains("EVLASWPOW")) { Cube<Float> fArray = cti.fparam(); Array<Float> selectedSwPow; for (auto& chan_slice : chansel) { Array<Float> swpowSlice; if (polnRatio_) { swpowSlice = fArray(parSlice1, chan_slice, Slice()) / fArray(parSlice2, chan_slice, Slice()); checkRatioArray(swpowSlice, chunk); } else { swpowSlice = fArray(parSlice1, chan_slice, Slice()); } ConcatArrays<casacore::Float>(selectedSwPow, swpowSlice); } *par_[chunk] = selectedSwPow; } else { throw(AipsError("SwPower has no meaning for this table")); } break; } case PMS::TSYS: { if ((!parsAreComplex()) && (calType_.contains("EVLASWPOW") || calType_.contains("TSYS"))) { Cube<Float> fArray = cti.fparam(); Array<Float> selectedTsys; for (auto& chan_slice : chansel) { Array<Float> tsysSlice; if (polnRatio_) { tsysSlice = fArray(parSlice1, chan_slice, Slice()) / fArray(parSlice2, chan_slice, Slice()); checkRatioArray(tsysSlice, chunk); } else { tsysSlice = fArray(parSlice1, chan_slice, Slice()); } ConcatArrays<casacore::Float>(selectedTsys, tsysSlice); } *par_[chunk] = selectedTsys; } else { throw(AipsError( "Tsys has no meaning for this table")); } break; } case PMS::SNR: { Cube<Float> fArray = cti.snr(); Array<Float> selectedSnr; for (auto& chan_slice : chansel) { Array<Float> snrSlice; if (polnRatio_) { snrSlice = fArray(parSlice1, chan_slice, Slice()) / fArray(parSlice2, chan_slice, Slice()); checkRatioArray(snrSlice, chunk); } else { snrSlice = fArray(parSlice1, chan_slice, Slice()); } ConcatArrays<casacore::Float>(selectedSnr, snrSlice); } *snr_[chunk] = selectedSnr; break; } case PMS::TEC: { if (!parsAreComplex() && (calType_[0] == 'F') && (calType_ != "Fringe Jones")) { Cube<Float> fArray = cti.fparam(); Array<Float> selectedTec; for (auto& chan_slice : chansel) { // No correlation selection (e.g. poln ratio) on TEC axis Array<Float> tecSlice = fArray(parSlice1, chan_slice, Slice()); ConcatArrays<casacore::Float>(selectedTec, tecSlice); } *tec_[chunk] = selectedTec / 1e+16; } else throw(AipsError( "TEC has no meaning for this table")); break; } case PMS::FLAG: { Cube<Bool> flagArray = cti.flag(); Array<Bool> selectedFlag; for (auto& chan_slice : chansel) { Array<Bool> flagSlice; if (polnRatio_) { flagSlice = flagArray(parSlice1, chan_slice, Slice()) | flagArray(parSlice2, chan_slice, Slice()); } else { flagSlice = flagArray(parSlice1, chan_slice, Slice()); } ConcatArrays<casacore::Bool>(selectedFlag, flagSlice); } *flag_[chunk] = selectedFlag; break; } case PMS::WT: { *wt_[chunk] = cti.wt(); break; } case PMS::AZ0: case PMS::EL0: { MDirection azelMDir = cti.azel0(cti.thisTime()); Vector<Double> azel = azelMDir.getAngle("deg").getValue(); az0_(chunk) = azel(0); el0_(chunk) = azel(1); break; } case PMS::HA0: ha0_(chunk) = cti.hourang(cti.thisTime()) * 12/C::pi; // in hours break; case PMS::PA0: { pa0_(chunk) = cti.parang0(cti.thisTime()) * 180.0/C::pi; // in degrees if (pa0_(chunk) < 0.0) pa0_(chunk) += 360.0; break; } case PMS::ANTENNA: { antenna_[chunk]->resize(nAnt_); indgen(*antenna_[chunk]); break; } /* case PMS::AZIMUTH: case PMS::ELEVATION: { Matrix<Double> azel; cti.azelMat(cti.time()(0),azel); *az_[chunk] = azel.row(0); *el_[chunk] = azel.row(1); break; } case PMS::PARANG: *parang_[chunk] = cti.feed_pa(cti.time()(0))*(180.0/C::pi); //degrees break; case PMS::ROW: { *row_[chunk] = cti.rowIds(); break; } */ case PMS::OBSERVATION: { (*obsid_[chunk]).resize(1); *obsid_[chunk] = cti.thisObs(); break; } case PMS::INTENT: { // metadata axis that always gets loaded so don't want to throw exception break; } case PMS::ATM: case PMS::TSKY: case PMS::IMAGESB: { casacore::Int spw = cti.thisSpw(); casacore::Int scan = cti.thisScan(); casacore::Vector<casacore::Double> freqsGHz = cti.freq()/1e9; casacore::Vector<casacore::Double> curve(1, 0.0); if (axis == PMS::ATM) { if (plotmsAtm_) { plotmsAtm_->calcAtmTskyCurve(curve, spw, scan, freqsGHz); } } else if (axis == PMS::TSKY) { if (plotmsAtm_) { plotmsAtm_->calcAtmTskyCurve(curve, spw, scan, freqsGHz); } } else { if (plotmsAtm_) { plotmsAtm_->calcImageCurve(curve, spw, scan, freqsGHz); } } // Apply channel selection casacore::Vector<casacore::Double> selectedCurve; for (auto& chan_slice : chansel) { casacore::Vector<casacore::Double> curveSlice = curve(chan_slice); ConcatArrays<casacore::Double>(selectedCurve, curveSlice); } if (axis == PMS::ATM) { *atm_[chunk] = selectedCurve; } else if (axis == PMS::TSKY) { *tsky_[chunk] = selectedCurve; } else { *imageSideband_[chunk] = selectedCurve; } break; } default: throw(AipsError("Axis choice not supported for Cal Tables")); break; } } void CalCache::flagToDisk(const PlotMSFlagging& flagging, Vector<Int>& flchunks, Vector<Int>& flrelids, Bool flag, PlotMSIndexer* indexer, int dataIndex ) { // Sort the flags by chunk: Sort sorter; sorter.sortKey(flchunks.data(),TpInt); sorter.sortKey(flrelids.data(),TpInt); Vector<uInt> order; uInt nflag; nflag = sorter.sort(order,flchunks.nelements()); stringstream ss; // Make the VisIterator writable, with selection revised as appropriate NewCalTable* ct = new NewCalTable(filename_, Table::Update, Table::Plain); NewCalTable* selct = new NewCalTable(); selection_.apply(*ct, *selct); casacore::Bool readonly(False); // write access for flagging setUpCalIter(*selct, averaging(), readonly); ci_p->reset(); Int iflag(0); for (Int ichk=0;ichk<nChunk_;++ichk) { if (ichk!=flchunks(order[iflag]) && !ci_p->pastEnd()) // nothing to flag this chunk, just advance ci_p->next(); else { // This chunk requires flag-setting Int ifl(iflag); // Get bits we need from the table Cube<Bool> ctflag; Vector<Int> channel,a1,a2; ci_p->flag(ctflag); ci_p->chan(channel); ci_p->antenna1(a1); ci_p->antenna2(a2); // Apply poln selection Int npar; String pol = selection_.corr(); if (pol=="" || pol=="RL" || pol=="XY" || pol=="/") { // both axes npar = ctflag.shape()(0); } else { // poln selection using calParSlice String paramAxis = toVisCalAxis(PMS::FLAG); npar = getParSlice(paramAxis, pol).length(); } Int nchan = channel.nelements(); Int nrow = ci_p->nrow(); if (True) { Int currChunk=flchunks(order[iflag]); Double time=getTime(currChunk,0); Double cttime=ci_p->time()(0); Int spw=Int(getSpw(currChunk,0)); Int ctspw=ci_p->thisSpw(); Int field=Int(getField(currChunk,0)); Int ctfld=ci_p->thisField(); ss << "Time diff: " << time-cttime << " " << time << " " << cttime << "\n"; ss << "Spw diff: " << spw-ctspw << " " << spw << " " << ctspw << "\n"; ss << "Field diff: " << field-ctfld << " " << field << " " << ctfld << "\n"; } // Apply all flags in this chunk to this VB ifl=iflag; while (ifl<Int(nflag) && flchunks(order[ifl])==ichk) { Int currChunk=flchunks(order[ifl]); Int irel=flrelids(order[ifl]); Slice par1,chan,bsln; Slice par2 = Slice(); // Set flag range on par axis: if (netAxesMask_[dataIndex](0) && !flagging.corrAll()) { // A specific single par if (pol=="" || pol=="RL" || pol=="XY") { // flag both axes Int ipar=indexer->getIndex1000(currChunk,irel); par1 = Slice(ipar,1,1); } else if (polnRatio_) { par1 = getParSlice(toVisCalAxis(PMS::AMP), "R"); par2 = getParSlice(toVisCalAxis(PMS::AMP), "L"); } else { par1 = getParSlice(toVisCalAxis(PMS::AMP), pol); } } else { // all on par axis par1 = Slice(0,npar,1); } // Set Flag range on channel axis: if (netAxesMask_[dataIndex](1) && !flagging.channel()) { // A single specific channel Int ichan=indexer->getIndex0100(currChunk,irel); chan=Slice(ichan,1,1); } else { // Extend to all channels chan=Slice(0,nchan,1); } // Set Flags on the baseline axis: Int thisA1=Int(getAnt1(currChunk,indexer->getIndex0010(currChunk,irel))); Int thisA2=Int(getAnt2(currChunk,indexer->getIndex0010(currChunk,irel))); if (netAxesMask_[dataIndex](2) && !flagging.antennaBaselinesBased() && thisA1>-1 ) { // i.e., if baseline is an explicit data axis, // full baseline extension is OFF // and the first antenna in the selected point is > -1 // Do some variety of detailed per-baseline flagging for (Int irow=0;irow<nrow;++irow) { if (thisA2>-1) { // match a baseline exactly if (a1(irow)==thisA1 && a2(irow)==thisA2) { ctflag(par1,chan,Slice(irow,1,1)) = flag; if (par2.length() > 0) ctflag(par2,chan,Slice(irow,1,1)) = flag; break; // found the one baseline, escape from for loop } } else { // either antenna matches the one specified antenna // (don't break because there will be more than one) if (a1(irow)==thisA1 || a2(irow)==thisA1) { ctflag(par1,chan,Slice(irow,1,1)) = flag; if (par2.length() > 0) ctflag(par2,chan,Slice(irow,1,1)) = flag; } } } } else { // Set flags for all baselines, because the plot // is ordinarily implicit in baseline, we've turned on baseline // extension, or we've avaraged over all baselines bsln=Slice(0,nrow,1); ctflag(par1,chan,bsln) = flag; if (par2.length() > 0) ctflag(par2,chan,bsln) = flag; } ++ifl; } // while // Put the flags back into the MS wci_p->setflag(ctflag); // Advance to the next vb if (!ci_p->pastEnd()) ci_p->next(); else // we are done, so escape chunk loop break; // step over the flags we've just done iflag=ifl; // Escape if we are already finished if (uInt(iflag)>=nflag) break; } // flaggable chunk } // ichk // Delete the NCTs and VisIter so lock is released if (ct != nullptr) { delete ct; ct = nullptr; } if (selct != nullptr) { delete selct; selct = nullptr; } if (wci_p != nullptr) { delete wci_p; wci_p = nullptr; } ci_p = nullptr; logFlag(ss.str()); } // ======================== end NewCalTable ========================== // ======================== CalTable ========================== void CalCache::countChunks(Int nchunks, vector<PMS::Axis>& loadAxes, vector<PMS::DataColumn>& loadData, ThreadCommunication* thread) { // for CalTable if (thread!=nullptr) { thread->setStatus("Establishing cache size. Please wait..."); thread->setAllowedOperations(false,false,false); } setCache(nchunks, loadAxes, loadData); } void CalCache::getNamesFromMS() { // Set antenna and field names for Locate. MSMetaInfoForCal msmeta(msname_); msmeta.antennaNames(antnames_); antstanames_ = antnames_; msmeta.fieldNames(fldnames_); nAnt_ = msmeta.nAnt(); } void CalCache::setUpLoad(ThreadCommunication* thread, Slice& parSlice) { // common setup for CalTables // permit cancel in progress meter: if(thread != nullptr) thread->setAllowedOperations(false,false,true); logLoad("Loading chunks......"); // initial chunk info chshapes_.resize(4,nChunk_); goodChunk_.resize(nChunk_); goodChunk_.set(False); // get selected npol String polSelection(selection_.corr()); String paramAxis(toVisCalAxis(PMS::AMP)); // getParSlice() checks for valid axis and pol sel if (polnRatio_) { // just pick one for length parSlice = getParSlice(paramAxis, "R"); } else { parSlice = getParSlice(paramAxis, polSelection); } } void CalCache::getCalDataAxis(PMS::Axis axis, Cube<Complex>& viscube, Int chunk) { // Get axes derived from calculated data cube; // poln selection (parSlice) already applied to viscube switch(axis) { case PMS::GAMP: case PMS::AMP: { Cube<Float> ampcube = amplitude(viscube); if (polnRatio_) checkRatioArray(ampcube, chunk); *amp_[chunk] = ampcube; break; } case PMS::GPHASE: case PMS::PHASE: { Cube<Float> phasecube = phase(viscube); if (polnRatio_) checkRatioArray(phasecube, chunk); *pha_[chunk] = phasecube; (*pha_[chunk]) *= Float(180.0/C::pi); break; } case PMS::GREAL: case PMS::REAL: { Cube<Float> realcube = real(viscube); if (polnRatio_) checkRatioArray(realcube, chunk); *real_[chunk] = realcube; break; } case PMS::GIMAG: case PMS::IMAG: { Cube<Float> imagcube = imag(viscube); if (polnRatio_) checkRatioArray(imagcube, chunk); *imag_[chunk] = imagcube; break; } default: throw(AipsError("Axis choice not supported for Cal Tables")); break; } } // ======================== end CalTable ========================== // ======================== BPOLY ========================== void CalCache::loadBPoly(vector<PMS::Axis>& loadAxes, vector<PMS::DataColumn>& loadData, ThreadCommunication* thread) { // Load cache for a BPoly cal table BJonesPolyTable pt = BJonesPolyTable(filename_, Table::Update); // Set ms-related data ROCalDescColumns calDescCol(pt); String msname(calDescCol.msName()(0)); setMSName(msname); // add path if (msname_.empty() || !Table::isReadable(msname_)) { throw(AipsError("Associated MS is not available, cannot plot solutions.")); } getNamesFromMS(); // field and antenna // Set number of antennas nAnt_ = pt.maxAntenna() + 1; // Create a B Jones table derived from BPOLY table NewCalTable* nct = virtualBPoly(pt); parsAreComplex_ = nct->isComplex(); // Apply selection to get selected cal table NewCalTable* selct = new NewCalTable(); Vector<Vector<Slice> > chansel; Vector<Vector<Slice> > corrsel; selection_.apply(*nct, *selct); // Use NewCalTable implementation casacore::Bool readonly(True); setUpCalIter(*selct, averaging(), readonly); countChunks(*ci_p, loadAxes, loadData, thread); loadCalChunks(*ci_p, loadAxes, thread); } NewCalTable* CalCache::virtualBPoly(BJonesPolyTable& polyTable) { // Returns B Jones NewCalTable derived from BPOLY; based on BPoly::loadMemCalTable. // Generate a NCT in memory to hold the BPOLY as a B NewCalTable* nct = new NewCalTable("BPolyAsB.tmp", VisCalEnum::COMPLEX, "B Jones", msname_, false); // Frequency info from MS Vector<Vector<Double>> mschanfreq; // [nspw, nchan] getChanFreqsFromMS(mschanfreq); // Ensure sort on TIME, so CalSet is filled in order Block <String> sortCol(3); sortCol[0]="CAL_DESC_ID"; sortCol[1]="TIME"; sortCol[2]="ANTENNA1"; polyTable.sort2(sortCol); Int nrows = polyTable.nRowMain(); Int nDesc = polyTable.nRowDesc(); // Spws to be calibrated by each caldesc Vector<Int> spwmap(nDesc,-1); for (Int idesc = 0; idesc < nDesc; ++idesc) { CalDescRecord calDescRec(polyTable.getRowDesc(idesc)); Vector<Int> spwId; calDescRec.getSpwId(spwId); Int currSpw = spwId(0); spwmap(idesc) = currSpw; // Set SPW subtable freqs Vector<Double> freq = mschanfreq(currSpw); nct->setSpwFreqs(currSpw, freq); } // Solve arrays, so we can fill them Cube<Complex> cparam; Cube<Bool> flag; Cube<Float> err, snr; bool cubeInit(false); // Attach a calibration table columns accessor BJonesPolyMCol maincol(polyTable); for (Int row = 0; row < nrows; row++) { // Extract the polynomial coefficients in amplitude and phase Int nAmp = maincol.nPolyAmp().asInt(row); Int nPhase = maincol.nPolyPhase().asInt(row); Array<Double> ampCoeffArray, phaseCoeffArray; maincol.polyCoeffAmp().get(row, ampCoeffArray); maincol.polyCoeffPhase().get(row, phaseCoeffArray); Matrix<Double> ampCoeff(nAmp, 2); IPosition ampPos = ampCoeffArray.shape(); ampPos = 0; for (Int k = 0; k < 2 * nAmp; k++) { ampPos.setLast(IPosition(1, k)); ampCoeff(k % nAmp, k / nAmp) = ampCoeffArray(ampPos); }; Matrix<Double> phaseCoeff(nPhase, 2); IPosition phasePos = phaseCoeffArray.shape(); phasePos = 0; for (Int k = 0; k < 2 * nPhase; k++) { phasePos.setLast(IPosition(1, k)); phaseCoeff(k % nPhase, k / nPhase) = phaseCoeffArray(phasePos); }; // Get frequencies for this spw Int nPol(2); Int thisDesc = maincol.calDescId().asInt(row); Int thisSpw = spwmap(thisDesc); Vector<Double> freq = mschanfreq(thisSpw); Int nChan = freq.nelements(); // Extract the valid domain for the polynomial Vector<Double> freqDomain(2); maincol.validDomain().get(row, freqDomain); Double x1 = freqDomain(0); Double x2 = freqDomain(1); Complex factor = maincol.scaleFactor().asComplex(row); Int thisAnt1 = maincol.antenna1().asInt(row); // Resize and initialize solve arrays if (!cubeInit) { cparam.resize(nPol, nChan, nAnt_); cparam.set(Complex(1.0)); flag.resize(nPol, nChan, nAnt_); flag.set(true); err.resize(nPol, nChan, nAnt_); err.set(0.0); snr.resize(nPol, nChan, nAnt_); snr.set(1.0); cubeInit = true; } for (Int pol = 0; pol < 2; pol++) { Vector<Double> ac(ampCoeff.column(pol)); Vector<Double> pc(phaseCoeff.column(pol)); // Only do calculation if coeffs are non-zero if (anyNE(ac, Double(0.0)) || anyNE(pc, Double(0.0)) ) { for (Int chan = 0; chan < nChan; ++chan) { Double ampval(1.0), phaseval(0.0); // Calculate Cheby if freq in domain Double thisFreq(freq(chan)); if ((thisFreq >= x1) && (thisFreq <= x2)) { ampval = getChebVal(ac, x1, x2, thisFreq); phaseval = getChebVal(pc, x1, x2, thisFreq); cparam(pol, chan, thisAnt1) = factor * Complex(exp(ampval)) * Complex(cos(phaseval), sin(phaseval)); flag(pol, chan, thisAnt1) = false; } else { // Unflagged unit calibration for now cparam(pol, chan, thisAnt1) = Complex(1.0); flag(pol, chan, thisAnt1) = false; } } // chan } } // pol // Every nAnt rows, store the result if ((row + 1) % nAnt_ == 0) { Double thisTime = maincol.time().asdouble(row); Double thisInterval = maincol.interval().asdouble(row); Int thisField = maincol.fieldId().asInt(row); Int thisObs = maincol.obsId().asInt(row); Array<Int> refant = maincol.refAnt().get(row); IPosition first(refant.shape().size(), 0); Int thisRefant = refant(first); Vector<Int> ant1list; // NewCalTable generates antenna ids based on nrows nct->fillAntBasedMainRows(nAnt_, thisTime, thisInterval, thisField, thisSpw, thisObs, ant1list, thisRefant, cparam, flag, err, snr); // reset arrays next loop cubeInit = false; } } // rows return nct; } void CalCache::getChanFreqsFromMS(Vector< Vector<Double> >& mschanfreqs) { // shape is (nspw, nchan) MeasurementSet ms(msname_); MSColumns mscol(ms); uInt nspw = mscol.spectralWindow().nrow(); mschanfreqs.resize(nspw); for (uInt spw=0; spw<nspw; ++spw) { mschanfreqs(spw) = mscol.spectralWindow().chanFreq().get(spw); } } Double CalCache::getChebVal(const Vector<Double>& coeff, const Double& xinit, const Double& xfinal, const Double& x) { // from synthesis/MeasurementComponents/BPoly.cc // Compute a Chebyshev polynomial value using the CLIC library // Input: // coeff const Vector<Double>& Chebyshev coefficients // xinit const Double& Domain start // xfinal const Double& Domain end // x const Double& x-ordinate // Output: // getChebVal Double Chebyshev polynomial value // // Re-scale x-ordinate Double xcap = ((x - xinit) - (xfinal - x)) / (xfinal - xinit); // Compute polynomial Int deg = coeff.shape().asVector()(0); Vector<Double> val(deg); Bool check; Int checkval; cheb(°, &xcap, val.getStorage(check), &checkval); Double soly(0.0); for (Int mm = 0; mm < deg; mm++){ soly += coeff[mm] * val[mm]; } return soly; } // ======================== end BPOLY ========================== // ======================== GSPLINE ========================== void CalCache::loadGSpline(vector<PMS::Axis>& loadAxes, vector<PMS::DataColumn>& loadData, ThreadCommunication* thread) { GJonesSplineTable ct = GJonesSplineTable(filename_); GJonesSplineTable selct(ct); if (!selection_.timerange().empty()) { logWarn("PlotMS::load_cache", "Time selection not supported for GSPLINE calibration tables"); selection_.setTimerange(""); } // chansel not applicable, corrsel done with parSlice Vector<Vector<Slice> > chansel, corrsel; selection_.apply(ct, selct, chansel, corrsel); Vector<Int> selectedAnts = selection_.getSelectedAntennas1(); ROGJonesSplineMCol mainCol(selct); ROCalDescColumns calDescCol(selct); String msname(calDescCol.msName()(0)); setMSName(msname); // add path if (msname_.empty() || !Table::isReadable(msname_)) { throw(AipsError("Associated MS is not available, cannot plot solutions.")); } getNamesFromMS(); // field and antenna // count and load chunks Int nsample(1000); // make time samples to load cache countChunks(nsample, loadAxes, loadData, thread); loadCalChunks(mainCol, calDescCol, nsample, loadAxes, selectedAnts, thread); } void CalCache::loadCalChunks(ROGJonesSplineMCol& mcol, ROCalDescColumns& dcol, Int nsample, const vector<PMS::Axis> loadAxes, Vector<int>& selectedAnts, ThreadCommunication* thread) { // GSPLINE does not load per chunk or row as in other cal tables. // Its "chunks" are per-time sample, generated from SPLINE_KNOTS_AMP/PHASE. // Therefore: // Scalar columns are read once and all values plotted per time sample. // Columns that normally have one value per chunk (field, spw, scan, time, // and interval) will use the value in the column if all are the same; // else multiple values result in setting the value to -1 (for locate). // The exception is field, which plots the field used for solutions. Slice parslice, parslice2; setUpLoad(thread, parslice); Int nPol = parslice.length(); // for chunk shapes if (polnRatio_) { parslice2 = getParSlice(toVisCalAxis(PMS::AMP), "L"); } // load main table metadata once; use vector/value for every timestamp // field Vector<Int> fields(mcol.fieldId().getColumn()); Int nItems = GenSort<Int>::sort(fields, Sort::Ascending, Sort::NoDuplicates); Int firstField(fields(0)); Int fieldForSolve = (firstField < 0 ? 0 : firstField); Int fieldForCache = (nItems>1 ? -1 : firstField); // spw Array<Int> spws(dcol.spwId().getColumn()); nItems = GenSort<Int>::sort(spws, Sort::Ascending, Sort::NoDuplicates); Int firstSpw(spws(IPosition(2,0,0))); Int spwForSolve = (firstSpw < 0 ? 0 : firstSpw); Int spwForCache = (nItems>1 ? -1 : firstSpw); // scan Vector<Int> scans(mcol.scanNo().getColumn()); nItems = GenSort<Int>::sort(scans, Sort::Ascending, Sort::NoDuplicates); Int firstScan(scans(0)); Int scanForSolve(firstScan < 0 ? 0 : firstScan); Int scanForCache = (nItems>1 ? -1 : firstScan); // antenna1 Vector<Int> ant1(selectedAnts); // use selected antenna1 Int nSelAnts = ant1.size(); Int nAnt(nSelAnts); if (nAnt == 0) { // no ant1 selection, get from main table ant1 = mcol.antenna1().getColumn(); nItems = GenSort<Int>::sort(ant1, Sort::Ascending, Sort::NoDuplicates); ant1.resize(nItems, true); nAnt = nItems; } // obsid Vector<Int> obsid(mcol.obsId().getColumn()); nItems = GenSort<Int>::sort(obsid, Sort::Ascending, Sort::NoDuplicates); obsid.resize(nItems, true); Int obsForSolve(obsid(0) < 0 ? 0 : obsid(0)); // feed1 Vector<Int> feed1(mcol.feed1().getColumn()); nItems = GenSort<Int>::sort(feed1, Sort::Ascending, Sort::NoDuplicates); feed1.resize(nItems, true); // interval Vector<Double> intervals(mcol.interval().getColumn()); nItems = GenSort<Double>::sort(intervals, Sort::Ascending, Sort::NoDuplicates); intervals.resize(nItems, true); Double interval = (nItems>1 ? -1 : intervals(0)); // poln Vector<Int> polns(parslice.length()); if (polnRatio_) { polns = -1; } else { Int start(parslice.start()), inc(parslice.inc()); indgen(polns, start, inc); } // set up gspline MSMetaInfoForCal msmeta(msname_); GJonesSpline* gspline = new GJonesSpline(msmeta); String mode(mcol.polyMode()(0)); Record rec; // for solving params rec.define("table", filename_); rec.define("apmode", mode); gspline->setApply(rec); // set up calbuffer for calcPar gspline->setSolve(rec); // set mode // get first row for chosen field uInt row; Vector<Int> fieldcol(mcol.fieldId().getColumn()); for (row=0; row<fieldcol.size(); ++row) if (fieldcol(row) == fieldForSolve) break; MFrequency refFreq = mcol.refFreqMeas()(row)(IPosition(3,0,0,0)); Double refFreqHz = refFreq.get("Hz").getValue(); Vector<Double> freq(1, refFreqHz); // for setMeta // Create 1000 time samples from spline knots: see PlotCal::virtualGSpline Vector<Double> splineKnots, times(nsample); // get splineKnots for that row if (mode.contains("AMP") || mode.contains("A&P")) { splineKnots = mcol.splineKnotsAmp()(row); } else if (mode.contains("PHAS") || mode.contains("A&P")) { splineKnots = mcol.splineKnotsPhase()(row); } // make samples based on spline knots Double dt((max(splineKnots)-min(splineKnots)) / Double(nsample)); Double mintime(splineKnots(0) + (dt/2.0)); for (Int sample=0; sample<nsample; sample++) { times(sample) = mintime + sample*dt; } // Now ready to load cache for (Int sample=0; sample<nsample; sample++) { // set time and field for calcPar Double time = times(sample); gspline->setMeta(obsForSolve, scanForSolve, time, spwForSolve, freq, fieldForSolve); gspline->doCalcPar(); // Cache the data shapes chshapes_(0,sample) = nPol; chshapes_(1,sample) = 1; // nChan chshapes_(2,sample) = nAnt; chshapes_(3,sample) = 1; goodChunk_(sample) = True; // load axes for each row for(unsigned int i = 0; i < loadAxes.size(); i++) { PMS::Axis axis = loadAxes[i]; // slice viscube for poln selection Slicer slicer1(Slicer(parslice, Slice(), Slice())); if (PMS::axisIsData(axis)) { // amp, phase, real, imag Cube<Complex> cpar, viscube; cpar = gspline->currCPar(); viscube = cpar(slicer1); if (polnRatio_) { Slicer slicer2(Slicer(parslice2, Slice(), Slice())); viscube /= cpar(slicer2); } if (nAnt == nSelAnts) // get selected rows getSelectedCube(viscube, selectedAnts); getCalDataAxis(axis, viscube, sample); } else { switch(axis) { case PMS::SCAN: scan_(sample) = scanForCache; break; case PMS::FIELD: field_(sample) = fieldForCache; break; case PMS::TIME: time_(sample) = time; break; case PMS::TIME_INTERVAL: timeIntr_[sample] = interval; break; case PMS::SPW: spw_(sample) = spwForCache; break; case PMS::CORR: { corr_[sample]->resize(nPol); String pol = selection_.corr(); if (pol=="" || pol=="RL" || pol=="XY") { indgen(*corr_[sample]); } else if (pol=="/") { corr_[sample]->resize(1); corr_[sample]->set(-1); } else { // R/X or L/Y Int poln = ((pol=="R" || pol=="X") ? 0 : 1); corr_[sample]->resize(1); corr_[sample]->set(poln); } break; } case PMS::ANTENNA1: *antenna1_[sample] = ant1; break; case PMS::ANTENNA: *antenna_[sample] = ant1; break; case PMS::FLAG: { Cube<Bool> parOK = gspline->currParOK(); // OK=true means flag=false Cube<Bool> flagcube(!parOK(slicer1)); if (polnRatio_) { Slicer slicer2(Slicer(parslice2, Slice(), Slice())); flagcube &= !parOK(slicer2); } if (nAnt == nSelAnts) // get selected rows getSelectedCube(flagcube, selectedAnts); *flag_[sample] = flagcube; break; } case PMS::ROW: { Vector<rownr_t> sampleRow(nAnt, sample); *row_[sample] = sampleRow; break; } case PMS::OBSERVATION: *obsid_[sample] = obsid; break; case PMS::FEED1: *feed1_[sample] = feed1; break; default: // invalid axes weeded out in checkAxes break; } } } // If a thread is given, update it. if(thread != nullptr) { double progress = ((double)sample) / nsample; thread->setProgress((unsigned int)((progress * 100) + 0.5)); } } if (divZero_) logWarn("CalCache::loadIt", "Caught divide-by-zero exception in ratio plots; result(s) set to 1.0 and flagged"); } void CalCache::checkAxes(const vector<PMS::Axis>& loadAxes) { // trap user-requested axes that are invalid for GSPLINE for(unsigned int i = 0; i < loadAxes.size(); i++) { PMS::Axis axis(loadAxes[i]); switch (axis) { case PMS::SCAN: case PMS::FIELD: case PMS::TIME: case PMS::SPW: case PMS::TIME_INTERVAL: case PMS::CORR: case PMS::ANTENNA1: case PMS::ANTENNA: // same as antenna1 case PMS::AMP: case PMS::GAMP: case PMS::PHASE: case PMS::GPHASE: case PMS::REAL: case PMS::GREAL: case PMS::IMAG: case PMS::GIMAG: case PMS::FLAG: case PMS::OBSERVATION: case PMS::FEED1: { // allowed break; } case PMS::CHANNEL: case PMS::FREQUENCY: case PMS::SNR: { String msg("GSPLINE plotting does not support " + PMS::axis(axis) + " axis"); if (axis==PMS::ANTENNA) msg += "; use ANTENNA1"; throw(AipsError(msg)); break; } case PMS::ANTENNA2: case PMS::BASELINE: case PMS::ROW: case PMS::DELAY: case PMS::DELAY_RATE: case PMS::DISP_DELAY: case PMS::OPAC: case PMS::SWP: case PMS::TSYS: case PMS::TEC: case PMS::INTENT: { String msg(PMS::axis(axis) + " has no meaning for this table"); throw(AipsError(msg)); break; } default: throw(AipsError("Axis choice not supported for Cal Tables")); break; } } } template<class T> void CalCache::getSelectedCube(Cube<T>& inputCube, const Vector<Int>& selectedRows) { // replaces input cube with cube selected by rows in vector Cube<T> selectedCube; for (uInt irow=0; irow<selectedRows.size(); ++irow) { Slice rowSlice = Slice(selectedRows(irow)); Slicer rowSlicer = Slicer(Slice(), Slice(), rowSlice); Cube<T> concatCube = concatenateArray(selectedCube, inputCube(rowSlicer)); selectedCube.resize(); selectedCube = concatCube; } inputCube.resize(); inputCube = selectedCube; } // ======================== end GSPLINE ========================== String CalCache::toVisCalAxis(PMS::Axis axis) { switch (axis) { // FLAG and SNR have same shape as AMP // and should be sliced the same way case PMS::AMP: case PMS::GAMP: case PMS::FLAG: case PMS::SNR: if (calType_.contains("EVLASWP")) return "GAINAMP"; if (calType_.contains("TSYS")) return "TSYS"; if (calType_[0] == 'K' && !calType_.startsWith("KAntPos")) return "DELAY"; if (calType_ == "F Jones") return "TEC"; if (calType_ == "Fringe Jones") return "DELAY"; if (calType_ == "TOpac") return "OPAC"; return "AMP"; break; case PMS::PHASE: case PMS::GPHASE: return "PHASE"; break; case PMS::REAL: case PMS::GREAL: return "REAL"; break; case PMS::IMAG: case PMS::GIMAG: return "IMAG"; break; case PMS::DELAY_RATE: return "RATE"; break; case PMS::DISP_DELAY: return "DISP"; break; default: return PMS::axis(axis); break; } } Slice CalCache::getParSlice(String axis, String polnSel) { Slice parSlice = Slice(); try { parSlice = viscal::calParSlice(filename_, axis, polnSel); } catch(AipsError& err) { if (err.getMesg().contains("Unsupported value type")) { // Message a bit vague at top level, add some explanation String errMsg = err.getMesg() + ". Invalid axis or polarization selection for cal table type."; throw (AipsError(errMsg)); } else if (calType_ == "M Mueller") { if (polnSel.empty()) { return Slice(0, 2, 1); // full selection } // include default "RL" when not set by user std::vector<casacore::String> valid_poln = {"R", "X", "RL"}; for (auto& poln : valid_poln) { if (poln == polnSel) { return Slice(0, 2, 1); // full selection } } throw (AipsError("Invalid polarization selection for cal table type.")); } else { // unsupported cal type throw(AipsError(err)); } } return parSlice; } void CalCache::checkRatioArray(Array<Float>& array, Int chunk) { // When array is result of division, check for division by zero (element is infinite) // Reset value to 1.0, flag, and set divZero_ for warning Cube<Float> ratioCube; ratioCube.reference(array); Cube<Bool> flags; flags.reference(*flag_[chunk]); IPosition cubeShape = ratioCube.shape(); for (uInt i=0; i<cubeShape[0]; ++i) { for (uInt j=0; j<cubeShape[1]; ++j) { for (uInt k=0; k<cubeShape[2]; ++k) { if (isInf(ratioCube(i,j,k))) { ratioCube(i,j,k) = 1.0; flags(i,j,k) = True; divZero_ = True; } } } } } } // namespace casa