Line data Source code
1 : //# Calibrater.cc: Implementation of Calibrater.h
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program 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 General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18 : //#
19 : //# Correspondence concerning AIPS++ should be addressed 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 : //# $Id: Calibrater.cc,v 19.37 2006/03/16 01:28:09 gmoellen Exp $
27 :
28 : #include <casacore/tables/Tables/Table.h>
29 : #include <casacore/tables/Tables/TableRecord.h>
30 : #include <casacore/tables/Tables/TableDesc.h>
31 : #include <casacore/tables/Tables/TableLock.h>
32 : #include <casacore/tables/Tables/TableUtil.h>
33 : #include <casacore/tables/TaQL/TableParse.h>
34 : #include <casacore/tables/Tables/ArrColDesc.h>
35 : #include <casacore/tables/DataMan/TiledShapeStMan.h>
36 :
37 : #include <casacore/casa/System/AipsrcValue.h>
38 : #include <casacore/casa/Arrays/ArrayUtil.h>
39 : #include <casacore/casa/Arrays/ArrayLogical.h>
40 : #include <casacore/casa/Arrays/ArrayPartMath.h>
41 : //#include <casa/Arrays/ArrayMath.h>
42 : #include <casacore/ms/MeasurementSets/MSColumns.h>
43 : #include <casacore/ms/MSSel/MSFieldIndex.h>
44 : #include <casacore/ms/MSSel/MSSelection.h>
45 : #include <casacore/ms/MSSel/MSSelectionTools.h>
46 : #include <casacore/ms/MSOper/MSMetaData.h>
47 : #include <casacore/casa/BasicSL/Constants.h>
48 : #include <casacore/casa/Exceptions/Error.h>
49 : #include <iostream>
50 : #include <sstream>
51 : #include <casacore/casa/OS/File.h>
52 : #include <synthesis/MeasurementComponents/Calibrater.h>
53 : #include <synthesis/CalTables/CLPatchPanel.h>
54 : #include <synthesis/MeasurementComponents/CalSolVi2Organizer.h>
55 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
56 : #include <synthesis/MeasurementComponents/VisCalSolver.h>
57 : #include <synthesis/MeasurementComponents/VisCalSolver2.h>
58 : #include <synthesis/MeasurementComponents/UVMod.h>
59 : #include <synthesis/MeasurementComponents/TsysGainCal.h>
60 : #include <msvis/MSVis/VisSetUtil.h>
61 : #include <msvis/MSVis/VisBuffAccumulator.h>
62 : #include <msvis/MSVis/VisibilityIterator2.h>
63 : #include <msvis/MSVis/VisBuffer2.h>
64 : #include <msvis/MSVis/ViFrequencySelection.h>
65 : #include <casacore/casa/Quanta/MVTime.h>
66 : #include <casacore/casa/Logging/LogMessage.h>
67 : #include <casacore/casa/Logging/LogIO.h>
68 : #include <casacore/casa/Utilities/Assert.h>
69 :
70 : #include <casacore/tables/Tables/SetupNewTab.h>
71 : #include <vector>
72 : using std::vector;
73 : #include <stdcasa/UtilJ.h>
74 :
75 : #ifdef _OPENMP
76 : #include <omp.h>
77 : #endif
78 :
79 : //#define REPORT_CAL_TIMING
80 :
81 : using namespace casacore;
82 : using namespace casa::utilj;
83 :
84 : using namespace casacore;
85 : using namespace casa::vpf;
86 :
87 : using namespace casacore;
88 : namespace casa { //# NAMESPACE CASA - BEGIN
89 :
90 31 : CalCounts::CalCounts():
91 31 : antennaMap_(),
92 31 : spwMap_(),
93 62 : totalMap_()
94 31 : {}
95 :
96 31 : void CalCounts::initCounts(Int NSpw, Int NAnt, Int NPol) {
97 31 : Vector<Int> polShape(NPol, 0);
98 31 : Vector<Int> singleVector(1, 0);
99 :
100 31 : totalMap_["expected"] = polShape;
101 31 : totalMap_["data_unflagged"] = polShape;
102 31 : totalMap_["above_minblperant"] = polShape;
103 31 : totalMap_["above_minsnr"] = polShape;
104 :
105 155 : for (Int spw=0; spw<NSpw; spw++) {
106 : // Init spw map keys
107 124 : spwMap_[spw]["expected"] = polShape;
108 124 : spwMap_[spw]["data_unflagged"] = polShape;
109 124 : spwMap_[spw]["above_minblperant"] = polShape;
110 124 : spwMap_[spw]["above_minsnr"] = polShape;
111 1388 : for (Int ant=0; ant<NAnt; ant++) {
112 : // Init antenna map keys
113 1264 : antennaMap_[spw][ant]["expected"] = polShape;
114 1264 : antennaMap_[spw][ant]["data_unflagged"] = polShape;
115 1264 : antennaMap_[spw][ant]["above_minblperant"] = polShape;
116 1264 : antennaMap_[spw][ant]["above_minsnr"] = polShape;
117 1264 : antennaMap_[spw][ant]["used_as_refant"] = singleVector;
118 : }
119 : }
120 31 : }
121 :
122 225 : void CalCounts::addAntennaCounts(Int spw, Int NAnt, Int NPol, std::map<Int, std::map<String, Vector<Int>>> resultMap) {
123 : // Add the results to the antenna Map
124 225 : Vector<Int> spwExp(NPol, 0), spwUnflag(NPol, 0), spwMinsnr(NPol, 0);
125 :
126 2457 : for (Int ant=0; ant<NAnt; ant++) {
127 2232 : if (resultMap.find(ant) == resultMap.end()) {continue;};
128 : // add if used as refant
129 2232 : antennaMap_[spw][ant]["used_as_refant"][0] += resultMap[ant]["used_as_refant"][0];
130 20088 : for (Int pol=0; pol<NPol; pol++) {
131 17856 : antennaMap_[spw][ant]["expected"][pol] += resultMap[ant]["expected"][pol];
132 17856 : antennaMap_[spw][ant]["data_unflagged"][pol] += resultMap[ant]["data_unflagged"][pol];
133 17856 : antennaMap_[spw][ant]["above_minblperant"][pol] += resultMap[ant]["above_minblperant"][pol];
134 17856 : antennaMap_[spw][ant]["above_minsnr"][pol] += resultMap[ant]["above_minsnr"][pol];
135 : // if any antenna value is good add to spw total
136 17856 : if (resultMap[ant]["expected"][pol] == 1) {spwExp[pol] = 1;};
137 17856 : if (resultMap[ant]["data_unflagged"][pol] == 1) {spwUnflag[pol] = 1;};
138 17856 : if (resultMap[ant]["above_minblperant"][pol] == 1) {spwMinsnr[pol] = 1;};
139 17856 : if (resultMap[ant]["above_minsnr"][pol] == 1) {spwMinsnr[pol] = 1;};
140 : }
141 : }
142 : // Add the totals to the spw Map
143 2025 : for (Int pol=0; pol<NPol; pol++) {
144 1800 : spwMap_[spw]["expected"][pol] += spwExp[pol];
145 1800 : spwMap_[spw]["data_unflagged"][pol] += spwUnflag[pol];
146 1800 : spwMap_[spw]["above_minblperant"][pol] += spwMinsnr[pol];
147 1800 : spwMap_[spw]["above_minsnr"][pol] += spwMinsnr[pol];
148 : // Acumulate total overall values
149 1800 : totalMap_["expected"][pol] += spwExp[pol];
150 1800 : totalMap_["data_unflagged"][pol] += spwUnflag[pol];
151 1800 : totalMap_["above_minblperant"][pol] += spwMinsnr[pol];
152 1800 : totalMap_["above_minsnr"][pol] += spwMinsnr[pol];
153 : }
154 :
155 225 : }
156 :
157 31 : void CalCounts::updateRefants(Int NSpw, Int NAnt, std::map<Int, std::map<Int, Int>> refantMap) {
158 155 : for (Int spw=0; spw<NSpw; spw++) {
159 1388 : for (Int ant=0; ant<NAnt; ant++) {
160 1264 : antennaMap_[spw][ant]["used_as_refant"][0] += refantMap[spw][ant];
161 : }
162 : }
163 31 : }
164 :
165 1264 : Vector<Int> CalCounts::antMapVal(Int spw, Int ant, String gate) {
166 1264 : return antennaMap_[spw][ant][gate];
167 : }
168 :
169 496 : Vector<Int> CalCounts::spwMapVal(Int spw, String gate) {
170 496 : return spwMap_[spw][gate];
171 : }
172 :
173 124 : Vector<Int> CalCounts::totalMapVal(String gate) {
174 124 : return totalMap_[gate];
175 : }
176 :
177 31 : Record CalCounts::makeRecord(Int NAnt, Int NPol) {
178 :
179 31 : Vector<Int> totExp(NPol, 0), totUnflag(NPol, 0), totMinsnr(NPol, 0);
180 :
181 31 : nSpw = spwMap_.size();
182 31 : Record containerRec = Record();
183 31 : Record resultRec = Record();
184 :
185 31 : resultRec.define("expected", totalMap_["expected"]);
186 31 : resultRec.define("data_unflagged", totalMap_["data_unflagged"]);
187 31 : resultRec.define("above_minblperant", totalMap_["above_minblperant"]);
188 31 : resultRec.define("above_minsnr", totalMap_["above_minsnr"]);
189 :
190 155 : for (Int spw=0; spw<nSpw; spw++) {
191 124 : Record spwRec = Record();
192 124 : spwRec.define("expected", spwMap_[spw]["expected"]);
193 124 : spwRec.define("data_unflagged", spwMap_[spw]["data_unflagged"]);
194 124 : spwRec.define("above_minblperant", spwMap_[spw]["above_minblperant"]);
195 124 : spwRec.define("above_minsnr", spwMap_[spw]["above_minsnr"]);
196 1388 : for (Int ant=0; ant<NAnt; ant++) {
197 1264 : Record subRec = Record();
198 1264 : subRec.define("expected", antennaMap_[spw][ant]["expected"]);
199 1264 : subRec.define("data_unflagged", antennaMap_[spw][ant]["data_unflagged"]);
200 1264 : subRec.define("above_minblperant", antennaMap_[spw][ant]["above_minblperant"]);
201 1264 : subRec.define("above_minsnr", antennaMap_[spw][ant]["above_minsnr"]);
202 1264 : subRec.define("used_as_refant", antennaMap_[spw][ant]["used_as_refant"]);
203 1264 : spwRec.defineRecord(RecordFieldId("ant"+to_string(ant)), subRec);
204 1264 : }
205 124 : resultRec.defineRecord(RecordFieldId("spw"+to_string(spw)), spwRec);
206 124 : }
207 31 : containerRec.defineRecord(RecordFieldId("solvestats"), resultRec);
208 :
209 62 : return containerRec;
210 31 : }
211 :
212 0 : CalCounts::~CalCounts() {}
213 :
214 242 : Calibrater::Calibrater():
215 242 : ms_p(0),
216 242 : mssel_p(0),
217 242 : mss_p(0),
218 242 : frequencySelections_p(nullptr),
219 242 : msmc_p(0),
220 242 : ve_p(0),
221 242 : vc_p(),
222 242 : svc_p(0),
223 242 : histLockCounter_p(),
224 242 : hist_p(0),
225 242 : usingCalLibrary_(false),
226 242 : corrDepFlags_(false), // default (==traditional behavior)
227 242 : actRec_(),
228 242 : resRec_(),
229 242 : simdata_p(false),
230 726 : ssvp_p()
231 : {
232 : // cout << "This is the NEW VI2-aware Calibrater" << endl;
233 242 : }
234 :
235 0 : Calibrater::Calibrater(String msname):
236 0 : msname_p(msname),
237 0 : ms_p(0),
238 0 : mssel_p(0),
239 0 : mss_p(0),
240 0 : frequencySelections_p(nullptr),
241 0 : msmc_p(0),
242 0 : ve_p(0),
243 0 : vc_p(),
244 0 : svc_p(0),
245 0 : histLockCounter_p(),
246 0 : hist_p(0),
247 0 : usingCalLibrary_(false),
248 0 : corrDepFlags_(false), // default (==traditional behavior)
249 0 : actRec_(),
250 0 : resRec_(),
251 0 : simdata_p(false),
252 0 : ssvp_p()
253 : {
254 :
255 :
256 0 : if (!Table::isReadable(msname))
257 0 : throw(AipsError("MS "+msname+" does not exist."));
258 :
259 :
260 0 : logSink() << LogOrigin("Calibrater","") << LogIO::NORMAL
261 0 : << "Arranging to calibrate MS: "+msname
262 0 : << LogIO::POST;
263 :
264 : // This is a bare Calibrater, intended to serve a VisEquation
265 :
266 : // We need very little of the usual stuff
267 :
268 : // A VisEquation
269 0 : ve_p = new VisEquation();
270 :
271 : // Reset the apply/solve VisCals
272 0 : reset();
273 :
274 0 : }
275 :
276 0 : Calibrater::Calibrater(const vi::SimpleSimVi2Parameters& ssvp):
277 0 : msname_p("<noms>"),
278 0 : ms_p(0),
279 0 : mssel_p(0),
280 0 : mss_p(0),
281 0 : frequencySelections_p(nullptr),
282 0 : msmc_p(0),
283 0 : ve_p(0),
284 0 : vc_p(),
285 0 : svc_p(0),
286 0 : histLockCounter_p(),
287 0 : hist_p(0),
288 0 : usingCalLibrary_(false),
289 0 : corrDepFlags_(false), // default (==traditional behavior)
290 0 : actRec_(),
291 0 : resRec_(),
292 0 : simdata_p(true),
293 0 : ssvp_p(ssvp)
294 : {
295 :
296 0 : logSink() << LogOrigin("Calibrater","") << LogIO::NORMAL
297 : << "Arranging SIMULATED MS data for testing!!!!"
298 0 : << LogIO::POST;
299 :
300 : // A VisEquation
301 0 : ve_p = new VisEquation();
302 :
303 : // Initialize the meta-info server that will be shared with VisCals
304 0 : if (msmc_p) delete msmc_p;
305 0 : msmc_p = new MSMetaInfoForCal(ssvp_p);
306 :
307 : // Reset the apply/solve VisCals
308 0 : reset();
309 :
310 0 : }
311 :
312 : /*
313 : Calibrater::Calibrater(const Calibrater & other)
314 : {
315 : operator=(other);
316 : }
317 :
318 : Calibrater& Calibrater::operator=(const Calibrater & other)
319 : {
320 : ms_p=other.ms_p;
321 : mssel_p=other.mssel_p;
322 : ve_p=other.ve_p;
323 : histLockCounter_p=other.histLockCounter_p;
324 : hist_p=other.hist_p;
325 : historytab_p=other.historytab_p;
326 :
327 : return *this;
328 : }
329 : */
330 :
331 726 : Calibrater::~Calibrater()
332 : {
333 242 : cleanup();
334 242 : if (msmc_p) delete msmc_p; msmc_p=0;
335 242 : if (ms_p) delete ms_p; ms_p=0;
336 242 : if (hist_p) delete hist_p; hist_p=0;
337 :
338 484 : }
339 :
340 242 : Calibrater* Calibrater::factory(Bool old) {
341 :
342 242 : Calibrater* cal(NULL);
343 :
344 242 : if (old)
345 0 : cal = new OldCalibrater();
346 : else
347 : //throw(AipsError("VI2-aware Calibrater not yet available..."));
348 242 : cal = new Calibrater();
349 :
350 242 : return cal;
351 :
352 : }
353 :
354 0 : Calibrater* Calibrater::factory(String msname, Bool old) {
355 :
356 0 : Calibrater* cal(NULL);
357 :
358 0 : if (old)
359 0 : cal = new OldCalibrater(msname);
360 : else
361 0 : throw(AipsError("VI2-aware Calibrater not yet available..."));
362 : //cal = new Calibrater(msname);
363 :
364 0 : return cal;
365 :
366 : }
367 :
368 :
369 3919 : LogIO& Calibrater::logSink() {return sink_p;};
370 :
371 0 : String Calibrater::timerString() {
372 0 : ostringstream o;
373 0 : o <<" [user: " << timer_p.user () <<
374 0 : " system: " << timer_p.system () <<
375 0 : " real: " << timer_p.real () << "]";
376 0 : timer_p.mark();
377 0 : return o;
378 0 : };
379 :
380 64 : Bool Calibrater::initialize(MeasurementSet& inputMS,
381 : Bool compress,
382 : Bool addScratch, Bool addModel) {
383 :
384 64 : logSink() << LogOrigin("Calibrater","") << LogIO::NORMAL3;
385 :
386 : try {
387 64 : timer_p.mark();
388 :
389 : // Set pointer ms_p from input MeasurementSet
390 64 : if (ms_p) {
391 0 : *ms_p=inputMS;
392 : } else {
393 64 : ms_p = new MeasurementSet(inputMS);
394 64 : AlwaysAssert(ms_p,AipsError);
395 : };
396 :
397 : // Remember the ms's name
398 64 : msname_p=ms_p->tableName();
399 :
400 : // Initialize the meta-info server that will be shared with VisCals
401 64 : if (msmc_p) delete msmc_p;
402 64 : msmc_p = new MSMetaInfoForCal(*ms_p);
403 :
404 : // Add/init scr cols, if requested (init is hard-wired)
405 64 : if (addScratch || addModel) {
406 0 : Bool alsoinit=true;
407 0 : VisSetUtil::addScrCols(*ms_p,addModel,addScratch,alsoinit,compress);
408 : }
409 :
410 : // Set the selected MeasurementSet to be the same initially
411 : // as the input MeasurementSet
412 : logSink() << LogIO::NORMAL
413 : << "Initializing nominal selection to the whole MS."
414 64 : << LogIO::POST;
415 :
416 64 : if (mssel_p) delete mssel_p;
417 64 : mssel_p=new MeasurementSet(*ms_p);
418 :
419 : // Create the associated VisEquation
420 : // TBD: move to ctor and make it non-pointer
421 64 : if (ve_p) {
422 0 : delete ve_p;
423 0 : ve_p=0;
424 : };
425 64 : ve_p=new VisEquation();
426 :
427 : // Reset the apply/solve VisCals
428 64 : reset(true,true);
429 :
430 64 : return true;
431 :
432 0 : } catch (AipsError x) {
433 0 : logSink() << LogOrigin("Calibrater","initialize",WHERE)
434 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
435 0 : << LogIO::POST;
436 0 : cleanup();
437 0 : if (msmc_p) delete msmc_p; msmc_p=NULL;
438 0 : if (ms_p) delete ms_p; ms_p=NULL;
439 0 : if (hist_p) delete hist_p; hist_p=NULL;
440 :
441 0 : throw(AipsError("Error in Calibrater::initialize()"));
442 : return false;
443 0 : }
444 : return false;
445 : }
446 :
447 :
448 : // Select data (using MSSelection syntax)
449 32 : void Calibrater::selectvis(const String& time,
450 : const String& spw,
451 : const String& scan,
452 : const String& field,
453 : const String& intent,
454 : const String& obsIDs,
455 : const String& baseline,
456 : const String& uvrange,
457 : const String& chanmode,
458 : const Int&,
459 : const Int&,
460 : const Int&,
461 : const MRadialVelocity&,
462 : const MRadialVelocity&,
463 : const String& msSelect)
464 : {
465 : // Define primary measurement set selection criteria
466 : // Inputs:
467 : // time
468 : // spw
469 : // scan
470 : // field
471 : // intent
472 : // obsIDs
473 : // baseline
474 : // uvrange
475 : // chanmode const String& Frequency/velocity selection mode
476 : // ("channel", "velocity" or
477 : // "opticalvelocity")
478 : // nchan const Int& No of channels to select
479 : // start const Int& Start channel to select
480 : // step const Int& Channel increment
481 : // mStart const MRadialVelocity& Start radial vel. to select
482 : // mStep const MRadialVelocity& Radial velocity increment
483 : // msSelect const String& MS selection string (TAQL)
484 : // Output to private data:
485 : //
486 32 : logSink() << LogOrigin("Calibrater","selectvis") << LogIO::NORMAL3;
487 :
488 : try {
489 :
490 : /*
491 : cout << "time = " << time << " " << time.length() <<endl;
492 : cout << "spw = " << spw << " " << spw.length() <<endl;
493 : cout << "scan = " << scan << " " << scan.length() <<endl;
494 : cout << "field = " << field << " " << field.length() <<endl;
495 : cout << "baseline = " << baseline << " " << baseline.length() << endl;
496 : cout << "uvrange = " << uvrange << " " << uvrange.length() << endl;
497 : */
498 :
499 32 : logSink() << "Selecting data" << LogIO::POST;
500 :
501 : // Apply selection to the original MeasurementSet
502 32 : logSink() << "Performing selection on MeasurementSet" << endl;
503 :
504 32 : if (mssel_p) {
505 32 : delete mssel_p;
506 32 : mssel_p=0;
507 : };
508 :
509 : // Report non-trivial user selections
510 32 : if (time!="")
511 2 : logSink() << " Selecting on time: '" << time << "'" << endl;
512 32 : if (spw!="")
513 13 : logSink() << " Selecting on spw: '" << spw << "'" << endl;
514 32 : if (scan!="")
515 2 : logSink() << " Selecting on scan: '" << scan << "'" << endl;
516 32 : if (field!="")
517 5 : logSink() << " Selecting on field: '" << field << "'" << endl;
518 32 : if (intent!="")
519 0 : logSink() << " Selecting on intent: '" << intent << "'" << endl;
520 32 : if(obsIDs != "")
521 0 : logSink() << " Selecting by observation IDs: '" << obsIDs << "'" << endl;
522 32 : if (baseline!="")
523 0 : logSink() << " Selecting on antenna/baseline: '" << baseline << "'" << endl;
524 32 : if (uvrange!="")
525 1 : logSink() << " Selecting on uvrange: '" << uvrange << "'" << endl;
526 32 : if (msSelect!="")
527 0 : logSink() << " Selecting with TaQL: '" << msSelect << "'" << endl;
528 32 : logSink() << LogIO::POST;
529 :
530 :
531 : // Assume no selection, for starters
532 32 : mssel_p = new MeasurementSet(*ms_p);
533 :
534 : // Apply user-supplied selection
535 32 : Bool nontrivsel=false;
536 :
537 : // Ensure use of a fresh MSSelection object
538 32 : if (mss_p) { delete mss_p; mss_p=NULL; }
539 32 : mss_p=new MSSelection();
540 64 : nontrivsel= mssSetData(*ms_p,
541 32 : *mssel_p,"",
542 : time,baseline,
543 : field,spw,
544 : uvrange,msSelect,
545 : "",scan,"",intent, obsIDs,mss_p);
546 :
547 : // Keep any MR status for the MS
548 32 : mssel_p->setMemoryResidentSubtables(ms_p->getMrsEligibility());
549 :
550 : // If non-trivial MSSelection invoked and nrow reduced:
551 32 : if(nontrivsel && mssel_p->nrow()<ms_p->nrow()) {
552 :
553 : // Escape if no rows selected
554 14 : if (mssel_p->nrow()==0)
555 0 : throw(AipsError("Specified selection selects zero rows!"));
556 :
557 : // ...otherwise report how many rows are selected
558 14 : logSink() << "By selection " << ms_p->nrow()
559 14 : << " rows are reduced to " << mssel_p->nrow()
560 14 : << LogIO::POST;
561 :
562 : // Revise the msmc_p so it contains only the _selected_ MS
563 14 : if (msmc_p) delete msmc_p;
564 14 : msmc_p = new MSMetaInfoForCal(*mssel_p, ms_p->tableName());
565 :
566 : }
567 : else {
568 : // Selection did nothing:
569 18 : logSink() << "Selection did not drop any rows" << LogIO::POST;
570 : }
571 :
572 : // Attempt to use MSSelection for channel selection
573 : // if user not using the old way
574 32 : if (chanmode=="none") {
575 32 : selectChannel(spw);
576 : }
577 : else {
578 : // Old-fashioned way now deprecated
579 : logSink() << LogIO::WARN
580 : << "You have used the old-fashioned mode parameter" << endl
581 : << "for channel selection." << endl
582 : << "Please begin using the new channel selection" << endl
583 0 : << "syntax in the spw parameter." << LogIO::POST;
584 0 : throw(AipsError("Old-fashioned chanmode selection is no longer supported!"));
585 : }
586 :
587 : }
588 0 : catch (MSSelectionError& x) {
589 : // Re-initialize with the existing MS
590 0 : logSink() << LogOrigin("Calibrater","selectvis",WHERE)
591 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
592 0 : << LogIO::POST;
593 : // jagonzal (CAS-4110): I guess it is not necessary to create these columns when the selection is empty
594 0 : initialize(*ms_p,false,false,false);
595 0 : throw(AipsError("Error in data selection specification: " + x.getMesg()));
596 0 : }
597 0 : catch (AipsError x) {
598 : // Re-initialize with the existing MS
599 0 : logSink() << LogOrigin("Calibrater","selectvis",WHERE)
600 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
601 0 : << LogIO::POST;
602 : // jagonzal (CAS-4110): I guess it is not necessary to create these columns when the selection is empty.
603 0 : initialize(*ms_p,false,false,false);
604 0 : throw(AipsError("Error in Calibrater::selectvis(): " + x.getMesg()));
605 0 : }
606 32 : };
607 :
608 :
609 7 : Bool Calibrater::setapply(const String& type,
610 : const Double& t,
611 : const String& table,
612 : const String& spw,
613 : const String& field,
614 : const String& interp,
615 : const Bool& calwt,
616 : const Vector<Int>& spwmap,
617 : const Vector<Double>& opacity)
618 : {
619 :
620 14 : logSink() << LogOrigin("Calibrater",
621 : "setapply(type, t, table, spw, field, interp, calwt, spwmap, opacity)")
622 14 : << LogIO::NORMAL;
623 :
624 :
625 : // cout << "Calibrater::setapply: field="<< field << endl;
626 :
627 : // Set record format for calibration table application information
628 7 : RecordDesc applyparDesc;
629 7 : applyparDesc.addField ("t", TpDouble);
630 7 : applyparDesc.addField ("table", TpString);
631 7 : applyparDesc.addField ("interp", TpString);
632 7 : applyparDesc.addField ("spw", TpArrayInt);
633 : // applyparDesc.addField ("field", TpArrayInt);
634 7 : applyparDesc.addField ("fieldstr", TpString);
635 7 : applyparDesc.addField ("calwt",TpBool);
636 7 : applyparDesc.addField ("spwmap",TpArrayInt);
637 7 : applyparDesc.addField ("opacity",TpArrayDouble);
638 :
639 : // Create record with the requisite field values
640 7 : Record applypar(applyparDesc);
641 7 : applypar.define ("t", t);
642 7 : applypar.define ("table", table);
643 :
644 : /*
645 : String cinterp=interp;
646 : // cinterp.erase(remove_if(cinterp.begin(), cinterp.end(), isspace), cinterp.end());
647 : cinterp.erase( remove( cinterp.begin(), cinterp.end(), ' ' ), cinterp.end() );
648 : */
649 :
650 7 : applypar.define ("interp", interp);
651 7 : applypar.define ("spw",getSpwIdx(spw));
652 : // applypar.define ("field",getFieldIdx(field));
653 7 : applypar.define ("fieldstr",field);
654 7 : applypar.define ("calwt",calwt);
655 7 : applypar.define ("spwmap",spwmap);
656 7 : applypar.define ("opacity", opacity);
657 :
658 7 : String upType=type;
659 7 : upType.upcase();
660 7 : if (upType=="")
661 : // Get type from table
662 5 : upType = calTableType(table);
663 :
664 14 : return setapply(upType,applypar);
665 :
666 7 : }
667 :
668 7 : Bool Calibrater::setapply (const String& type,
669 : const Record& applypar)
670 : {
671 7 : logSink() << LogOrigin("Calibrater", "setapply(type, applypar)");
672 :
673 : // First try to create the requested VisCal object
674 7 : VisCal *vc(NULL);
675 :
676 : try {
677 :
678 7 : if(!ok())
679 0 : throw(AipsError("Calibrater not prepared for setapply."));
680 :
681 7 : String upType=type;
682 7 : upType.upcase();
683 :
684 : logSink() << LogIO::NORMAL
685 : << "Arranging to APPLY:"
686 7 : << LogIO::POST;
687 :
688 : // Add a new VisCal to the apply list
689 7 : vc = createVisCal(upType,*msmc_p);
690 :
691 7 : vc->setApply(applypar);
692 :
693 : logSink() << LogIO::NORMAL << ". "
694 7 : << vc->applyinfo()
695 7 : << LogIO::POST;
696 :
697 7 : } catch (AipsError x) {
698 : logSink() << LogIO::SEVERE << x.getMesg()
699 : << " Check inputs and try again."
700 0 : << LogIO::POST;
701 0 : if (vc) delete vc;
702 0 : throw(AipsError("Error in Calibrater::setapply."));
703 : return false;
704 0 : }
705 :
706 : // Creation apparently successful, so add to the apply list
707 : // TBD: consolidate with above?
708 : try {
709 :
710 7 : uInt napp=vc_p.nelements();
711 7 : vc_p.resize(napp+1,false,true);
712 7 : vc_p[napp] = vc;
713 7 : vc=NULL;
714 :
715 : // Maintain sort of apply list
716 7 : ve_p->setapply(vc_p);
717 :
718 7 : return true;
719 :
720 0 : } catch (AipsError x) {
721 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
722 0 : << LogIO::POST;
723 0 : if (vc) delete vc;
724 0 : throw(AipsError("Error in Calibrater::setapply."));
725 : return false;
726 0 : }
727 : return false;
728 : }
729 :
730 :
731 : // Validate a Cal Library record
732 0 : Bool Calibrater::validatecallib(Record callib) {
733 :
734 0 : uInt ntab=callib.nfields();
735 0 : for (uInt itab=0;itab<ntab;++itab) {
736 :
737 0 : String tabname=callib.name(itab);
738 0 : Record thistabrec=callib.asRecord(itab);
739 0 : uInt ncl=thistabrec.nfields();
740 :
741 0 : for (uInt icl=0;icl<ncl;++icl) {
742 :
743 0 : if (thistabrec.dataType(icl)!=TpRecord)
744 0 : continue;
745 :
746 0 : Record thisicl=thistabrec.asRecord(icl);
747 : try {
748 0 : CalLibSlice::validateCLS(thisicl);
749 : }
750 0 : catch ( AipsError x) {
751 : logSink() << LogIO::SEVERE
752 : << "Caltable " << tabname
753 : << " is missing the following fields: "
754 : << x.getMesg()
755 0 : << LogIO::POST;
756 0 : }
757 :
758 0 : }
759 0 : }
760 0 : return true;
761 : }
762 :
763 :
764 : // Set up apply-able calibration via a Cal Library
765 2 : Bool Calibrater::setcallib2(Record callib, const MeasurementSet* ms) {
766 :
767 2 : logSink() << LogOrigin("Calibrater", "setcallib2(callib)");
768 :
769 : // cout << "Calibrater::setcallib2(callib) : " << boolalpha << callib << endl;
770 :
771 2 : uInt ntab=callib.nfields();
772 :
773 : // cout << "callib.nfields() = " << ntab << endl;
774 :
775 : // Do some preliminary per-table verification
776 4 : for (uInt itab=0;itab<ntab;++itab) {
777 :
778 2 : String tabname=callib.name(itab);
779 :
780 : // Trap parang
781 : // TBD...
782 : // if (tabname=="<parang>")
783 : // continue;
784 :
785 : // Insist that the table exists on disk
786 2 : if (!Table::isReadable(tabname))
787 0 : throw(AipsError("Caltable "+tabname+" does not exist."));
788 :
789 2 : }
790 :
791 : // Tables exist, so deploy them...
792 :
793 : // Local MS object for callib parsing (only)
794 : // MeasurementSet lms(msname_p,Table::Update);
795 : // TBD: Use selected MS instead (not yet available in OTF plotms context!)
796 : //const MeasurementSet lms(*mssel_p);
797 : //MeasurementSet lms(msname_p);
798 :
799 : // Local const MS object for callib parsing (only)
800 2 : const MeasurementSet *lmsp(0);
801 2 : if (ms) {
802 : // Use supplied MS (from outside), if specified...
803 : // TBD: should we verify same base MS as ms_p/mssel_p?
804 : //cout << "Using externally-specified MS!!" << endl;
805 0 : lmsp=ms;
806 : }
807 : else {
808 : // ...use internal one instead
809 : //cout << "Using internal MS (mssel_p)!!" << endl;
810 2 : lmsp=mssel_p;
811 : }
812 : // Reference for use below
813 2 : const MeasurementSet &lms(*lmsp);
814 :
815 :
816 4 : for (uInt itab=0;itab<ntab;++itab) {
817 :
818 2 : String tabname=callib.name(itab);
819 :
820 : // Get the type from the table
821 2 : String upType=calTableType(tabname);
822 2 : upType.upcase();
823 :
824 : // Add table name to the record
825 2 : Record thistabrec=callib.asrwRecord(itab);
826 2 : thistabrec.define("tablename",tabname);
827 :
828 : // First try to create the requested VisCal object
829 2 : VisCal *vc(NULL);
830 :
831 : try {
832 :
833 : // if(!ok())
834 : // throw(AipsError("Calibrater not prepared for setapply."));
835 :
836 : logSink() << LogIO::NORMAL
837 : << "Arranging to APPLY:"
838 2 : << LogIO::POST;
839 :
840 : // Add a new VisCal to the apply list
841 2 : vc = createVisCal(upType,*msmc_p);
842 :
843 : // ingest this table according to its callib
844 2 : vc->setCallib(thistabrec,lms);
845 :
846 0 : } catch (AipsError x) {
847 : logSink() << LogIO::SEVERE << x.getMesg()
848 : << " Check inputs and try again."
849 0 : << LogIO::POST;
850 0 : if (vc) delete vc;
851 0 : throw(AipsError("Error in Calibrater::callib2."));
852 : return false;
853 0 : }
854 :
855 : // Creation apparently successful, so add to the apply list
856 : // TBD: consolidate with above?
857 : try {
858 :
859 2 : uInt napp=vc_p.nelements();
860 2 : vc_p.resize(napp+1,false,true);
861 2 : vc_p[napp] = vc;
862 2 : vc=NULL;
863 :
864 : // Maintain sort of apply list
865 2 : ve_p->setapply(vc_p);
866 :
867 0 : } catch (AipsError x) {
868 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
869 0 : << LogIO::POST;
870 0 : if (vc) delete vc;
871 0 : throw(AipsError("Error in Calibrater::setapply."));
872 : return false;
873 0 : }
874 2 : }
875 :
876 : // Signal use of CalLibrary
877 2 : usingCalLibrary_=true;
878 :
879 : // All ok, if we get this far!
880 2 : return true;
881 :
882 : }
883 :
884 0 : Bool Calibrater::setmodel(const String& modelImage)
885 : {
886 0 : if (!svc_p)
887 0 : throw(AipsError("Calibrater::setmodel() called before Calibrater::setsolve()"));
888 0 : svc_p->setModel(modelImage);
889 0 : return true;
890 : }
891 :
892 0 : Bool Calibrater::setModel(const Vector<Double>& stokes) {
893 :
894 0 : if (ve_p) {
895 0 : Vector<Float> fstokes(stokes.shape());
896 0 : convertArray(fstokes,stokes);
897 0 : ve_p->setModel(fstokes);
898 0 : }
899 : else
900 0 : throw(AipsError("Error in Calibrater::setModel: no VisEquation."));
901 :
902 0 : return true;
903 :
904 : }
905 :
906 32 : Bool Calibrater::setsolve (const String& type,
907 : const String& solint,
908 : const String& table,
909 : const Bool append,
910 : const Double preavg,
911 : const String& apmode,
912 : const Int minblperant,
913 : const String& refant,
914 : const String& refantmode,
915 : const Bool solnorm,
916 : const String& normtype,
917 : const Float minsnr,
918 : const String& combine,
919 : const Int fillgaps,
920 : const String& cfcache,
921 : const Double painc,
922 : const Int fitorder,
923 : const Float fraction,
924 : const Int numedge,
925 : const String& radius,
926 : const Bool smooth,
927 : const Bool zerorates,
928 : const Bool globalsolve,
929 : const Int niter,
930 : const String& corrcomb,
931 : const Vector<Double>& delaywindow,
932 : const Vector<Double>& ratewindow,
933 : const Vector<Bool>& paramactive,
934 : const Bool concatspws,
935 : const String& solmode,
936 : const Vector<Double>& rmsthresh
937 : )
938 : {
939 :
940 32 : logSink() << LogOrigin("Calibrater","setsolve") << LogIO::NORMAL3;
941 :
942 : // Create a record description containing the solver parameters
943 32 : RecordDesc solveparDesc;
944 32 : solveparDesc.addField ("solint", TpString);
945 32 : solveparDesc.addField ("preavg", TpDouble);
946 32 : solveparDesc.addField ("apmode", TpString);
947 32 : solveparDesc.addField ("refant", TpArrayInt);
948 32 : solveparDesc.addField ("refantmode", TpString);
949 32 : solveparDesc.addField ("minblperant", TpInt);
950 32 : solveparDesc.addField ("table", TpString);
951 32 : solveparDesc.addField ("append", TpBool);
952 32 : solveparDesc.addField ("solnorm", TpBool);
953 32 : solveparDesc.addField ("normtype", TpString);
954 32 : solveparDesc.addField ("type", TpString);
955 32 : solveparDesc.addField ("combine", TpString);
956 32 : solveparDesc.addField ("maxgap", TpInt);
957 32 : solveparDesc.addField ("cfcache", TpString);
958 32 : solveparDesc.addField ("painc", TpDouble);
959 32 : solveparDesc.addField ("fitorder", TpInt);
960 32 : solveparDesc.addField ("solmode", TpString);
961 32 : solveparDesc.addField ("rmsthresh", TpArrayDouble);
962 :
963 : // fringe-fit specific fields
964 32 : solveparDesc.addField ("zerorates", TpBool);
965 32 : solveparDesc.addField ("minsnr", TpFloat);
966 32 : solveparDesc.addField ("globalsolve", TpBool);
967 32 : solveparDesc.addField ("delaywindow", TpArrayDouble);
968 32 : solveparDesc.addField ("ratewindow", TpArrayDouble);
969 32 : solveparDesc.addField ("niter", TpInt);
970 32 : solveparDesc.addField ("corrcomb", TpString);
971 32 : solveparDesc.addField ("paramactive", TpArrayBool);
972 32 : solveparDesc.addField ("concatspws", TpBool);
973 :
974 : // single dish specific fields
975 32 : solveparDesc.addField ("fraction", TpFloat);
976 32 : solveparDesc.addField ("numedge", TpInt);
977 32 : solveparDesc.addField ("radius", TpString);
978 32 : solveparDesc.addField ("smooth", TpBool);
979 :
980 :
981 : // Create a solver record with the requisite field values
982 32 : Record solvepar(solveparDesc);
983 32 : solvepar.define ("solint", solint);
984 32 : solvepar.define ("preavg", preavg);
985 32 : String upmode=apmode;
986 32 : upmode.upcase();
987 32 : solvepar.define ("apmode", upmode);
988 32 : solvepar.define ("refant", getRefantIdxList(refant));
989 32 : solvepar.define ("refantmode", refantmode);
990 32 : solvepar.define ("minblperant", minblperant);
991 32 : solvepar.define ("table", table);
992 32 : solvepar.define ("append", append);
993 32 : solvepar.define ("solnorm", solnorm);
994 32 : solvepar.define ("normtype", normtype);
995 : // Fringe-fit specific
996 32 : solvepar.define ("minsnr", minsnr);
997 32 : solvepar.define ("zerorates", zerorates);
998 32 : solvepar.define ("globalsolve", globalsolve);
999 32 : solvepar.define ("niter", niter);
1000 32 : solvepar.define ("corrcomb", corrcomb);
1001 32 : solvepar.define ("delaywindow", delaywindow);
1002 32 : solvepar.define ("ratewindow", ratewindow);
1003 32 : solvepar.define ("solmode", solmode);
1004 32 : solvepar.define ("rmsthresh", rmsthresh);
1005 32 : solvepar.define ("paramactive", paramactive);
1006 32 : solvepar.define ("concatspws", concatspws);
1007 :
1008 32 : String uptype=type;
1009 32 : uptype.upcase();
1010 32 : solvepar.define ("type", uptype);
1011 :
1012 32 : String upcomb=combine;
1013 32 : upcomb.upcase();
1014 32 : solvepar.define("combine",upcomb);
1015 32 : solvepar.define("maxgap",fillgaps);
1016 32 : solvepar.define ("cfcache", cfcache);
1017 32 : solvepar.define ("painc", painc);
1018 32 : solvepar.define("fitorder", fitorder);
1019 :
1020 : // single dish specific
1021 32 : solvepar.define("fraction", fraction);
1022 32 : solvepar.define("numedge", numedge);
1023 32 : solvepar.define("radius", radius);
1024 32 : solvepar.define("smooth", smooth);
1025 :
1026 64 : return setsolve(type,solvepar);
1027 32 : }
1028 :
1029 0 : Bool Calibrater::setsolvebandpoly(const String& table,
1030 : const Bool& append,
1031 : const String& solint,
1032 : const String& combine,
1033 : const Vector<Int>& degree,
1034 : const Bool& visnorm,
1035 : const Bool& solnorm,
1036 : const Int& maskcenter,
1037 : const Float& maskedge,
1038 : const String& refant) {
1039 :
1040 0 : logSink() << LogOrigin("Calibrater","setsolvebandpoly") << LogIO::NORMAL3;
1041 :
1042 : // TBD: support solution interval!
1043 :
1044 : // Create a record description containing the solver parameters
1045 0 : RecordDesc solveparDesc;
1046 0 : solveparDesc.addField ("table", TpString);
1047 0 : solveparDesc.addField ("append", TpBool);
1048 0 : solveparDesc.addField ("solint", TpString);
1049 0 : solveparDesc.addField ("combine", TpString);
1050 0 : solveparDesc.addField ("degree", TpArrayInt);
1051 0 : solveparDesc.addField ("visnorm", TpBool);
1052 0 : solveparDesc.addField ("solnorm", TpBool);
1053 0 : solveparDesc.addField ("maskcenter", TpInt);
1054 0 : solveparDesc.addField ("maskedge", TpFloat);
1055 0 : solveparDesc.addField ("refant", TpArrayInt);
1056 :
1057 : // solveparDesc.addField ("preavg", TpDouble);
1058 : // solveparDesc.addField ("phaseonly", TpBool);
1059 :
1060 : // Create a solver record with the requisite field values
1061 0 : Record solvepar(solveparDesc);
1062 0 : solvepar.define ("table", table);
1063 0 : solvepar.define ("append", append);
1064 0 : solvepar.define ("solint",solint);
1065 0 : String upcomb=combine;
1066 0 : upcomb.upcase();
1067 0 : solvepar.define ("combine",combine);
1068 0 : solvepar.define ("degree", degree);
1069 0 : solvepar.define ("visnorm", visnorm);
1070 0 : solvepar.define ("solnorm", solnorm);
1071 0 : solvepar.define ("maskcenter", maskcenter);
1072 0 : solvepar.define ("maskedge", maskedge);
1073 0 : solvepar.define ("refant", getRefantIdxList(refant));
1074 :
1075 :
1076 : // solvepar.define ("t", t);
1077 : // solvepar.define ("preavg", preavg);
1078 : // solvepar.define ("phaseonly", phaseonly);
1079 :
1080 :
1081 0 : return setsolve("BPOLY",solvepar);
1082 :
1083 0 : }
1084 :
1085 0 : Bool Calibrater::setsolvegainspline(const String& table,
1086 : const Bool& append,
1087 : const String& apmode,
1088 : const Double& splinetime,
1089 : const Double& preavg,
1090 : const Int& numpoint,
1091 : const Double& phasewrap,
1092 : const String& refant) {
1093 :
1094 0 : logSink() << LogOrigin("Calibrater","setsolvegainspline") << LogIO::NORMAL3;
1095 :
1096 : // Create a record description containing the solver parameters
1097 0 : RecordDesc solveparDesc;
1098 0 : solveparDesc.addField ("table", TpString);
1099 0 : solveparDesc.addField ("append", TpBool);
1100 0 : solveparDesc.addField ("apmode", TpString);
1101 0 : solveparDesc.addField ("splinetime", TpDouble);
1102 0 : solveparDesc.addField ("preavg", TpDouble);
1103 0 : solveparDesc.addField ("refant", TpArrayInt);
1104 0 : solveparDesc.addField ("numpoint", TpInt);
1105 0 : solveparDesc.addField ("phasewrap", TpDouble);
1106 :
1107 : // Create a solver record with the requisite field values
1108 0 : Record solvepar(solveparDesc);
1109 0 : solvepar.define ("table", table);
1110 0 : solvepar.define ("append", append);
1111 0 : String upMode=apmode;
1112 0 : upMode.upcase();
1113 0 : solvepar.define ("apmode", upMode);
1114 0 : solvepar.define ("splinetime",splinetime);
1115 0 : solvepar.define ("preavg", preavg);
1116 0 : solvepar.define ("refant", getRefantIdxList(refant));
1117 0 : solvepar.define ("numpoint",numpoint);
1118 0 : solvepar.define ("phasewrap",phasewrap);
1119 :
1120 0 : return setsolve("GSPLINE",solvepar);
1121 :
1122 0 : }
1123 :
1124 32 : Bool Calibrater::setsolve (const String& type,
1125 : const Record& solvepar) {
1126 :
1127 : // Attempt to create the solvable object
1128 32 : SolvableVisCal *svc(NULL);
1129 : try {
1130 :
1131 32 : if(!ok())
1132 0 : throw(AipsError("Calibrater not prepared for setsolve."));
1133 :
1134 32 : String upType = type;
1135 32 : upType.upcase();
1136 :
1137 : // Clean out any old solve that was lying around
1138 32 : unsetsolve();
1139 :
1140 : logSink() << LogIO::NORMAL
1141 : << "Arranging to SOLVE:"
1142 32 : << LogIO::POST;
1143 :
1144 : // Create the new SolvableVisCal
1145 32 : svc = createSolvableVisCal(upType,*msmc_p);
1146 32 : svc->setSolve(solvepar);
1147 :
1148 : logSink() << LogIO::NORMAL << ". "
1149 32 : << svc->solveinfo()
1150 32 : << LogIO::POST;
1151 :
1152 : // Creation apparently successful, keep it
1153 32 : svc_p=svc;
1154 32 : svc=NULL;
1155 :
1156 : // if calibration specific data filter is necessary
1157 : // keep configuration parameter as a record
1158 32 : setCalFilterConfiguration(upType, solvepar);
1159 :
1160 32 : return true;
1161 :
1162 32 : } catch (AipsError x) {
1163 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1164 0 : << LogIO::POST;
1165 0 : unsetsolve();
1166 0 : if (svc) delete svc;
1167 0 : throw(AipsError("Error in Calibrater::setsolve."));
1168 : return false;
1169 0 : }
1170 : return false;
1171 : }
1172 :
1173 0 : Vector<Int> Calibrater::convertSetToVector(set<Int> selset){
1174 0 : Vector<Int> vtmpint(selset.begin(), selset.size(),0);
1175 0 : return vtmpint;
1176 : }
1177 :
1178 0 : Vector<Int> Calibrater::getSelectedSpws() {
1179 0 : MSMetaData msmdtmp(mssel_p, 50.0);
1180 0 : set<uInt> selspwset = msmdtmp.getSpwIDs();
1181 0 : Vector<Int> res(selspwset.begin(), selspwset.size(),0);
1182 0 : return res;
1183 0 : }
1184 :
1185 0 : Vector<String> Calibrater::getSelectedIntents() {
1186 0 : MSMetaData msmdtmp(mssel_p, 50.0);
1187 0 : set<String> selintentset = msmdtmp.getIntents();
1188 0 : Vector<String> res(selintentset.begin(), selintentset.size(),0);
1189 0 : return res;
1190 0 : }
1191 :
1192 0 : Vector<String> Calibrater::getApplyTables() {
1193 0 : set<String> applytableset;
1194 0 : Int numtables = vc_p.nelements();
1195 0 : string space_delimiter = " ";
1196 0 : vector<string> words{};
1197 0 : size_t pos = 0;
1198 0 : string applyinf;
1199 :
1200 : // parse through the apply info
1201 : // Grab all the tables and put them in a vector
1202 0 : if (numtables > 0){
1203 0 : for (Int i=0;i<numtables;i++){
1204 0 : applyinf = vc_p[i]->applyinfo();
1205 0 : while ((pos=applyinf.find(space_delimiter))!=string::npos) {
1206 0 : words.push_back(applyinf.substr(0, pos));
1207 0 : applyinf.erase(0, pos+space_delimiter.length());
1208 : }
1209 0 : for (const auto &str : words){
1210 0 : if (str.rfind("table=",0) == 0){
1211 0 : applytableset.insert(str.substr(6));
1212 : }
1213 : }
1214 : }
1215 : }
1216 0 : Vector<String> res(applytableset.begin(),applytableset.size(),0);
1217 0 : return res;
1218 0 : }
1219 :
1220 0 : Vector<String> Calibrater::getSolveTable() {
1221 0 : string solveinf;
1222 0 : set<String> solvetableset;
1223 0 : string space_delimiter = " ";
1224 0 : vector<string> words{};
1225 0 : size_t pos = 0;
1226 :
1227 : // Parse through the solve info
1228 : // and extract the applied solve table into a vector
1229 0 : if (svc_p){
1230 0 : solveinf = svc_p->solveinfo();
1231 0 : while ((pos=solveinf.find(space_delimiter)) !=string::npos) {
1232 0 : words.push_back(solveinf.substr(0, pos));
1233 0 : solveinf.erase(0, pos+space_delimiter.length());
1234 : }
1235 0 : for (const auto &str : words) {
1236 0 : if (str.rfind("table=",0) == 0){
1237 0 : solvetableset.insert(str.substr(6));
1238 : }
1239 : }
1240 : }
1241 0 : Vector<String> res(solvetableset.begin(), solvetableset.size(),0);
1242 0 : return res;
1243 0 : }
1244 :
1245 0 : Bool Calibrater::getIteratorSelection(Vector<Int>* observationlist, Vector<Int>* scanlist, Vector<Int>* fieldlist, Vector<Int>* antennalist) {
1246 : // get the observation, scan, field, and antenna from the iterator
1247 0 : set<Int> selobsset;
1248 0 : set<Int> selscanset;
1249 0 : set<Int> selfieldset;
1250 0 : set<Int> selantset;
1251 :
1252 0 : vi::VisibilityIterator2 vi(*mssel_p);
1253 0 : vi.originChunks();
1254 0 : vi.origin();
1255 0 : vi::VisBuffer2 *vb = vi.getVisBuffer();
1256 :
1257 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()){
1258 0 : for (vi.origin (); vi.more(); vi.next()){
1259 : // get obs ids
1260 0 : selobsset.insert(vb->observationId()(0));
1261 : // get scan ids
1262 0 : selscanset.insert(vb->scan()(0));
1263 : // get field ids
1264 0 : selfieldset.insert(vb->fieldId()(0));
1265 : }
1266 : }
1267 0 : for (auto & ant : vb->antenna1()){
1268 0 : selantset.insert(ant);
1269 0 : }
1270 0 : for (auto & ant : vb->antenna2()){
1271 0 : selantset.insert(ant);
1272 0 : }
1273 :
1274 0 : *observationlist = convertSetToVector(selobsset);
1275 0 : *scanlist = convertSetToVector(selscanset);
1276 0 : *fieldlist = convertSetToVector(selfieldset);
1277 0 : *antennalist = convertSetToVector(selantset);
1278 :
1279 0 : return true;
1280 0 : }
1281 :
1282 0 : Record Calibrater::returndict()
1283 : {
1284 : // Create record and subrecrds for selection data
1285 0 : Record rec;
1286 0 : Record selectVisRec;
1287 :
1288 0 : Vector<Int> selscanlist;
1289 0 : Vector<Int> selfieldlist;
1290 0 : Vector<Int> selantlist;
1291 0 : Vector<Int> selobslist;
1292 :
1293 0 : getIteratorSelection(&selobslist, &selscanlist, &selfieldlist, &selantlist);
1294 :
1295 : // Define sub-record for selection parameters
1296 0 : selectVisRec.define("antennas", selantlist);
1297 0 : selectVisRec.define("field", selfieldlist);
1298 0 : selectVisRec.define("spw", getSelectedSpws());
1299 0 : selectVisRec.define("scan", selscanlist);
1300 0 : selectVisRec.define("observation", selobslist);
1301 0 : selectVisRec.define("intents", getSelectedIntents());
1302 :
1303 : // Create a record with current calibrater state information
1304 : //rec.define("antennas", selantlist);
1305 : //rec.define("field", selfieldlist);
1306 : //rec.define("spw", getSelectedSpws());
1307 : //rec.define("scan", selscanlist);
1308 : //rec.define("observation", selobslist);
1309 : //rec.define("intents", getSelectedIntents());
1310 0 : rec.defineRecord("selectvis", selectVisRec);
1311 0 : rec.define("apply_tables", getApplyTables());
1312 0 : rec.define("solve_tables", getSolveTable());
1313 :
1314 0 : rec.merge(resRec_);
1315 :
1316 0 : return rec;
1317 0 : }
1318 :
1319 0 : Bool Calibrater::state() {
1320 :
1321 0 : logSink() << LogOrigin("Calibrater","state") << LogIO::NORMAL3;
1322 :
1323 0 : applystate();
1324 0 : solvestate();
1325 :
1326 0 : return true;
1327 :
1328 : }
1329 :
1330 32 : Bool Calibrater::applystate() {
1331 :
1332 :
1333 : // logSink() << LogOrigin("Calibrater","applystate") << LogIO::NORMAL;
1334 :
1335 : logSink() << LogIO::NORMAL
1336 : << "The following calibration terms are arranged for apply:"
1337 32 : << LogIO::POST;
1338 :
1339 32 : Int napp(vc_p.nelements());
1340 32 : if (napp>0)
1341 17 : for (Int iapp=0;iapp<napp;++iapp)
1342 : logSink() << LogIO::NORMAL << ". "
1343 18 : << vc_p[iapp]->applyinfo()
1344 18 : << LogIO::POST;
1345 : else
1346 : logSink() << LogIO::NORMAL << ". "
1347 : << "(None)"
1348 24 : << LogIO::POST;
1349 :
1350 32 : return true;
1351 :
1352 : }
1353 :
1354 :
1355 32 : Bool Calibrater::solvestate() {
1356 :
1357 : // logSink() << LogOrigin("Calibrater","solvestate") << LogIO::NORMAL;
1358 :
1359 : logSink() << LogIO::NORMAL
1360 : << "The following calibration term is arranged for solve:"
1361 32 : << LogIO::POST;
1362 :
1363 32 : if (svc_p)
1364 : logSink() << LogIO::NORMAL << ". "
1365 32 : << svc_p->solveinfo()
1366 32 : << LogIO::POST;
1367 : else
1368 : logSink() << LogIO::NORMAL << ". "
1369 : << "(None)"
1370 0 : << LogIO::POST;
1371 :
1372 32 : return true;
1373 : }
1374 :
1375 339 : Bool Calibrater::reset(const Bool& apply, const Bool& solve) {
1376 :
1377 : // logSink() << LogOrigin("Calibrater","reset") << LogIO::NORMAL;
1378 :
1379 : // Delete the VisCal apply list
1380 339 : if (apply)
1381 339 : unsetapply();
1382 :
1383 : // Delete the VisCal solve object
1384 339 : if (solve)
1385 339 : unsetsolve();
1386 :
1387 339 : return true;
1388 : }
1389 :
1390 : // Delete all (default) or one VisCal in apply list
1391 339 : Bool Calibrater::unsetapply(const Int& which) {
1392 :
1393 : // logSink() << LogOrigin("Calibrater","unsetapply") << LogIO::NORMAL;
1394 :
1395 : try {
1396 339 : if (which<0) {
1397 348 : for (uInt i=0;i<vc_p.nelements();i++)
1398 9 : if (vc_p[i]) delete vc_p[i];
1399 339 : vc_p.resize(0,true);
1400 : } else {
1401 0 : if (vc_p[which]) delete vc_p[which];
1402 0 : vc_p.remove(which);
1403 : }
1404 :
1405 : // Maintain size/sort of apply list
1406 339 : if(ve_p) ve_p->setapply(vc_p);
1407 :
1408 339 : return true;
1409 0 : } catch (AipsError x) {
1410 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1411 0 : << LogIO::POST;
1412 0 : throw(AipsError("Error in Calibrater::unsetapply."));
1413 :
1414 : return false;
1415 0 : }
1416 : return false;
1417 : }
1418 :
1419 : // Delete solve VisCal
1420 371 : Bool Calibrater::unsetsolve() {
1421 :
1422 : // logSink() << LogOrigin("Calibrater","unsetsolve") << LogIO::NORMAL;
1423 :
1424 : try {
1425 371 : if (svc_p) delete svc_p;
1426 371 : svc_p=NULL;
1427 :
1428 371 : if(ve_p) ve_p->setsolve(*svc_p);
1429 :
1430 371 : return true;
1431 :
1432 0 : } catch (AipsError x) {
1433 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1434 0 : << LogIO::POST;
1435 0 : throw(AipsError("Error in Calibrater::unsetsolve."));
1436 : return false;
1437 0 : }
1438 : return false;
1439 : }
1440 :
1441 : Bool
1442 5 : Calibrater::setCorrDepFlags(const Bool& corrDepFlags)
1443 : {
1444 :
1445 5 : logSink() << LogOrigin("Calibrater","setCorrDepFlags") << LogIO::NORMAL;
1446 :
1447 : // Set it
1448 5 : corrDepFlags_=corrDepFlags;
1449 :
1450 5 : logSink() << "Setting correlation dependent flags = " << (corrDepFlags_ ? "True" : "False") << LogIO::POST;
1451 :
1452 5 : return true;
1453 :
1454 : }
1455 :
1456 : Bool
1457 0 : Calibrater::setCorrcomb(const String& corrcomb)
1458 : {
1459 :
1460 0 : logSink() << LogOrigin("Calibrater", "setCorrcomb") << LogIO::NORMAL;
1461 :
1462 0 : corrcomb_= corrcomb;
1463 :
1464 0 : logSink() << "Setting correlation combination = " << corrcomb << LogIO::POST;
1465 :
1466 0 : return true;
1467 :
1468 : }
1469 :
1470 :
1471 : Bool
1472 0 : Calibrater::correct2(String mode)
1473 : {
1474 0 : logSink() << LogOrigin("Calibrater","correct2 (VI2/VB2)") << LogIO::NORMAL;
1475 :
1476 : //cout << "Artificial STOP!" << endl;
1477 : //return false;
1478 :
1479 :
1480 0 : Bool retval = true;
1481 :
1482 : try {
1483 :
1484 : // Ensure apply list non-zero and properly sorted
1485 0 : ve_p->setapply(vc_p);
1486 :
1487 0 : bool forceOldVIByEnv(false);
1488 0 : forceOldVIByEnv = (getenv("VI1CAL")!=NULL);
1489 0 : if (forceOldVIByEnv && anyEQ(ve_p->listTypes(),VisCal::A)) {
1490 0 : logSink() << LogIO::WARN << "Using VI2 calibration apply. AMueller (uvcontsub) no longer requires VI1 for apply." << LogIO::POST;
1491 : }
1492 :
1493 : /* CAS-12434 (2019Jun07, gmoellen): AMueller works with VB2 in _apply_ (only) context now
1494 : // Trap uvcontsub case, since it does not yet handlg VB2
1495 : if (anyEQ(ve_p->listTypes(),VisCal::A)) {
1496 :
1497 : // Only uvcontsub, nothing else
1498 : if (ve_p->nTerms()==1) {
1499 : // Use old method (which doesn't need WEIGHT_SPECTRUM support in this context)
1500 : return this->correct(mode);
1501 : }
1502 : else
1503 : throw(AipsError("Cannot handle AMueller (uvcontsub) and other types simultaneously."));
1504 : }
1505 : */
1506 :
1507 : // Report the types that will be applied
1508 0 : applystate();
1509 :
1510 : // make mode all-caps
1511 0 : String upmode=mode;
1512 0 : upmode.upcase();
1513 :
1514 : // If trialmode=T, only the flags will be set
1515 : // (and only written if not TRIAL)
1516 0 : Bool trialmode=(upmode.contains("TRIAL") ||
1517 0 : upmode.contains("FLAGONLY"));
1518 :
1519 :
1520 : // Arrange for iteration over data
1521 0 : Block<Int> columns;
1522 : // include scan iteration
1523 0 : columns.resize(5);
1524 0 : columns[0]=MS::ARRAY_ID;
1525 0 : columns[1]=MS::SCAN_NUMBER;
1526 0 : columns[2]=MS::FIELD_ID;
1527 0 : columns[3]=MS::DATA_DESC_ID;
1528 0 : columns[4]=MS::TIME;
1529 :
1530 0 : vi::SortColumns sc(columns);
1531 0 : vi::VisibilityIterator2 vi(*mssel_p,sc,true);
1532 :
1533 : // Apply channel selection (see selectChannel(spw))
1534 0 : if (frequencySelections_p)
1535 0 : vi.setFrequencySelection(*frequencySelections_p);
1536 :
1537 0 : vi::VisBuffer2 *vb = vi.getVisBuffer();
1538 :
1539 : // Detect if we will be setting WEIGHT_SPECTRUM, and arrange for this
1540 0 : vi.originChunks(); // required for wSExists() in next line to work
1541 0 : vi.origin();
1542 0 : Bool doWtSp=vi.weightSpectrumExists(); // Exists non-trivially
1543 :
1544 0 : if (doWtSp && calWt())
1545 0 : logSink() << "Found valid WEIGHT_SPECTRUM, correcting it." << LogIO::POST;
1546 :
1547 : // Pass each timestamp (VisBuffer) to VisEquation for correction
1548 :
1549 0 : Vector<Bool> uncalspw(vi.nSpectralWindows()); // Used to accumulate error messages
1550 0 : uncalspw.set(false); // instead of bombing the user
1551 :
1552 0 : uInt nvb(0);
1553 :
1554 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
1555 :
1556 0 : for (vi.origin(); vi.more(); vi.next()) {
1557 :
1558 0 : uInt spw = vb->spectralWindows()(0);
1559 : //if (ve_p->spwOK(spw)){
1560 : //if ( (usingCalLibrary_ && ve_p->VBOKforCalApply(*vb)) // CalLibrary case
1561 : // || (!usingCalLibrary_ && ve_p->spwOK(spw)) // old-fashioned case
1562 : // ) {
1563 0 : if ( ve_p->VBOKforCalApply(*vb) ) { // Handles old and new (CL) contexts
1564 :
1565 : // Re-initialize weight info from sigma info
1566 : // This is smart wrt spectral weights, etc.
1567 : // (this makes W and WS, if present, "dirty" in the vb)
1568 : // TBD: only do this if !trial (else: avoid the I/O)
1569 0 : vb->resetWeightsUsingSigma();
1570 :
1571 : // Arrange for _in-place_ apply on CORRECTED_DATA (init from DATA)
1572 : // (this makes CD "dirty" in the vb)
1573 : // TBD: only do this if !trial (else: avoid the I/O)
1574 0 : vb->setVisCubeCorrected(vb->visCube());
1575 :
1576 : // Make flagcube dirty in the vb
1577 : // NB: we must _always_ do this I/O (even trial mode)
1578 0 : vb->setFlagCube(vb->flagCube());
1579 :
1580 : // Make all vb "not dirty"; we'll carefully arrange the writeback below
1581 0 : vb->dirtyComponentsClear();
1582 :
1583 : // throws exception if nothing to apply
1584 0 : ve_p->correct2(*vb,trialmode,doWtSp);
1585 :
1586 : // Only if not a trial run, manage writes to disk
1587 0 : if (upmode!="TRIAL") {
1588 :
1589 0 : if (upmode.contains("CAL")) {
1590 0 : vb->setVisCubeCorrected(vb->visCubeCorrected());
1591 :
1592 0 : if (calWt()) {
1593 : // Set weights dirty only if calwt=T
1594 0 : if (doWtSp) {
1595 0 : vb->setWeightSpectrum(vb->weightSpectrum());
1596 : // If WS was calibrated, set W to its channel-axis median
1597 0 : vb->setWeight(partialMedians(vb->weightSpectrum(),IPosition(1,1)));
1598 : }
1599 : else
1600 0 : vb->setWeight(vb->weight());
1601 : }
1602 : }
1603 :
1604 0 : if (upmode.contains("FLAG"))
1605 0 : vb->setFlagCube(vb->flagCube());
1606 :
1607 : // Push the calibrated data etc. back to the MS
1608 0 : vb->writeChangesBack();
1609 0 : nvb+=1; // increment vb counter
1610 : }
1611 :
1612 : }
1613 : else{
1614 0 : uncalspw[spw] = true;
1615 :
1616 : // set the flags, if we are being strict
1617 : // (don't touch the data/weights, which are initialized)
1618 0 : if (upmode.contains("STRICT")) {
1619 :
1620 : // reference (to avoid copy) and set the flags
1621 0 : Cube<Bool> fC(vb->flagCube()); // reference
1622 0 : fC.set(true);
1623 :
1624 : // make dirty for writeChangesBack (does this do an actual copy?)
1625 0 : vb->setFlagCube(vb->flagCube());
1626 :
1627 : // write back to
1628 0 : vb->writeChangesBack();
1629 :
1630 0 : }
1631 : else {
1632 :
1633 : // Break out of inner VI2 loop, if possible
1634 0 : if (! vi.isAsynchronous()){
1635 :
1636 : // Asynchronous I/O doesn't have a way to skip
1637 : // VisBuffers, so only break out when not using
1638 : // async i/o.
1639 0 : break;
1640 :
1641 : }
1642 :
1643 : }
1644 : }
1645 : }
1646 : }
1647 :
1648 : // Now that we're out of the loop, summarize any errors.
1649 :
1650 0 : retval = summarize_uncalspws(uncalspw, "correct",
1651 0 : upmode.contains("STRICT"));
1652 :
1653 0 : actRec_=Record();
1654 0 : actRec_.define("origin","Calibrater::correct");
1655 0 : actRec_.defineRecord("VisEquation",ve_p->actionRec());
1656 :
1657 :
1658 : //cout << "Number of VisBuffers corrected: " << nvb << endl;
1659 :
1660 0 : }
1661 0 : catch (AipsError x) {
1662 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1663 0 : << LogIO::POST;
1664 :
1665 0 : logSink() << "Resetting all calibration application settings." << LogIO::POST;
1666 0 : unsetapply();
1667 :
1668 0 : throw(AipsError("Error in Calibrater::correct."));
1669 : retval = false; // Not that it ever gets here...
1670 0 : }
1671 0 : return retval;
1672 : }
1673 :
1674 0 : Bool Calibrater::corrupt2()
1675 : {
1676 0 : logSink() << LogOrigin("Calibrater","corrupt2 (VI2/VB2)") << LogIO::NORMAL;
1677 0 : Bool retval = true;
1678 :
1679 : try {
1680 :
1681 0 : if (!ok())
1682 0 : throw(AipsError("Calibrater not prepared for corrupt2!"));
1683 :
1684 : // MODEL_DATA column must be physically present!
1685 0 : if (!ms_p->tableDesc().isColumn("MODEL_DATA"))
1686 0 : throw(AipsError("MODEL_DATA column unexpectedly absent. Cannot corrupt."));
1687 :
1688 : // Ensure apply list non-zero and properly sorted
1689 0 : ve_p->setapply(vc_p);
1690 :
1691 : /*
1692 : // Trap uvcontsub case, since it does not yet handlg VB2
1693 : if (anyEQ(ve_p->listTypes(),VisCal::A)) {
1694 :
1695 : // Only uvcontsub, nothing else
1696 : if (ve_p->nTerms()==1) {
1697 : // Use old method (which doesn't need WEIGHT_SPECTRUM support in this context)
1698 : return this->correct(mode);
1699 : }
1700 : else
1701 : throw(AipsError("Cannot handle AMueller (uvcontsub) and other types simultaneously."));
1702 : }
1703 : */
1704 :
1705 : // Report the types that will be applied
1706 0 : applystate();
1707 :
1708 : // Arrange for iteration over data
1709 0 : Block<Int> columns;
1710 : // include scan iteration
1711 0 : columns.resize(5);
1712 0 : columns[0]=MS::ARRAY_ID;
1713 0 : columns[1]=MS::SCAN_NUMBER;
1714 0 : columns[2]=MS::FIELD_ID;
1715 0 : columns[3]=MS::DATA_DESC_ID;
1716 0 : columns[4]=MS::TIME;
1717 :
1718 0 : vi::SortColumns sc(columns);
1719 0 : vi::VisibilityIterator2 vi(*mssel_p,sc,true);
1720 :
1721 : // Apply channel selection (see selectChannel(spw))
1722 0 : if (frequencySelections_p)
1723 0 : vi.setFrequencySelection(*frequencySelections_p);
1724 :
1725 0 : vi::VisBuffer2 *vb = vi.getVisBuffer();
1726 :
1727 : // Detect if we will be setting WEIGHT_SPECTRUM, and arrange for this
1728 0 : vi.originChunks(); // required for wSExists() in next line to work
1729 0 : vi.origin();
1730 :
1731 0 : Vector<Bool> uncalspw(vi.nSpectralWindows()); // Used to accumulate error messages
1732 0 : uncalspw.set(false); // instead of bombing the user
1733 :
1734 : // Pass each timestamp (VisBuffer) to VisEquation for correction
1735 0 : uInt nvb(0);
1736 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
1737 0 : for (vi.origin(); vi.more(); vi.next()) {
1738 0 : uInt spw = vb->spectralWindows()(0);
1739 : //if (ve_p->spwOK(spw)){
1740 0 : if (usingCalLibrary_ || ve_p->spwOK(spw)){
1741 :
1742 : // Make all vb "not dirty", for safety
1743 0 : vb->dirtyComponentsClear();
1744 :
1745 : // throws exception if nothing to apply
1746 0 : ve_p->corrupt2(*vb);
1747 :
1748 : // make visCubeModel _ONLY_ dirty for writeChangesBack
1749 0 : vb->setVisCubeModel(vb->visCubeModel());
1750 :
1751 : // Push the corrupted model back to the MS
1752 0 : vb->writeChangesBack();
1753 0 : nvb+=1; // increment vb counter
1754 : }
1755 : else{
1756 0 : uncalspw[spw] = true;
1757 : }
1758 : }
1759 : }
1760 :
1761 : // Now that we're out of the loop, summarize any errors.
1762 0 : retval = summarize_uncalspws(uncalspw, "corrupt2",false); // (didn't flag them)
1763 :
1764 0 : actRec_=Record();
1765 0 : actRec_.define("origin","Calibrater::corrupt2");
1766 0 : actRec_.defineRecord("VisEquation",ve_p->actionRec());
1767 :
1768 : //cout << "Number of VisBuffers corrupted: " << nvb << endl;
1769 :
1770 0 : }
1771 0 : catch (AipsError x) {
1772 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
1773 0 : << LogIO::POST;
1774 :
1775 0 : logSink() << "Resetting all calibration application settings." << LogIO::POST;
1776 0 : unsetapply();
1777 :
1778 0 : throw(AipsError("Error in Calibrater::corrupt2."));
1779 : retval = false; // Not that it ever gets here...
1780 0 : }
1781 0 : return retval;
1782 : }
1783 :
1784 0 : Bool Calibrater::initWeightsWithTsys(String wtmode, Bool dowtsp,
1785 : String tsystable, String gainfield, String interp, Vector<Int> spwmap) {
1786 :
1787 0 : logSink() << LogOrigin("Calibrater", "initWeightsWithTsys")
1788 0 : << LogIO::NORMAL;
1789 0 : Bool retval = true;
1790 :
1791 : try {
1792 :
1793 0 : if (!ok())
1794 0 : throw(AipsError("Calibrater not prepared for initWeights!"));
1795 :
1796 0 : String uptype = calTableType(tsystable);
1797 0 : if (!uptype.contains("TSYS")) {
1798 0 : throw(AipsError(
1799 0 : "Invalid calibration table type for Tsys weighting."));
1800 : }
1801 : // Set record format for calibration table application information
1802 0 : RecordDesc applyparDesc;
1803 0 : applyparDesc.addField("t", TpDouble);
1804 0 : applyparDesc.addField("table", TpString);
1805 0 : applyparDesc.addField("interp", TpString);
1806 0 : applyparDesc.addField("spw", TpArrayInt);
1807 0 : applyparDesc.addField("fieldstr", TpString);
1808 0 : applyparDesc.addField("calwt", TpBool);
1809 0 : applyparDesc.addField("spwmap", TpArrayInt);
1810 0 : applyparDesc.addField("opacity", TpArrayDouble);
1811 :
1812 : // Create record with the requisite field values
1813 0 : Record applypar(applyparDesc);
1814 0 : applypar.define("t", 0.0);
1815 0 : applypar.define("table", tsystable);
1816 0 : applypar.define("interp", interp);
1817 0 : applypar.define("spw", getSpwIdx(""));
1818 0 : applypar.define("fieldstr", gainfield);
1819 0 : applypar.define("calwt", true);
1820 0 : applypar.define("spwmap", spwmap);
1821 0 : applypar.define("opacity", Vector<Double>(1, 0.0));
1822 :
1823 0 : if (vc_p.nelements() > 0) {
1824 0 : logSink() << LogIO::WARN << "Resetting all calibration application settings." << LogIO::POST;
1825 0 : unsetapply();
1826 : }
1827 0 : logSink() << LogIO::NORMAL << "Weight initialization does not support selection. Resetting MS selection." << LogIO::POST;
1828 0 : selectvis();
1829 0 : StandardTsys vc = StandardTsys(*msmc_p);
1830 0 : vc.setApply(applypar);
1831 :
1832 0 : logSink() << LogIO::NORMAL << ". " << vc.applyinfo() << LogIO::POST;
1833 0 : PtrBlock<VisCal*> vcb(1, &vc);
1834 : // Maintain sort of apply list
1835 0 : ve_p->setapply(vcb);
1836 :
1837 : // Detect WEIGHT_SPECTRUM and SIGMA_SPECTRUM
1838 0 : TableDesc mstd = ms_p->actualTableDesc();
1839 0 : String colWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
1840 0 : Bool wtspexists = mstd.isColumn(colWtSp);
1841 0 : String colSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
1842 0 : Bool sigspexists = mstd.isColumn(colSigSp);
1843 0 : Bool addsigsp = (dowtsp && !sigspexists);
1844 :
1845 : // Some log info
1846 0 : bool use_exposure = false;
1847 0 : if (wtmode == "tsys") {
1848 : logSink()
1849 : << "Initializing SIGMA and WEIGHT according to channel bandwidth and Tsys. NOTE this is an expert mode."
1850 0 : << LogIO::WARN << LogIO::POST;
1851 0 : } else if (wtmode == "tinttsys") {
1852 : logSink()
1853 : << "Initializing SIGMA and WEIGHT according to channel bandwidth, integration time, and Tsys. NOTE this is an expert mode."
1854 0 : << LogIO::WARN << LogIO::POST;
1855 0 : use_exposure = true;
1856 : } else {
1857 0 : throw(AipsError("Unrecognized wtmode specified: " + wtmode));
1858 : }
1859 :
1860 : // Force dowtsp if the column already exists
1861 0 : if (wtspexists && !dowtsp) {
1862 : logSink() << "Found WEIGHT_SPECTRUM; will force its initialization."
1863 0 : << LogIO::POST;
1864 0 : dowtsp = true;
1865 : }
1866 :
1867 : // Report that we are initializing the WEIGHT_SPECTRUM, and prepare to do so.
1868 0 : if (dowtsp) {
1869 :
1870 : // Ensure WEIGHT_SPECTRUM really exists at all
1871 : // (often it exists but is empty)
1872 0 : if (!wtspexists) {
1873 0 : logSink() << "Creating WEIGHT_SPECTRUM." << LogIO::POST;
1874 :
1875 : // Nominal defaulttileshape
1876 0 : IPosition dts(3, 4, 32, 1024);
1877 :
1878 : // Discern DATA's default tile shape and use it
1879 0 : const Record dminfo = ms_p->dataManagerInfo();
1880 0 : for (uInt i = 0; i < dminfo.nfields(); ++i) {
1881 0 : Record col = dminfo.asRecord(i);
1882 : //if (upcase(col.asString("NAME"))=="TILEDDATA") {
1883 0 : if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
1884 0 : dts = IPosition(
1885 0 : col.asRecord("SPEC").asArrayInt(
1886 0 : "DEFAULTTILESHAPE"));
1887 : //cout << "Found DATA's default tile: " << dts << endl;
1888 0 : break;
1889 : }
1890 0 : }
1891 :
1892 : // Add the column
1893 0 : String colWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
1894 0 : TableDesc tdWtSp;
1895 0 : tdWtSp.addColumn(
1896 0 : ArrayColumnDesc<Float>(colWtSp, "weight spectrum", 2));
1897 0 : TiledShapeStMan wtSpStMan("TiledWgtSpectrum", dts);
1898 0 : ms_p->addColumn(tdWtSp, wtSpStMan);
1899 0 : } else
1900 0 : logSink() << "Found WEIGHT_SPECTRUM." << LogIO::POST;
1901 : // Ensure WEIGHT_SPECTRUM really exists at all
1902 : // (often it exists but is empty)
1903 0 : if (!sigspexists) {
1904 0 : logSink() << "Creating SIGMA_SPECTRUM." << LogIO::POST;
1905 :
1906 : // Nominal defaulttileshape
1907 0 : IPosition dts(3, 4, 32, 1024);
1908 :
1909 : // Discern DATA's default tile shape and use it
1910 0 : const Record dminfo = ms_p->dataManagerInfo();
1911 0 : for (uInt i = 0; i < dminfo.nfields(); ++i) {
1912 0 : Record col = dminfo.asRecord(i);
1913 : //if (upcase(col.asString("NAME"))=="TILEDDATA") {
1914 0 : if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
1915 0 : dts = IPosition(
1916 0 : col.asRecord("SPEC").asArrayInt(
1917 0 : "DEFAULTTILESHAPE"));
1918 : //cout << "Found DATA's default tile: " << dts << endl;
1919 0 : break;
1920 : }
1921 0 : }
1922 :
1923 : // Add the column
1924 0 : String colSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
1925 0 : TableDesc tdSigSp;
1926 0 : tdSigSp.addColumn(
1927 0 : ArrayColumnDesc<Float>(colSigSp, "sigma spectrum", 2));
1928 0 : TiledShapeStMan sigSpStMan("TiledSigtSpectrum", dts);
1929 0 : ms_p->addColumn(tdSigSp, sigSpStMan);
1930 : {
1931 0 : TableDesc loctd = ms_p->actualTableDesc();
1932 0 : String loccolSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
1933 0 : AlwaysAssert(loctd.isColumn(loccolSigSp),AipsError);
1934 0 : }
1935 0 : }
1936 : }
1937 : else {
1938 0 : if (sigspexists) {
1939 0 : logSink() << "Removing SIGMA_SPECTRUM for non-channelized weight." << LogIO::POST;
1940 : if (true || ms_p->canRemoveColumn(colSigSp)) {
1941 0 : ms_p->removeColumn(colSigSp);
1942 : }
1943 : else
1944 : logSink() << LogIO::WARN << "Failed to remove SIGMA_SPECTRUM column. Values in SIGMA and SIGMA_SPECTRUM columns may be inconsistent after the operation." << LogIO::POST;
1945 : }
1946 : }
1947 :
1948 : // Arrange for iteration over data
1949 : // TBD: Be sure this sort is optimal for creating WS?
1950 0 : Block<Int> columns;
1951 : // include scan iteration
1952 0 : columns.resize(5);
1953 0 : columns[0] = MS::ARRAY_ID;
1954 0 : columns[1] = MS::SCAN_NUMBER;
1955 0 : columns[2] = MS::FIELD_ID;
1956 0 : columns[3] = MS::DATA_DESC_ID;
1957 0 : columns[4] = MS::TIME;
1958 :
1959 0 : vi::SortColumns sc(columns);
1960 0 : vi::VisibilityIterator2 vi2(*ms_p, sc, true);
1961 0 : vi::VisBuffer2 *vb = vi2.getVisBuffer();
1962 :
1963 0 : MSColumns mscol(*ms_p);
1964 0 : const MSSpWindowColumns& msspw(mscol.spectralWindow());
1965 0 : uInt nSpw = msspw.nrow();
1966 0 : Vector<Double> effChBw(nSpw, 0.0);
1967 0 : for (uInt ispw = 0; ispw < nSpw; ++ispw) {
1968 0 : effChBw[ispw] = mean(msspw.effectiveBW()(ispw));
1969 : }
1970 :
1971 0 : Int ivb(0);
1972 0 : for (vi2.originChunks(); vi2.moreChunks(); vi2.nextChunk()) {
1973 :
1974 0 : for (vi2.origin(); vi2.more(); vi2.next(), ++ivb) {
1975 :
1976 0 : Int spw = vb->spectralWindows()(0);
1977 :
1978 0 : auto nrow = vb->nRows();
1979 0 : Int nchan = vb->nChannels();
1980 0 : Int ncor = vb->nCorrelations();
1981 :
1982 : // Prepare for WEIGHT_SPECTRUM and SIGMA_SPECTRUM, if nec.
1983 0 : Cube<Float> newwtsp(0, 0, 0), newsigsp(0, 0, 0);
1984 0 : if (dowtsp) {
1985 0 : newwtsp.resize(ncor, nchan, nrow);
1986 0 : newwtsp.set(1.0);
1987 0 : newsigsp.resize(ncor, nchan, nrow);
1988 0 : newsigsp.set(1.0);
1989 : }
1990 :
1991 0 : if (ve_p->spwOK(spw)) {
1992 :
1993 : // Re-initialize weight info from sigma info
1994 : // This is smart wrt spectral weights, etc.
1995 : // (this makes W and WS, if present, "dirty" in the vb)
1996 : // TBD: only do this if !trial (else: avoid the I/O)
1997 : // vb->resetWeightsUsingSigma();
1998 : // Handle non-trivial modes
1999 : // Init WEIGHT, SIGMA from bandwidth & time
2000 0 : Matrix<Float> newwt(ncor, nrow), newsig(ncor, nrow);
2001 0 : newwt.set(1.0);
2002 0 : newsig.set(1.0);
2003 :
2004 : // Detect ACs
2005 0 : const Vector<Int> a1(vb->antenna1());
2006 0 : const Vector<Int> a2(vb->antenna2());
2007 0 : Vector<Bool> ac(a1 == a2);
2008 :
2009 : // XCs need an extra factor of 2
2010 0 : Vector<Float> xcfactor(nrow, 2.0);
2011 0 : xcfactor(ac) = 1.0; // (but not ACs)
2012 :
2013 : // The row-wise integration time
2014 0 : Vector<Float> expo(nrow);
2015 0 : convertArray(expo, vb->exposure());
2016 :
2017 : // Set weights to channel bandwidth first.
2018 0 : newwt.set(Float(effChBw(spw)));
2019 :
2020 : // For each correlation, apply exposure and xcfactor
2021 0 : for (Int icor = 0; icor < ncor; ++icor) {
2022 :
2023 0 : Vector<Float> wt(newwt.row(icor));
2024 0 : if (use_exposure) {
2025 0 : wt *= expo;
2026 : }
2027 0 : wt *= xcfactor;
2028 0 : if (dowtsp) {
2029 0 : for (Int ich = 0; ich < nchan; ++ich) {
2030 : Vector<Float> wtspi(
2031 0 : newwtsp(Slice(icor, 1, 1),
2032 0 : Slice(ich, 1, 1), Slice()).nonDegenerate(
2033 0 : IPosition(1, 2)));
2034 0 : wtspi = wt;
2035 0 : }
2036 : }
2037 0 : }
2038 : // Handle SIGMA_SPECTRUM
2039 0 : if (dowtsp) {
2040 0 : newsigsp = 1.0f / sqrt(newwtsp);
2041 : }
2042 : // sig from wt is inverse sqrt
2043 0 : newsig = 1.0f / sqrt(newwt);
2044 :
2045 : // Arrange write-back of both SIGMA and WEIGHT
2046 0 : vb->setSigma(newsig);
2047 0 : vb->setWeight(newwt);
2048 0 : if (dowtsp) {
2049 0 : vb->initWeightSpectrum(newwtsp);
2050 0 : vb->initSigmaSpectrum(newsigsp);
2051 : }
2052 : // Force writeback to disk (need to initialize weight/sigma before applying cal table)
2053 0 : vb->writeChangesBack();
2054 :
2055 : // Arrange for _in-place_ apply on CORRECTED_DATA (init from DATA)
2056 : // (this makes CD "dirty" in the vb)
2057 : // TBD: only do this if !trial (else: avoid the I/O)
2058 0 : vb->setVisCubeCorrected(vb->visCube());
2059 :
2060 : // Make flagcube dirty in the vb
2061 : // NB: we must _always_ do this I/O (even trial mode)
2062 0 : vb->setFlagCube(vb->flagCube());
2063 :
2064 : // Make all vb "not dirty"; we'll carefully arrange the writeback below
2065 0 : vb->dirtyComponentsClear();
2066 :
2067 : // throws exception if nothing to apply
2068 0 : ve_p->correct2(*vb, false, dowtsp);
2069 :
2070 0 : if (dowtsp) {
2071 0 : vb->setWeightSpectrum(vb->weightSpectrum());
2072 : // If WS was calibrated, set W to its channel-axis median
2073 0 : vb->setWeight( partialMedians(vb->weightSpectrum(), IPosition(1, 1)) );
2074 0 : newsigsp = 1.0f / sqrt(vb->weightSpectrum());
2075 0 : vb->initSigmaSpectrum(newsigsp);
2076 0 : vb->setSigma( partialMedians(newsigsp, IPosition(1, 1)) );
2077 : } else {
2078 0 : vb->setWeight(vb->weight());
2079 0 : newsig = 1.0f / sqrt(vb->weight());
2080 0 : vb->setSigma(newsig);
2081 : }
2082 : // Force writeback to disk
2083 0 : vb->writeChangesBack();
2084 :
2085 0 : } else {//Not calibrating the spw
2086 0 : if (dowtsp && !wtspexists) {
2087 : // newly created WS Need to initialize
2088 0 : vb->initWeightSpectrum(newwtsp);
2089 : }
2090 0 : if (addsigsp) {
2091 : // newly created SS Need to initialize
2092 0 : vb->initSigmaSpectrum(newsigsp);
2093 0 : vb->writeChangesBack();
2094 : }
2095 : }
2096 0 : }
2097 : }
2098 : // clear-up Tsys caltable from list of apply
2099 0 : unsetapply();
2100 :
2101 0 : } catch (AipsError x) {
2102 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
2103 0 : << LogIO::POST;
2104 :
2105 : logSink() << "Resetting all calibration application settings."
2106 0 : << LogIO::POST;
2107 0 : unsetapply();
2108 :
2109 0 : throw(AipsError("Error in Calibrater::initWeights."));
2110 : retval = false; // Not that it ever gets here...
2111 0 : }
2112 0 : return retval;
2113 : }
2114 :
2115 :
2116 0 : Bool Calibrater::initWeights(String wtmode, Bool dowtsp) {
2117 :
2118 0 : logSink() << LogOrigin("Calibrater","initWeights") << LogIO::NORMAL;
2119 0 : Bool retval = true;
2120 :
2121 : try {
2122 :
2123 0 : if (!ok())
2124 0 : throw(AipsError("Calibrater not prepared for initWeights!"));
2125 :
2126 : // An enum for the requested wtmode
2127 : enum wtInitModeEnum {
2128 : NONE, // ambiguous, will complain
2129 : ONES, // SIGMA=1.0, propagate to WEIGHT, WEIGHT_SEPCTRUM
2130 : NYQ, // SIGMA=f(df,dt), propagate to WEIGHT, WEIGHT_SEPCTRUM
2131 : SIGMA, // SIGMA as is, propagate to WEIGHT, WEIGHT_SEPCTRUM
2132 : WEIGHT, // SIGMA, WEIGHT as are, propagate to WEIGHT_SEPCTRUM (if requested)
2133 : DELWTSP, // just delete WEIGHT_SPECTRUM if it exists
2134 : DELSIGSP, // just delete SIGMA_SPECTRUM if it exists
2135 : };
2136 :
2137 : // Translate mode String to enum value
2138 0 : wtInitModeEnum initmode(NONE);
2139 0 : if (wtmode=="ones") initmode=ONES;
2140 0 : else if (wtmode=="nyq") initmode=NYQ;
2141 0 : else if (wtmode=="sigma") initmode=SIGMA;
2142 0 : else if (wtmode=="weight") initmode=WEIGHT;
2143 0 : else if (wtmode=="delwtsp") initmode=DELWTSP;
2144 0 : else if (wtmode=="delsigsp") initmode=DELSIGSP;
2145 :
2146 : // Detect WEIGHT_SPECTRUM
2147 0 : TableDesc mstd = ms_p->actualTableDesc();
2148 0 : String colWtSp=MS::columnName(MS::WEIGHT_SPECTRUM);
2149 0 : Bool wtspexists=mstd.isColumn(colWtSp);
2150 0 : String colSigSp=MS::columnName(MS::SIGMA_SPECTRUM);
2151 0 : Bool sigspexists=mstd.isColumn(colSigSp);
2152 :
2153 : // Some log info
2154 0 : switch (initmode) {
2155 0 : case DELWTSP: {
2156 0 : if (wtspexists) {
2157 : if (true || ms_p->canRemoveColumn(colWtSp)) {
2158 0 : logSink() << "Removing WEIGHT_SPECTRUM." << LogIO::POST;
2159 0 : ms_p->removeColumn(colWtSp);
2160 : }
2161 : else
2162 : logSink() << "Sorry, WEIGHT_SPECTRUM is not removable." << LogIO::POST;
2163 : }
2164 : else
2165 0 : logSink() << "WEIGHT_SPECTRUM already absent." << LogIO::POST;
2166 :
2167 : // Nothing more to do here
2168 0 : return true;
2169 : break;
2170 : }
2171 0 : case DELSIGSP: {
2172 0 : if (sigspexists) {
2173 : if (true || ms_p->canRemoveColumn(colSigSp)) {
2174 0 : logSink() << "Removing SIGMA_SPECTRUM." << LogIO::POST;
2175 0 : ms_p->removeColumn(colSigSp);
2176 : }
2177 : else
2178 : logSink() << "Sorry, SIGMA_SPECTRUM is not removable." << LogIO::POST;
2179 : }
2180 : else
2181 0 : logSink() << "SIGMA_SPECTRUM already absent." << LogIO::POST;
2182 :
2183 : // Nothing more to do here
2184 0 : return true;
2185 : break;
2186 : }
2187 0 : case NONE: {
2188 0 : throw(AipsError("Unrecognized wtmode specified: "+wtmode));
2189 : break;
2190 : }
2191 0 : case ONES: {
2192 0 : logSink() << "Initializing SIGMA and WEIGHT with 1.0" << LogIO::POST;
2193 0 : break;
2194 : }
2195 0 : case NYQ: {
2196 : logSink() << "Initializing SIGMA and WEIGHT according to channel bandwidth and integration time."
2197 0 : << LogIO::POST;
2198 0 : break;
2199 : }
2200 0 : case SIGMA: {
2201 0 : logSink() << "Initializing WEIGHT according to existing SIGMA." << LogIO::POST;
2202 0 : break;
2203 : }
2204 0 : case WEIGHT: {
2205 : // Complain if dowtsp==F, because otherwise we have nothing to do
2206 0 : if (!dowtsp)
2207 0 : throw(AipsError("Specified wtmode requires dowtsp=T"));
2208 0 : break;
2209 : }
2210 : }
2211 :
2212 : // Force dowtsp if the column already exists
2213 0 : if (wtspexists && !dowtsp) {
2214 0 : logSink() << "Found WEIGHT_SPECTRUM; will force its initialization." << LogIO::POST;
2215 0 : dowtsp=true;
2216 : }
2217 : // remove SIGMA_SPECTRUM column for non-channelized weight
2218 0 : if (sigspexists) {
2219 0 : logSink() << "Removing SIGMA_SPECTRUM for non-channelized weight." << LogIO::POST;
2220 : if (true || ms_p->canRemoveColumn(colSigSp)) {
2221 0 : ms_p->removeColumn(colSigSp);
2222 : }
2223 : else
2224 : logSink() << LogIO::WARN << "Failed to remove SIGMA_SPECTRUM column. Values in SIGMA and SIGMA_SPECTRUM columns may be inconsistent after the operation." << LogIO::POST;
2225 : }
2226 :
2227 : // Report that we are initializing the WEIGHT_SPECTRUM, and prepare to do so.
2228 0 : if (dowtsp) {
2229 :
2230 : // Ensure WEIGHT_SPECTRUM really exists at all
2231 : // (often it exists but is empty)
2232 0 : if (!wtspexists) {
2233 0 : logSink() << "Creating WEIGHT_SPECTRUM." << LogIO::POST;
2234 :
2235 : // Nominal defaulttileshape
2236 0 : IPosition dts(3,4,32,1024);
2237 :
2238 : // Discern DATA's default tile shape and use it
2239 0 : const Record dminfo=ms_p->dataManagerInfo();
2240 0 : for (uInt i=0;i<dminfo.nfields();++i) {
2241 0 : Record col=dminfo.asRecord(i);
2242 : //if (upcase(col.asString("NAME"))=="TILEDDATA") {
2243 0 : if (anyEQ(col.asArrayString("COLUMNS"),String("DATA"))) {
2244 0 : dts=IPosition(col.asRecord("SPEC").asArrayInt("DEFAULTTILESHAPE"));
2245 : //cout << "Found DATA's default tile: " << dts << endl;
2246 0 : break;
2247 : }
2248 0 : }
2249 :
2250 : // Add the column
2251 0 : String colWtSp=MS::columnName(MS::WEIGHT_SPECTRUM);
2252 0 : TableDesc tdWtSp;
2253 0 : tdWtSp.addColumn(ArrayColumnDesc<Float>(colWtSp,"weight spectrum", 2));
2254 0 : TiledShapeStMan wtSpStMan("TiledWgtSpectrum",dts);
2255 0 : ms_p->addColumn(tdWtSp,wtSpStMan);
2256 0 : }
2257 : else
2258 0 : logSink() << "Found WEIGHT_SPECTRUM." << LogIO::POST;
2259 :
2260 :
2261 0 : logSink() << "Initializing WEIGHT_SPECTRUM uniformly in channel (==WEIGHT)." << LogIO::POST;
2262 :
2263 : }
2264 :
2265 : // Arrange for iteration over data
2266 : // TBD: Be sure this sort is optimal for creating WS?
2267 0 : Block<Int> columns;
2268 : // include scan iteration
2269 0 : columns.resize(5);
2270 0 : columns[0]=MS::ARRAY_ID;
2271 0 : columns[1]=MS::SCAN_NUMBER;
2272 0 : columns[2]=MS::FIELD_ID;
2273 0 : columns[3]=MS::DATA_DESC_ID;
2274 0 : columns[4]=MS::TIME;
2275 :
2276 0 : vi::SortColumns sc(columns);
2277 0 : vi::VisibilityIterator2 vi2(*ms_p,sc,true);
2278 0 : vi::VisBuffer2 *vb = vi2.getVisBuffer();
2279 :
2280 0 : MSColumns mscol(*ms_p);
2281 0 : const MSSpWindowColumns& msspw(mscol.spectralWindow());
2282 0 : uInt nSpw=msspw.nrow();
2283 0 : Vector<Double> effChBw(nSpw,0.0);
2284 0 : for (uInt ispw=0;ispw<nSpw;++ispw) {
2285 0 : effChBw[ispw]=mean(msspw.effectiveBW()(ispw));
2286 : }
2287 :
2288 0 : Int ivb(0);
2289 0 : for (vi2.originChunks(); vi2.moreChunks(); vi2.nextChunk()) {
2290 :
2291 0 : for (vi2.origin(); vi2.more(); vi2.next(),++ivb) {
2292 :
2293 0 : Int spw = vb->spectralWindows()(0);
2294 :
2295 0 : auto nrow=vb->nRows();
2296 0 : Int nchan=vb->nChannels();
2297 0 : Int ncor=vb->nCorrelations();
2298 :
2299 : // Vars for new sigma/weight info
2300 :
2301 : // Prepare for WEIGHT_SPECTRUM, if nec.
2302 0 : Cube<Float> newwtsp(0,0,0);
2303 0 : if (dowtsp) {
2304 0 : newwtsp.resize(ncor,nchan,nrow);
2305 0 : newwtsp.set(1.0);
2306 : }
2307 :
2308 : // Handle non-trivial modes
2309 0 : switch (initmode) {
2310 : // Init WEIGHT, SIGMA with ones
2311 0 : case ONES: {
2312 :
2313 0 : Matrix<Float> newwt(ncor,nrow),newsig(ncor,nrow);
2314 0 : newwt.set(1.0);
2315 0 : newsig.set(1.0);
2316 :
2317 : // Arrange write-back of both SIGMA and WEIGHT
2318 0 : vb->setSigma(newsig);
2319 0 : vb->setWeight(newwt);
2320 :
2321 : // NB: newwtsp already ready
2322 :
2323 0 : break;
2324 0 : }
2325 :
2326 : // Init WEIGHT, SIGMA from bandwidth & time
2327 0 : case NYQ: {
2328 :
2329 0 : Matrix<Float> newwt(ncor,nrow),newsig(ncor,nrow);
2330 0 : newwt.set(1.0);
2331 0 : newsig.set(1.0);
2332 :
2333 : // Detect ACs
2334 0 : const Vector<Int> a1(vb->antenna1());
2335 0 : const Vector<Int> a2(vb->antenna2());
2336 0 : Vector<Bool> ac(a1==a2);
2337 :
2338 : // XCs need an extra factor of 2
2339 0 : Vector<Float> xcfactor(nrow,2.0);
2340 0 : xcfactor(ac)=1.0; // (but not ACs)
2341 :
2342 : // The row-wise integration time
2343 0 : Vector<Float> expo(nrow);
2344 0 : convertArray(expo,vb->exposure());
2345 :
2346 : // Set weights to channel bandwidth first.
2347 0 : newwt.set(Float(effChBw(spw)));
2348 :
2349 : // For each correlation, apply exposure and xcfactor
2350 0 : for (Int icor=0;icor<ncor;++icor) {
2351 :
2352 0 : Vector<Float> wt(newwt.row(icor));
2353 0 : wt*=expo;
2354 0 : wt*=xcfactor;
2355 0 : if (dowtsp) {
2356 0 : for (Int ich=0;ich<nchan;++ich) {
2357 0 : Vector<Float> wtspi(newwtsp(Slice(icor,1,1),Slice(ich,1,1),Slice()).nonDegenerate(IPosition(1,2)));
2358 0 : wtspi=wt;
2359 0 : }
2360 : }
2361 0 : }
2362 :
2363 : // sig from wt is inverse sqrt
2364 0 : newsig=1.0f/sqrt(newwt);
2365 :
2366 : // Arrange write-back of both SIGMA and WEIGHT
2367 0 : vb->setSigma(newsig);
2368 0 : vb->setWeight(newwt);
2369 :
2370 0 : break;
2371 0 : }
2372 : // Init WEIGHT from SIGMA
2373 0 : case SIGMA: {
2374 :
2375 0 : Matrix<Float> newwt(ncor,nrow);
2376 0 : newwt.set(FLT_EPSILON); // effectively zero, but strictly not zero
2377 :
2378 0 : Matrix<Float> sig;
2379 0 : sig.assign(vb->sigma()); // access existing SIGMA
2380 :
2381 0 : LogicalArray mask(sig==0.0f); // mask out unphysical SIGMA
2382 0 : sig(mask)=1.0f;
2383 :
2384 0 : newwt=1.0f/square(sig);
2385 0 : newwt(mask)=FLT_EPSILON; // effectively zero
2386 :
2387 : // Writeback WEIGHT
2388 0 : vb->setWeight(newwt);
2389 :
2390 0 : if (dowtsp) {
2391 0 : for (Int ich=0;ich<nchan;++ich) {
2392 0 : Matrix<Float> wtspi(newwtsp(Slice(),Slice(ich,1,1),Slice()).nonDegenerate(IPosition(2,0,2)));
2393 0 : wtspi=newwt;
2394 0 : }
2395 : }
2396 :
2397 0 : break;
2398 0 : }
2399 : // Init WEIGHT_SPECTRUM from WEIGHT
2400 0 : case WEIGHT: {
2401 0 : const Matrix<Float> wt(vb->weight()); // access existing WEIGHT
2402 0 : for (Int ich=0;ich<nchan;++ich) {
2403 0 : Matrix<Float> wtspi(newwtsp(Slice(),Slice(ich,1,1),Slice()).nonDegenerate(IPosition(2,0,2)));
2404 0 : wtspi=wt;
2405 0 : }
2406 0 : break;
2407 0 : }
2408 0 : default: {
2409 0 : throw(AipsError("Problem in weight initialization loop."));
2410 : }
2411 : }
2412 :
2413 : // Force writeback to disk
2414 0 : vb->writeChangesBack();
2415 :
2416 : // Handle WEIGHT_SPECTRUM
2417 0 : if (dowtsp)
2418 0 : vb->initWeightSpectrum(newwtsp);
2419 :
2420 0 : }
2421 : }
2422 :
2423 0 : }
2424 0 : catch (AipsError x) {
2425 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
2426 0 : << LogIO::POST;
2427 :
2428 0 : logSink() << "Resetting all calibration application settings." << LogIO::POST;
2429 0 : unsetapply();
2430 :
2431 0 : throw(AipsError("Error in Calibrater::initWeights."));
2432 : retval = false; // Not that it ever gets here...
2433 0 : }
2434 0 : return retval;
2435 : }
2436 :
2437 :
2438 0 : Bool Calibrater::initWeights(Bool doBT, Bool dowtsp) {
2439 :
2440 0 : logSink() << LogOrigin("Calibrater","initWeights") << LogIO::NORMAL;
2441 0 : Bool retval = true;
2442 :
2443 : try {
2444 :
2445 0 : if (!ok())
2446 0 : throw(AipsError("Calibrater not prepared for initWeights!"));
2447 :
2448 : // Log that we are beginning...
2449 0 : if (doBT)
2450 0 : logSink() << "Initializing SIGMA and WEIGHT according to channel bandwidth and integration time." << LogIO::POST;
2451 : else
2452 0 : logSink() << "Initializing SIGMA and WEIGHT to unity." << LogIO::POST;
2453 :
2454 0 : if (dowtsp) {
2455 0 : logSink() << "Also initializing WEIGHT_SPECTRUM uniformly in channel (==WEIGHT)." << LogIO::POST;
2456 :
2457 : // Ensure WEIGHT_SPECTRUM really exists at all
2458 : // (often it exists but is empty)
2459 0 : TableDesc mstd = ms_p->actualTableDesc();
2460 0 : if (!mstd.isColumn("WEIGHT_SPECTRUM")) {
2461 0 : cout << "Creating WEIGHT_SPECTRUM!" << endl;
2462 :
2463 : // Nominal defaulttileshape
2464 0 : IPosition dts(3,4,32,1024);
2465 :
2466 : // Discern DATA's default tile shape and use it
2467 0 : const Record dminfo=ms_p->dataManagerInfo();
2468 0 : for (uInt i=0;i<dminfo.nfields();++i) {
2469 0 : Record col=dminfo.asRecord(i);
2470 : //if (upcase(col.asString("NAME"))=="TILEDDATA") {
2471 0 : if (anyEQ(col.asArrayString("COLUMNS"),String("DATA"))) {
2472 0 : dts=IPosition(col.asRecord("SPEC").asArrayInt("DEFAULTTILESHAPE"));
2473 0 : cout << "Found DATA's default tile: " << dts << endl;
2474 0 : break;
2475 : }
2476 0 : }
2477 :
2478 : // Add the column
2479 0 : String colWtSp=MS::columnName(MS::WEIGHT_SPECTRUM);
2480 0 : TableDesc tdWtSp;
2481 0 : tdWtSp.addColumn(ArrayColumnDesc<Float>(colWtSp,"weight spectrum", 2));
2482 0 : TiledShapeStMan wtSpStMan("TiledWgtSpectrum",dts);
2483 0 : ms_p->addColumn(tdWtSp,wtSpStMan);
2484 0 : }
2485 0 : }
2486 :
2487 : // Arrange for iteration over data
2488 0 : Block<Int> columns;
2489 : // include scan iteration
2490 0 : columns.resize(5);
2491 0 : columns[0]=MS::ARRAY_ID;
2492 0 : columns[1]=MS::SCAN_NUMBER;
2493 0 : columns[2]=MS::FIELD_ID;
2494 0 : columns[3]=MS::DATA_DESC_ID;
2495 0 : columns[4]=MS::TIME;
2496 :
2497 0 : vi::SortColumns sc(columns);
2498 0 : vi::VisibilityIterator2 vi2(*ms_p,sc,true);
2499 0 : vi::VisBuffer2 *vb = vi2.getVisBuffer();
2500 :
2501 0 : MSColumns mscol(*ms_p);
2502 0 : const MSSpWindowColumns& msspw(mscol.spectralWindow());
2503 0 : uInt nSpw=msspw.nrow();
2504 0 : Vector<Double> effChBw(nSpw,0.0);
2505 0 : for (uInt ispw=0;ispw<nSpw;++ispw) {
2506 0 : effChBw[ispw]=mean(msspw.effectiveBW()(ispw));
2507 : }
2508 :
2509 0 : Int ivb(0);
2510 0 : for (vi2.originChunks(); vi2.moreChunks(); vi2.nextChunk()) {
2511 :
2512 0 : for (vi2.origin(); vi2.more(); vi2.next()) {
2513 :
2514 0 : Int spw = vb->spectralWindows()(0);
2515 :
2516 0 : auto nrow=vb->nRows();
2517 0 : Int nchan=vb->nChannels();
2518 0 : Int ncor=vb->nCorrelations();
2519 :
2520 0 : Matrix<Float> newwt(ncor,nrow),newsig(ncor,nrow);
2521 0 : newwt.set(1.0);
2522 0 : newsig.set(1.0);
2523 :
2524 0 : Cube<Float> newwtsp(0,0,0);
2525 0 : if (dowtsp) {
2526 0 : newwtsp.resize(ncor,nchan,nrow);
2527 0 : newwtsp.set(1.0);
2528 : }
2529 :
2530 : // If requested, set weights according to bandwidth and integration time
2531 0 : if (doBT) {
2532 :
2533 : // Detect ACs
2534 0 : const Vector<Int> a1(vb->antenna1());
2535 0 : const Vector<Int> a2(vb->antenna2());
2536 0 : Vector<Bool> ac(a1==a2);
2537 :
2538 : // XCs need an extra factor of 2
2539 0 : Vector<Float> xcfactor(nrow,2.0);
2540 0 : xcfactor(ac)=1.0; // (but not ACs)
2541 :
2542 : // The row-wise integration time
2543 0 : Vector<Float> expo(nrow);
2544 0 : convertArray(expo,vb->exposure());
2545 :
2546 : // Set weights to channel bandwidth first.
2547 0 : newwt.set(Float(effChBw(spw)));
2548 :
2549 : // For each correlation, apply exposure and xcfactor
2550 0 : for (Int icor=0;icor<ncor;++icor) {
2551 :
2552 0 : Vector<Float> wt(newwt.row(icor));
2553 0 : wt*=expo;
2554 0 : wt*=xcfactor;
2555 0 : if (dowtsp) {
2556 0 : for (Int ich=0;ich<nchan;++ich) {
2557 0 : Vector<Float> wtsp(newwtsp(Slice(icor,1,1),Slice(ich,1,1),Slice()));
2558 0 : wtsp=wt;
2559 0 : }
2560 : }
2561 :
2562 0 : }
2563 :
2564 : // sig from wt is inverse sqrt
2565 0 : newsig=newsig/sqrt(newwt);
2566 :
2567 0 : }
2568 :
2569 : /*
2570 : cout << ivb << " "
2571 : << ncor << " " << nchan << " " << nrow << " "
2572 : << expo(0) << " "
2573 : << newwt(0,0) << " "
2574 : << newsig(0,0) << " "
2575 : << endl;
2576 : */
2577 0 : ++ivb;
2578 :
2579 : // Set in vb, and writeback
2580 0 : vb->setWeight(newwt);
2581 0 : vb->setSigma(newsig);
2582 0 : vb->writeChangesBack();
2583 :
2584 0 : if (dowtsp)
2585 0 : vb->initWeightSpectrum(newwtsp);
2586 :
2587 0 : }
2588 : }
2589 :
2590 0 : }
2591 0 : catch (AipsError x) {
2592 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
2593 0 : << LogIO::POST;
2594 :
2595 0 : logSink() << "Resetting all calibration application settings." << LogIO::POST;
2596 0 : unsetapply();
2597 :
2598 0 : throw(AipsError("Error in Calibrater::initWeightsOLD."));
2599 : retval = false; // Not that it ever gets here...
2600 0 : }
2601 0 : return retval;
2602 : }
2603 :
2604 :
2605 32 : Bool Calibrater::solve() {
2606 :
2607 32 : logSink() << LogOrigin("Calibrater","solve") << LogIO::NORMAL3;
2608 :
2609 : try {
2610 :
2611 32 : if (!ok())
2612 0 : throw(AipsError("Calibrater not prepared for solve."));
2613 :
2614 : // Handle nothing-to-solve-for case
2615 32 : if (!svc_p)
2616 0 : throw(AipsError("Please run setsolve before attempting to solve."));
2617 :
2618 : // Handle specified caltable
2619 32 : if (svc_p) {
2620 :
2621 : // If table exists (readable) and not deletable
2622 : // we have to abort (append=T requires deletable)
2623 46 : if ( Table::isReadable(svc_p->calTableName()) &&
2624 14 : !TableUtil::canDeleteTable(svc_p->calTableName()) ) {
2625 0 : throw(AipsError("Specified caltable ("+svc_p->calTableName()+") exists and\n cannot be replaced (or appended to) because it appears to be open somewhere (Quit plotcal?)."));
2626 : }
2627 : }
2628 :
2629 : // Arrange VisEquation for solve
2630 32 : ve_p->setsolve(*svc_p);
2631 :
2632 : // Ensure apply list properly sorted w.r.t. solvable term
2633 32 : ve_p->setapply(vc_p);
2634 :
2635 : // Report what is being applied and solved-for
2636 32 : applystate();
2637 32 : solvestate();
2638 :
2639 : // Report correct/corrupt apply order
2640 : // ve_p->state();
2641 :
2642 : // Generally use standard solver
2643 32 : if (svc_p->useGenericGatherForSolve())
2644 32 : genericGatherAndSolve(); // VI2-driven, SDBs
2645 : else {
2646 : //cout << "Fully self-directed data gather and solve" << endl;
2647 : // Fully self-directed data gather and solve
2648 0 : throw(AipsError("Can't do selfGatherAndSolve for "+svc_p->typeName()));
2649 : //svc_p->selfGatherAndSolve(*vs_p,*ve_p);
2650 : }
2651 :
2652 1 : } catch (AipsError x) {
2653 1 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2654 :
2655 1 : logSink() << "Reseting entire solve/apply state." << LogIO::POST;
2656 1 : reset();
2657 :
2658 1 : throw(AipsError("Error in Calibrater::solve."));
2659 : return false;
2660 2 : }
2661 :
2662 31 : return true;
2663 :
2664 : }
2665 :
2666 :
2667 31 : Bool Calibrater::summarize_uncalspws(const Vector<Bool>& uncalspw,
2668 : const String& origin,
2669 : Bool strictflag)
2670 : {
2671 31 : Bool hadprob = false;
2672 31 : uInt totNspw = uncalspw.nelements();
2673 :
2674 155 : for(uInt i = 0; i < totNspw; ++i){
2675 124 : if(uncalspw[i]){
2676 0 : hadprob = true;
2677 0 : break;
2678 : }
2679 : }
2680 31 : if(hadprob){
2681 : logSink() << LogIO::WARN
2682 0 : << "Spectral window(s) ";
2683 0 : for(uInt i = 0; i < totNspw; ++i){
2684 0 : if(uncalspw[i]){
2685 0 : logSink() << i << ", ";
2686 : }
2687 : }
2688 : logSink() << "\n could not be " << origin << "ed due to missing (pre-)calibration\n"
2689 0 : << " in one or more of the specified tables.\n";
2690 0 : if (strictflag)
2691 0 : logSink() << " These spws have been flagged!";
2692 : else
2693 0 : logSink() << " Please check your results carefully!";
2694 :
2695 0 : logSink() << LogIO::POST;
2696 : }
2697 31 : return !hadprob;
2698 : }
2699 :
2700 29 : void Calibrater::fluxscale(const String& infile,
2701 : const String& outfile,
2702 : const String& refFields,
2703 : const Vector<Int>& refSpwMap,
2704 : const String& tranFields,
2705 : const Bool& append,
2706 : const Float& inGainThres,
2707 : const String& antSel,
2708 : const String& timerangeSel,
2709 : const String& scanSel,
2710 : SolvableVisCal::fluxScaleStruct& oFluxScaleFactor,
2711 : Vector<Int>& tranidx,
2712 : const String& oListFile,
2713 : const Bool& incremental,
2714 : const Int& fitorder,
2715 : const Bool& display) {
2716 :
2717 : // TBD: Permit more flexible matching on specified field names
2718 : // (Currently, exact matches are required.)
2719 :
2720 29 : logSink() << LogOrigin("Calibrater","fluxscale") << LogIO::NORMAL3;
2721 :
2722 : //outfile check
2723 29 : if (outfile=="")
2724 0 : throw(AipsError("output fluxscaled caltable name must be specified!"));
2725 : else {
2726 29 : if (File(outfile).exists() && !append)
2727 0 : throw(AipsError("output caltable name, "+outfile+" exists. Please specify a different caltable name"));
2728 : }
2729 :
2730 : // Convert refFields/transFields to index lists
2731 29 : Vector<Int> refidx(0);
2732 :
2733 29 : if (refFields.length()>0)
2734 29 : refidx=getFieldIdx(refFields);
2735 : else
2736 0 : throw(AipsError("A reference field must be specified!"));
2737 :
2738 29 : if (tranFields.length()>0)
2739 12 : tranidx=getFieldIdx(tranFields);
2740 :
2741 : // Call Vector<Int> version:
2742 29 : fluxscale(infile,outfile,refidx,refSpwMap,tranidx,append,inGainThres,antSel,timerangeSel,
2743 : scanSel,oFluxScaleFactor,oListFile,incremental,fitorder,display);
2744 :
2745 29 : }
2746 :
2747 29 : void Calibrater::fluxscale(const String& infile,
2748 : const String& outfile,
2749 : const Vector<Int>& refField,
2750 : const Vector<Int>& refSpwMap,
2751 : const Vector<Int>& tranField,
2752 : const Bool& append,
2753 : const Float& inGainThres,
2754 : const String& antSel,
2755 : const String& timerangeSel,
2756 : const String& scanSel,
2757 : SolvableVisCal::fluxScaleStruct& oFluxScaleFactor,
2758 : const String& oListFile,
2759 : const Bool& incremental,
2760 : const Int& fitorder,
2761 : const Bool& display) {
2762 :
2763 : // throw(AipsError("Method 'fluxscale' is temporarily disabled."));
2764 :
2765 : // TBD: write inputs to MSHistory
2766 29 : logSink() << LogOrigin("Calibrater","fluxscale") << LogIO::NORMAL3;
2767 :
2768 29 : SolvableVisCal *fsvj_(NULL);
2769 : try {
2770 : // If infile is Calibration table
2771 58 : if (Table::isReadable(infile) &&
2772 58 : TableUtil::tableInfo(infile).type()=="Calibration") {
2773 :
2774 : // get calibration type
2775 29 : String caltype;
2776 29 : caltype = TableUtil::tableInfo(infile).subType();
2777 : logSink() << "Table " << infile
2778 : << " is of type: "<< caltype
2779 29 : << LogIO::POST;
2780 58 : String message="Table "+infile+" is of type: "+caltype;
2781 29 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
2782 :
2783 : // form selection
2784 29 : String select="";
2785 : // Selection is empty for case of no tran specification
2786 29 : if (tranField.nelements()>0) {
2787 :
2788 : // All selected fields
2789 12 : Vector<Int> allflds = concatenateArray(refField,tranField);
2790 :
2791 : // Assemble TaQL
2792 12 : ostringstream selectstr;
2793 12 : selectstr << "FIELD_ID IN [";
2794 42 : for (Int iFld=0; iFld<allflds.shape(); iFld++) {
2795 30 : if (iFld>0) selectstr << ", ";
2796 30 : selectstr << allflds(iFld);
2797 : }
2798 12 : selectstr << "]";
2799 12 : select=selectstr.str();
2800 12 : }
2801 :
2802 : // Construct proper SVC object
2803 29 : if (caltype == "G Jones") {
2804 29 : fsvj_ = createSolvableVisCal("G",*msmc_p);
2805 0 : } else if (caltype == "T Jones") {
2806 0 : fsvj_ = createSolvableVisCal("T",*msmc_p);
2807 : } else {
2808 : // Can't process other than G and T (add B?)
2809 0 : ostringstream typeErr;
2810 0 : typeErr << "Type " << caltype
2811 0 : << " not supported in fluxscale.";
2812 :
2813 0 : throw(AipsError(typeErr.str()));
2814 0 : }
2815 :
2816 : // fill table with selection
2817 29 : RecordDesc applyparDesc;
2818 29 : applyparDesc.addField ("table", TpString);
2819 29 : applyparDesc.addField ("select", TpString);
2820 29 : Record applypar(applyparDesc);
2821 29 : applypar.define ("table", infile);
2822 29 : applypar.define ("select", select);
2823 29 : fsvj_->setApply(applypar);
2824 :
2825 : //Bool incremental=false;
2826 : // Make fluxscale calculation
2827 29 : Vector<String> fldnames(MSFieldColumns(ms_p->field()).name().getColumn());
2828 : //fsvj_->fluxscale(refField,tranField,refSpwMap,fldnames,oFluxScaleFactor,
2829 29 : fsvj_->fluxscale(outfile,refField,tranField,refSpwMap,fldnames,inGainThres,antSel,
2830 : timerangeSel,scanSel,oFluxScaleFactor, oListFile,incremental,fitorder,display);
2831 : // oListFile);
2832 :
2833 : // If no outfile specified, use infile (overwrite!)
2834 29 : String out(outfile);
2835 29 : if (out.length()==0)
2836 0 : out = infile;
2837 :
2838 : // Store result
2839 29 : if (append) {
2840 1 : logSink() << "Appending result to " << out << LogIO::POST;
2841 1 : String message="Appending result to "+out;
2842 1 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
2843 1 : } else {
2844 28 : logSink() << "Storing result in " << out << LogIO::POST;
2845 28 : String message="Storing result in "+out;
2846 28 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
2847 28 : }
2848 29 : fsvj_->storeNCT(out,append);
2849 :
2850 : // Clean up
2851 29 : delete fsvj_;
2852 :
2853 29 : } else {
2854 : // Table not found/unreadable, or not Calibration
2855 0 : ostringstream tabErr;
2856 0 : tabErr << "File " << infile
2857 0 : << " does not exist or is not a Calibration Table.";
2858 :
2859 0 : throw(AipsError(tabErr.str()));
2860 :
2861 0 : }
2862 0 : } catch (AipsError x) {
2863 :
2864 : logSink() << LogIO::SEVERE
2865 : << "Caught Exception: "
2866 : << x.getMesg()
2867 0 : << LogIO::POST;
2868 :
2869 : // Clean up
2870 0 : if (fsvj_) delete fsvj_;
2871 :
2872 : // Write to MS History table
2873 : // String message="Caught Exception: "+x.getMesg();
2874 : // MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
2875 :
2876 0 : throw(AipsError("Error in Calibrater::fluxscale."));
2877 :
2878 : return;
2879 :
2880 0 : }
2881 29 : return;
2882 :
2883 :
2884 : }
2885 :
2886 0 : void Calibrater::specifycal(const String& type,
2887 : const String& caltable,
2888 : const String& time,
2889 : const String& spw,
2890 : const String& antenna,
2891 : const String& pol,
2892 : const Vector<Double>& parameter,
2893 : const String& infile,
2894 : const Bool& uniform) {
2895 :
2896 0 : logSink() << LogOrigin("Calibrater","specifycal") << LogIO::NORMAL;
2897 :
2898 : // SVJ objects:
2899 0 : SolvableVisCal *cal_(NULL);
2900 :
2901 : try {
2902 :
2903 : // Set record format for calibration table application information
2904 0 : RecordDesc specifyDesc;
2905 0 : specifyDesc.addField ("caltable", TpString);
2906 0 : specifyDesc.addField ("time", TpString);
2907 0 : specifyDesc.addField ("spw", TpArrayInt);
2908 0 : specifyDesc.addField ("antenna", TpArrayInt);
2909 0 : specifyDesc.addField ("pol", TpString);
2910 0 : specifyDesc.addField ("parameter", TpArrayDouble);
2911 0 : specifyDesc.addField ("caltype",TpString);
2912 0 : specifyDesc.addField ("infile",TpString);
2913 0 : specifyDesc.addField ("uniform",TpBool);
2914 :
2915 : // Create record with the requisite field values
2916 0 : Record specify(specifyDesc);
2917 0 : specify.define ("caltable", caltable);
2918 0 : specify.define ("time", time);
2919 0 : if (spw=="*")
2920 0 : specify.define ("spw",Vector<Int>(1,-1));
2921 : else
2922 0 : specify.define ("spw",getSpwIdx(spw));
2923 0 : if (antenna=="*")
2924 0 : specify.define ("antenna",Vector<Int>(1,-1) );
2925 : else
2926 0 : specify.define ("antenna",getAntIdx(antenna));
2927 0 : specify.define ("pol",pol);
2928 0 : specify.define ("parameter",parameter);
2929 0 : specify.define ("caltype",type);
2930 0 : specify.define ("infile",infile);
2931 0 : specify.define ("uniform",uniform);
2932 :
2933 : // Now do it
2934 0 : String utype=upcase(type);
2935 0 : if (utype=="G" || utype.contains("AMP") || utype.contains("PH"))
2936 0 : cal_ = createSolvableVisCal("G",*msmc_p);
2937 0 : else if (utype=='K' || utype.contains("SBD") || utype.contains("DELAY"))
2938 0 : cal_ = createSolvableVisCal("K",*msmc_p);
2939 0 : else if (utype.contains("MBD"))
2940 0 : cal_ = createSolvableVisCal("K",*msmc_p); // As of 5.3, MBD is ordinary K
2941 0 : else if (utype.contains("ANTPOS"))
2942 0 : cal_ = createSolvableVisCal("KANTPOS",*msmc_p);
2943 0 : else if (utype.contains("TSYS"))
2944 0 : cal_ = createSolvableVisCal("TSYS",*msmc_p);
2945 0 : else if (utype.contains("EVLAGAIN") ||
2946 0 : utype.contains("SWP") ||
2947 0 : utype.contains("RQ") ||
2948 0 : utype.contains("WTS"))
2949 0 : cal_ = createSolvableVisCal("EVLASWP",*msmc_p);
2950 0 : else if (utype.contains("OPAC"))
2951 0 : cal_ = createSolvableVisCal("TOPAC",*msmc_p);
2952 0 : else if (utype.contains("GC") && ms_p && ms_p->keywordSet().isDefined("GAIN_CURVE"))
2953 0 : cal_ = createSolvableVisCal("POWERCURVE",*msmc_p);
2954 0 : else if (utype.contains("GC") || utype.contains("EFF"))
2955 0 : cal_ = createSolvableVisCal("GAINCURVE",*msmc_p);
2956 0 : else if (utype.contains("TEC"))
2957 0 : cal_ = createSolvableVisCal("TEC",*msmc_p);
2958 0 : else if (utype.contains("SDSKY_PS"))
2959 0 : cal_ = createSolvableVisCal("SDSKY_PS",*msmc_p);
2960 0 : else if (utype.contains("SDSKY_RASTER"))
2961 0 : cal_ = createSolvableVisCal("SDSKY_RASTER",*msmc_p);
2962 0 : else if (utype.contains("SDSKY_OTF"))
2963 0 : cal_ = createSolvableVisCal("SDSKY_OTF",*msmc_p);
2964 : else
2965 0 : throw(AipsError("Unrecognized caltype."));
2966 :
2967 : // set up for specification (set up the CalSet)
2968 0 : cal_->setSpecify(specify);
2969 :
2970 : // fill with specified values
2971 0 : cal_->specify(specify);
2972 :
2973 : // Store result
2974 0 : cal_->storeNCT();
2975 :
2976 0 : delete cal_;
2977 :
2978 0 : } catch (AipsError x) {
2979 : logSink() << LogIO::SEVERE
2980 : << "Caught Exception: "
2981 : << x.getMesg()
2982 0 : << LogIO::POST;
2983 :
2984 0 : if (cal_) delete cal_;
2985 :
2986 0 : throw(AipsError("Error in Calibrater::specifycal."));
2987 : return;
2988 0 : }
2989 0 : return;
2990 :
2991 : }
2992 :
2993 :
2994 0 : Bool Calibrater::smooth(const String& infile,
2995 : String& outfile, // const Bool& append,
2996 : const String& smoothtype,
2997 : const Double& smoothtime,
2998 : const String& fields,
2999 : const bool ratesmooth)
3000 : {
3001 :
3002 : // TBD: support append?
3003 : // TBD: spw selection?
3004 :
3005 0 : logSink() << LogOrigin("Calibrater","smooth") << LogIO::NORMAL;
3006 :
3007 0 : logSink() << "Beginning smoothing/interpolating method." << LogIO::POST;
3008 :
3009 :
3010 : // A pointer to an SVC
3011 0 : SolvableVisCal *svc(NULL);
3012 :
3013 : try {
3014 :
3015 : // Handle no in file
3016 0 : if (infile=="")
3017 0 : throw(AipsError("Please specify an input calibration table."));
3018 :
3019 : // Handle bad smoothtype
3020 0 : if (smoothtype!="mean" && smoothtype!="median")
3021 0 : throw(AipsError("Unrecognized smooth type!"));
3022 :
3023 : // Handle bad smoothtime
3024 0 : if (smoothtime<=0)
3025 0 : throw(AipsError("Please specify a strictly positive smoothtime."));
3026 :
3027 : // Handle no outfile
3028 0 : if (outfile=="") {
3029 0 : outfile=infile;
3030 : logSink() << "Will overwrite input file with smoothing result."
3031 0 : << LogIO::POST;
3032 : }
3033 :
3034 :
3035 0 : svc = createSolvableVisCal(calTableType(infile),*msmc_p);
3036 :
3037 0 : if (svc->smoothable()) {
3038 :
3039 : // Fill calibration table using setApply
3040 0 : RecordDesc applyparDesc;
3041 0 : applyparDesc.addField ("table", TpString);
3042 0 : Record applypar(applyparDesc);
3043 0 : applypar.define ("table", infile);
3044 0 : svc->setApply(applypar);
3045 :
3046 : // Convert refFields/transFields to index lists
3047 0 : Vector<Int> fldidx(0);
3048 0 : if (fields.length()>0)
3049 0 : fldidx=getFieldIdx(fields);
3050 :
3051 : // Delegate to SVC
3052 0 : svc->smooth(fldidx,smoothtype,smoothtime,ratesmooth);
3053 :
3054 : // Store the result on disk
3055 : // if (append) logSink() << "Appending result to " << outfile << LogIO::POST;
3056 : //else
3057 0 : logSink() << "Storing result in " << outfile << LogIO::POST;
3058 :
3059 :
3060 0 : if (outfile != "")
3061 0 : svc->calTableName()=outfile;
3062 0 : svc->storeNCT();
3063 :
3064 : // Clean up
3065 0 : if (svc) delete svc; svc=NULL;
3066 :
3067 : // Apparently, it worked
3068 0 : return true;
3069 :
3070 0 : }
3071 : else
3072 0 : throw(AipsError("This type ("+svc->typeName()+") does not support smoothing."));
3073 :
3074 0 : } catch (AipsError x) {
3075 :
3076 : logSink() << LogIO::SEVERE
3077 : << "Caught Exception: "
3078 : << x.getMesg()
3079 0 : << LogIO::POST;
3080 : // Clean up
3081 0 : if (svc) delete svc; svc=NULL;
3082 :
3083 0 : throw(AipsError("Error in Calibrater::smooth."));
3084 :
3085 : return false;
3086 0 : }
3087 : return false;
3088 : }
3089 :
3090 : // Apply new reference antenna to calibration
3091 0 : Bool Calibrater::reRefant(const casacore::String& infile,
3092 : casacore::String& outfile,
3093 : const casacore::String& refantmode,
3094 : const casacore::String& refant)
3095 : {
3096 :
3097 0 : logSink() << LogOrigin("Calibrater","reRefant") << LogIO::NORMAL;
3098 :
3099 : // logSink() << "Beginning smoothing/interpolating method." << LogIO::POST;
3100 :
3101 :
3102 : // A pointer to an SVC
3103 0 : SolvableVisJones *svj(NULL);
3104 :
3105 : try {
3106 :
3107 : // Handle no in file
3108 0 : if (infile=="")
3109 0 : throw(AipsError("Please specify an input calibration table."));
3110 :
3111 : // Handle bad refantmode
3112 0 : if (refantmode!="strict" &&
3113 0 : refantmode!="flex")
3114 0 : throw(AipsError("Unrecognized refantmode!"));
3115 :
3116 : // Handle no outfile
3117 0 : if (outfile=="") {
3118 0 : outfile=infile;
3119 : logSink() << "Will overwrite input file with rerefant result."
3120 0 : << LogIO::POST;
3121 : }
3122 :
3123 :
3124 0 : svj = (SolvableVisJones*) createSolvableVisCal(calTableType(infile),*msmc_p);
3125 :
3126 : // Fill calibration table using setApply
3127 0 : RecordDesc applyparDesc;
3128 0 : applyparDesc.addField ("table", TpString);
3129 0 : Record applypar(applyparDesc);
3130 0 : applypar.define ("table", infile);
3131 0 : svj->setApply(applypar);
3132 :
3133 : // Do the work
3134 0 : svj->refantmode() = refantmode;
3135 0 : svj->refantlist().reference(getRefantIdxList(refant)); // replaces the default list
3136 0 : svj->applyRefAnt();
3137 :
3138 : // Store the result on disk
3139 0 : logSink() << "Storing result in " << outfile << LogIO::POST;
3140 :
3141 0 : if (outfile != "")
3142 0 : svj->calTableName()=outfile;
3143 0 : svj->storeNCT();
3144 :
3145 : // Clean up
3146 0 : if (svj) delete svj; svj=NULL;
3147 :
3148 : // Apparently, it worked
3149 0 : return true;
3150 :
3151 0 : } catch (AipsError x) {
3152 :
3153 : logSink() << LogIO::SEVERE
3154 : << "Caught Exception: "
3155 : << x.getMesg()
3156 0 : << LogIO::POST;
3157 : // Clean up
3158 0 : if (svj) delete svj; svj=NULL;
3159 :
3160 0 : throw(AipsError("Error in Calibrater::reRefant."));
3161 :
3162 : return false;
3163 0 : }
3164 : return false;
3165 : }
3166 :
3167 :
3168 : // List a calibration table
3169 0 : Bool Calibrater::listCal(const String& infile,
3170 : const String& field,
3171 : const String& antenna,
3172 : const String& spw,
3173 : const String& listfile,
3174 : const Int& pagerows) {
3175 :
3176 0 : SolvableVisCal *svc(NULL);
3177 0 : logSink() << LogOrigin("Calibrater","listCal");
3178 :
3179 : try {
3180 :
3181 : // Trap (currently) unsupported types
3182 0 : if (upcase(calTableType(infile))=="GSPLINE" ||
3183 0 : upcase(calTableType(infile))=="BPOLY")
3184 0 : throw(AipsError("GSPLINE and BPOLY tables cannot currently be listed."));
3185 :
3186 : // Get user's selected fields, ants
3187 0 : Vector<Int> ufldids=getFieldIdx(field);
3188 0 : Vector<Int> uantids=getAntIdx(antenna);
3189 :
3190 0 : String newSpw = spw;
3191 0 : Bool defaultSelect = false;
3192 0 : if (spw.empty()) { // list all channels (default)
3193 0 : defaultSelect = true;
3194 0 : newSpw = "*";
3195 : logSink() << LogIO::NORMAL1 << "Spws selected: ALL" << endl
3196 0 : << "Channels selected: ALL" << LogIO::POST;
3197 : }
3198 : // Get user's selected spw and channels
3199 0 : Vector<Int> uspwids=getSpwIdx(newSpw);
3200 0 : Matrix<Int> uchanids=getChanIdx(newSpw);
3201 0 : if (!defaultSelect) {
3202 : logSink() << LogIO::NORMAL1 << "Spw and Channel selection matrix: "
3203 : << endl << "Each rows shows: [ Spw , Start Chan , Stop Chan , Chan Step ]"
3204 0 : << endl << uchanids << LogIO::POST;
3205 : }
3206 : logSink() << LogIO::DEBUG2
3207 : << "uspwids = " << uspwids << endl
3208 0 : << "uchanids = " << uchanids << LogIO::POST;
3209 :
3210 : // By default, do first spw, first chan
3211 0 : if (uspwids.nelements()==0) {
3212 0 : uchanids.resize(1,4);
3213 0 : uchanids=0;
3214 : }
3215 :
3216 : // Set record format for calibration table application information
3217 0 : RecordDesc applyparDesc;
3218 0 : applyparDesc.addField ("table", TpString);
3219 :
3220 : // Create record with the requisite field values
3221 0 : Record applypar(applyparDesc);
3222 0 : applypar.define ("table", infile);
3223 :
3224 : // Generate the VisCal to be listed
3225 0 : svc = createSolvableVisCal(calTableType(infile),*msmc_p);
3226 0 : svc->setApply(applypar);
3227 :
3228 : // list it
3229 0 : svc->listCal(ufldids,uantids,uchanids, //uchanids(0,0),uchanids(0,1),
3230 : listfile,pagerows);
3231 :
3232 0 : if (svc) delete svc; svc=NULL;
3233 :
3234 0 : return true;
3235 :
3236 0 : } catch (AipsError x) {
3237 :
3238 : logSink() << LogIO::SEVERE
3239 : << "Caught Exception: "
3240 : << x.getMesg()
3241 0 : << LogIO::POST;
3242 : // Clean up
3243 0 : if (svc) delete svc; svc=NULL;
3244 :
3245 0 : throw(AipsError("Error in Calibrater::listCal."));
3246 :
3247 : return false;
3248 0 : }
3249 : return false;
3250 :
3251 : }
3252 :
3253 :
3254 :
3255 0 : Bool Calibrater::updateCalTable(const String& caltable) {
3256 :
3257 : // Call the SVC method that knows how
3258 0 : return NewCalTable::CTBackCompat(caltable);
3259 :
3260 : }
3261 :
3262 32 : void Calibrater::selectChannel(const String& spw) {
3263 :
3264 32 : if (mss_p && mssel_p) {
3265 :
3266 : // Refresh the frequencySelections object to feed to VI2, if relevant
3267 32 : frequencySelections_p.reset(new vi::FrequencySelections());
3268 :
3269 32 : vi::FrequencySelectionUsingChannels usingChannels;
3270 32 : usingChannels.add(*mss_p,mssel_p);
3271 32 : frequencySelections_p->add(usingChannels);
3272 :
3273 : // cout << usingChannels.toString() << endl;
3274 : // cout << "FS.size() = " << frequencySelections_p->size() << endl;
3275 :
3276 32 : }
3277 :
3278 : // TBD: Does frequencySelections_p support string info for reporting to logger?
3279 :
3280 32 : Matrix<Int> chansel = getChanIdx(spw);
3281 32 : uInt nselspw=chansel.nrow();
3282 :
3283 32 : if (nselspw==0)
3284 : logSink() << "Frequency selection: Selecting all channels in all spws."
3285 19 : << LogIO::POST;
3286 : else {
3287 :
3288 13 : logSink() << "Frequency selection: " << LogIO::POST;
3289 :
3290 : // Trap non-unit step (for now)
3291 13 : if (ntrue(chansel.column(3)==1)!=nselspw) {
3292 : logSink() << LogIO::WARN
3293 : << "Calibration does not support non-unit channel stepping; "
3294 : << "using step=1."
3295 0 : << LogIO::POST;
3296 0 : chansel.column(3)=1;
3297 : }
3298 :
3299 13 : Vector<Int> uspw(chansel.column(0));
3300 13 : Vector<Int> ustart(chansel.column(1));
3301 13 : Vector<Int> uend(chansel.column(2));
3302 13 : Vector<Int> ustep(chansel.column(3));
3303 :
3304 13 : logSink() << LogIO::NORMAL;
3305 47 : for (uInt i=0;i<nselspw;++i) {
3306 :
3307 68 : logSink() << ". Spw " << uspw(i) << ":"
3308 136 : << ustart(i) << "~" << uend(i)
3309 34 : << " (" << uend(i)-ustart(i)+1 << " channels,"
3310 68 : << " step by " << ustep(i) << ")"
3311 68 : << endl;
3312 :
3313 : } // i
3314 13 : logSink() << LogIO::POST;
3315 :
3316 13 : } // non-triv spw selection
3317 :
3318 32 : }
3319 :
3320 :
3321 242 : Bool Calibrater::cleanup() {
3322 :
3323 : // logSink() << LogOrigin("Calibrater","cleanup") << LogIO::NORMAL;
3324 :
3325 : // Delete the VisCals
3326 242 : reset();
3327 :
3328 : // Delete derived dataset stuff
3329 242 : if(mssel_p) delete mssel_p; mssel_p=0;
3330 242 : if(mss_p) delete mss_p; mss_p=0;
3331 242 : frequencySelections_p.reset();
3332 :
3333 : // Delete the current VisEquation
3334 242 : if(ve_p) delete ve_p; ve_p=0;
3335 :
3336 242 : return true;
3337 :
3338 : }
3339 :
3340 : // Parse refant specification
3341 32 : Vector<Int> Calibrater::getRefantIdxList(const String& refant) {
3342 :
3343 32 : Vector<Int> irefant;
3344 32 : if (refant.length()==0) {
3345 : // Nothing specified, return -1 in single-element vector
3346 0 : irefant.resize(1);
3347 0 : irefant(0)=-1;
3348 : }
3349 : else {
3350 : // parse the specification
3351 32 : MSSelection msselect;
3352 32 : msselect.setAntennaExpr(refant);
3353 32 : Vector<Int> iref=msselect.getAntenna1List(mssel_p);
3354 32 : if (anyLT(iref,0))
3355 0 : cout << "Negated selection (via '!') not yet implemented for refant," << endl << " and will be ignored." << endl;
3356 32 : irefant=iref(iref>-1).getCompressedArray();
3357 32 : if (irefant.nelements()==0) {
3358 0 : irefant.resize(1);
3359 0 : irefant(0)=-1;
3360 : }
3361 32 : }
3362 :
3363 32 : return irefant;
3364 0 : }
3365 :
3366 :
3367 : // Interpret refant *index*
3368 0 : Vector<Int> Calibrater::getAntIdx(const String& antenna) {
3369 :
3370 0 : MSSelection msselect;
3371 0 : msselect.setAntennaExpr(antenna);
3372 0 : return msselect.getAntenna1List(mssel_p);
3373 :
3374 0 : }
3375 :
3376 : // Interpret field indices (MSSelection)
3377 41 : Vector<Int> Calibrater::getFieldIdx(const String& fields) {
3378 :
3379 41 : MSSelection mssel;
3380 41 : mssel.setFieldExpr(fields);
3381 82 : return mssel.getFieldList(mssel_p);
3382 :
3383 41 : }
3384 :
3385 : // Interpret spw indices (MSSelection)
3386 7 : Vector<Int> Calibrater::getSpwIdx(const String& spws) {
3387 :
3388 7 : MSSelection mssel;
3389 7 : mssel.setSpwExpr(spws);
3390 :
3391 : /*
3392 : cout << "mssel.getSpwList(mssel_p) = " << mssel.getSpwList(mssel_p) << endl;
3393 : cout << "mssel.getChanList(mssel_p) = " << mssel.getChanList(mssel_p) << endl; cout << "Vector<Int>() = " << Vector<Int>() << endl;
3394 : */
3395 :
3396 : // Use getChanList (column 0 for spw) because it is
3397 : // more reliable about the number and order of specified spws.
3398 : // return mssel.getSpwList(mssel_p);
3399 7 : Matrix<Int> chanmat=mssel.getChanList(mssel_p);
3400 7 : if (chanmat.nelements()>0)
3401 0 : return chanmat.column(0);
3402 : else
3403 7 : return Vector<Int>();
3404 :
3405 7 : }
3406 :
3407 32 : Matrix<Int> Calibrater::getChanIdx(const String& spw) {
3408 :
3409 32 : MSSelection mssel;
3410 32 : mssel.setSpwExpr(spw);
3411 :
3412 64 : return mssel.getChanList(mssel_p);
3413 :
3414 32 : }
3415 :
3416 :
3417 : // Query apply types if we are calibrating the weights
3418 0 : Bool Calibrater::calWt() {
3419 :
3420 0 : Int napp(vc_p.nelements());
3421 : // Return true as soon as we find a type which is cal'ing wts
3422 0 : for (Int iapp=0;iapp<napp;++iapp)
3423 0 : if (vc_p[iapp] && vc_p[iapp]->calWt())
3424 0 : return true;
3425 :
3426 : // None cal'd weights, so return false
3427 0 : return false;
3428 :
3429 : }
3430 :
3431 71 : Bool Calibrater::ok() {
3432 :
3433 71 : if ((simdata_p ||
3434 71 : (ms_p && mssel_p))
3435 71 : && ve_p) {
3436 71 : return true;
3437 : }
3438 : else {
3439 0 : logSink() << "Calibrater is not yet initialized" << LogIO::POST;
3440 0 : return false;
3441 : }
3442 : }
3443 : // The standard solving mechanism
3444 32 : casacore::Bool Calibrater::genericGatherAndSolve()
3445 : {
3446 :
3447 : #ifdef _OPENMP
3448 32 : Double Tsetup(0.0),Tgather(0.0),Tsolve(0.0),Tadd(0.0);
3449 32 : Double time0=omp_get_wtime();
3450 : #endif
3451 :
3452 : // Create record to store results
3453 32 : resRec_ = Record();
3454 32 : Record attemptRec;
3455 32 : attemptRec = Record();
3456 :
3457 : // Condition solint values
3458 32 : svc_p->reParseSolintForVI2();
3459 :
3460 : // Organize VI2 layers for solving:
3461 32 : CalSolVi2Organizer vi2org;
3462 :
3463 : //-------------------------------------------------
3464 : // NB: Populating the vi2org should probably be delegated to the svc_p
3465 : // since that is where most of the required (non-data) info is,
3466 : // and then some calibration types may introduce layers that are
3467 : // very specific to them (e.g., SD data selection, etc.)
3468 : // e.g., replace the following with something like:
3469 : //
3470 : // svc_p->formVi2LayerFactories(vi2org,mssel_p,ve_p);
3471 : // ssvp_p // a simdata version
3472 : //
3473 :
3474 : // Add the (bottom) data access layer
3475 32 : if (simdata_p)
3476 : // Simulated data (used for testing)
3477 0 : vi2org.addSimIO(ssvp_p);
3478 : else
3479 : {
3480 : // Real (selected) data in an MS (channel selection handled by addDiskIO method)
3481 : // The iteration time-interval is the solution interval
3482 32 : vi2org.addDiskIO(mssel_p,svc_p->solTimeInterval(),
3483 32 : svc_p->combobs(),svc_p->combscan(),
3484 32 : svc_p->combfld(),svc_p->combspw(),
3485 : true, // use MSIter2
3486 32 : frequencySelections_p); // Tell VI2 factory about the freq selection!
3487 : }
3488 :
3489 : // Add ad hoc SD section layer (e.g., OTF select of raster boundaries, etc.)
3490 : // only double circle gain calibration is implemented
3491 32 : bool SD = svc_p->longTypeName().startsWith("SDGAIN_OTFD");
3492 32 : if (SD) {
3493 0 : vi2org.addCalFilter(calFilterConfig_p);
3494 : }
3495 :
3496 : // Add pre-cal layer, using the VisEquation
3497 : // (include control for corr-dep flags)
3498 32 : vi2org.addCalForSolving(*ve_p,corrDepFlags_);
3499 :
3500 :
3501 : // Add the freq-averaging layer, if needed
3502 : //cout << "svc_p->fintervalChV() = " << svc_p->fintervalChV() << endl;
3503 : //cout << "svc_p->fsolint() = " << svc_p->fsolint() << endl;
3504 : // TBD: improve the following logic with new method in SVC...
3505 64 : if (!svc_p->freqDepMat() || // entirely unchannelized cal OR
3506 32 : (svc_p->freqDepPar() && // channelized par and
3507 0 : svc_p->fsolint()!="none" && // some partial channel averaging
3508 32 : anyGT(svc_p->fintervalChV(),1.0)) // explicity specified (non-trivially)
3509 : ) {
3510 0 : Vector<Int> chanbin(svc_p->fintervalChV().nelements());
3511 0 : convertArray(chanbin,svc_p->fintervalChV());
3512 0 : vi2org.addChanAve(chanbin);
3513 0 : }
3514 :
3515 : // Add the time-averaging layer, if needed
3516 : // NB: There is some evidence that doing time-averaging _after_
3517 : // channel averaging may be more expensive then doing it _before_,
3518 : // but this ensures that the meta-info (e.g. time, timeCentroid, etc.)
3519 : // appear correctly in the VB2 accessed at this scope.
3520 : // The problem is that AveragingTvi2 has not explicitly
3521 : // implemented all of the relevant data accessors; it just
3522 : // assumes---incorrectly---that you are going to use its VB2.
3523 : // Making AveragingTvi2 the top layer ensures you do use it.
3524 : // (Hard-wired use of AveragingTvi2 in MSTransform set it up
3525 : // as the top [last] layer, and has been well-tested...)
3526 :
3527 : //cout << "svc_p->solint() = " << svc_p->solint() << endl;
3528 : //cout << "svc_p->solTimeInterval() = " << svc_p->solTimeInterval() << endl;
3529 : //cout << "svc_p->preavg() = " << svc_p->preavg() << endl;
3530 64 : if (!svc_p->solint().contains("int") && // i.e., not "int")
3531 32 : svc_p->preavg() != -DBL_MAX) {
3532 0 : Float avetime(svc_p->solTimeInterval()); // avetime is solint, nominally
3533 : // Use preavg instead, if...
3534 0 : if (svc_p->preavg()>FLT_EPSILON && // ...set meaningfully
3535 0 : svc_p->preavg()<svc_p->solTimeInterval()) // ...and less than solint
3536 0 : avetime=svc_p->preavg();
3537 0 : vi2org.addTimeAve(avetime); // use min of solint and preavg here!
3538 : }
3539 : // small@jive.eu (2021-12-17): Currently we only handle the corrcomb
3540 : // cases of "all" and "none"; there will be a need to extend this, but the
3541 : // PolAverageTVILayerFactory that underlies this feature will also need to be extended
3542 : // to make that possible
3543 : // gmoellen (2024-08-06): We now support 'stokes' (forms Stokes I) and 'parallel' (direct
3544 : // statistical combo of parallel hands). If one parallel-hand is flagged, 'stokes'
3545 : // effectively flags both; 'parallel' will retain the unflagged parallel hand.
3546 61 : if (svc_p->corrcomb().contains("stokes") ||
3547 29 : svc_p->corrcomb().contains("parallel") ) {
3548 : //cerr << "Calibrater::genericGatherAndSolve(): Combining correlations: "
3549 : // << svc_p->corrcomb()
3550 : // << endl;
3551 5 : vi2org.addCorrCombine(svc_p->corrcomb());
3552 : }
3553 : //else {
3554 : // cerr << "Calibrater::genericGatherAndSolve(): Not combining correlations!" << endl;
3555 : //}
3556 :
3557 : // vi2org should be fully configured at this point
3558 : //-------------------------------------------------
3559 :
3560 :
3561 : // Form the VI2 to drive data iteration below
3562 32 : vi::VisibilityIterator2& vi(vi2org.makeFullVI());
3563 : // cout << "VI Layers: " << vi.ViiType() << endl;
3564 :
3565 : // Establish output freq meta data prior to solving (CAS-12735 and related)
3566 :
3567 : // Discern (even implicitly-)selected spw subset
3568 32 : set<uInt> selspwset;
3569 32 : if (!simdata_p) {
3570 :
3571 : // We can't proceed rationally if we have more than one freqSelection (multiple MSs)
3572 : // (this is a sort-of out-of-the-blue test here....)
3573 32 : AlwaysAssert(frequencySelections_p->size()<2,AipsError);
3574 :
3575 : // TBD can we use the msmc_p? (or maybe it is for the whole MS, not just the selected one?)
3576 32 : MSMetaData msmdtmp(mssel_p,50.0);
3577 32 : selspwset=msmdtmp.getSpwIDs();
3578 32 : }
3579 :
3580 : // Extract selected spw list as a Vector<uInt>
3581 32 : Vector<uInt> selspwlist;
3582 32 : if (!simdata_p && selspwset.size()>0) {
3583 32 : Vector<uInt> vtmp(selspwset.begin(),selspwset.size(),0);
3584 32 : selspwlist.reference(vtmp);
3585 32 : } else {
3586 : // All spws in the MS (this is a last resort!)
3587 : // NB: This is used in simdata_p=True context!
3588 0 : selspwlist.resize(msmc_p->nSpw());
3589 0 : indgen(selspwlist);
3590 : }
3591 : // cout << "selspwlist = " << selspwlist << endl;
3592 :
3593 : // Delegate freq meta calculation to the SVC (wherein this info is stored)
3594 32 : svc_p->discernAndSetSolnFrequencies(vi,selspwlist);
3595 :
3596 :
3597 : // Access to the net VB2 for forming SolveDataBuffers in the SDBList below
3598 : // NB: this assumes that the internal address of the VI2's VB2 will never change!
3599 31 : vi::VisBuffer2 *vb = vi.getImpl()->getVisBuffer();
3600 :
3601 : // VI2 should now be ready....
3602 :
3603 : // Discern number of Solutions and how many VI2 chunks per each
3604 31 : Int nSol(0);
3605 31 : Vector<Int> nChPSol(1000,0); // nominal size; will enlarge as needed, and trim at end
3606 :
3607 : {
3608 : //cout << "Counting chunks and solutions..." << endl;
3609 31 : uInt isol(0);
3610 : // uInt ichk(0);
3611 268 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
3612 : // NB: this loop NEVER touches the the subChunks, since
3613 : // this would cause data I/O, cal, averaging, etc.
3614 :
3615 : // int frame=vi.getImpl()->getReportingFrameOfReference();
3616 : // cout << "ichk=" << ichk << " freqs= " << vi.getImpl()->getFrequencies(-1.0e0,frame,selspwlist(ispw),0);
3617 :
3618 :
3619 : // Enlarge nChPSol Vector, if necessary
3620 237 : if (isol>nChPSol.nelements()-1) {
3621 : //cout << " ...At isol=" << isol << ", resizing nchPSol: " << nChPSol.nelements() << "-->";
3622 0 : nChPSol.resize(isol*2,True); // double length, copy existing elements
3623 0 : nChPSol(Slice(isol,isol,1)).set(0); // init new elements
3624 : //cout << nChPSol.nelements() << endl;
3625 : //cout << "sol counting: " << isol << " " << nChPSol << " (" << ntrue(nChPSol>0) << "/" << nChPSol.nelements() << ")" << endl;
3626 : }
3627 :
3628 : // incr chunk counter for current solution
3629 237 : ++nChPSol(isol);
3630 : //cout << "sol counting: " << isol << " " << nChPSol(Slice(0,min(isol+2,nChPSol.nelements()),1))
3631 : // << " (" << ntrue(nChPSol>0) << "/" << nChPSol.nelements() << ")"
3632 : // << " keyCh=" << vi.keyChange()
3633 : // << endl;
3634 :
3635 : // next chunk is nominally next solution
3636 237 : ++isol;
3637 :
3638 : // If combining spw or field, and the next chunk changes by
3639 : // DATA_DESC_ID or FIELD_ID (and _not_ TIME, which is
3640 : // changing slower than these when combine is on),
3641 : // then the next chunk should be counted in the same soln
3642 699 : if ( (svc_p->combspw() && vi.keyChange()=="DATA_DESC_ID") ||
3643 462 : (svc_p->combfld() && vi.keyChange()=="FIELD_ID") )
3644 12 : --isol; // next one is the same solution interval
3645 :
3646 : }
3647 : // Trim list to realized length, which is the number of solutions we'll
3648 : // attempt below
3649 31 : nSol=ntrue(nChPSol>0);
3650 31 : nChPSol.resize(nSol,true);
3651 : //cout << "nSol =" << nSol << endl;
3652 : //cout << "nChPSol=" << nChPSol << endl;
3653 :
3654 : /* from SVC...
3655 : if (svc_p->combobs())
3656 : logSink() << "Combining observation Ids." << LogIO::POST;
3657 : if (svc_p->combscan())
3658 : logSink() << "Combining scans." << LogIO::POST;
3659 : if (svc_p->combspw())
3660 : logSink() << "Combining spws: " << spwlist << " -> " << spwlab << LogIO::POST;
3661 : if (svc_p->combfld())
3662 : logSink() << "Combining fields." << LogIO::POST;
3663 : */
3664 :
3665 31 : logSink() << "For solint = " << svc_p->solint() << ", found "
3666 : << nSol << " solution intervals."
3667 31 : << LogIO::POST;
3668 : }
3669 :
3670 : // throw(AipsError("EARLY ESCAPE!!"));
3671 :
3672 : // Create the output caltable
3673 : // (this version doesn't revise frequencies)
3674 31 : svc_p->createMemCalTable2();
3675 :
3676 31 : Vector<Float> spwwts(msmc_p->nSpw(),-1.0);
3677 31 : Vector<Int64> nexp(msmc_p->nSpw(),0), natt(msmc_p->nSpw(),0),nsuc(msmc_p->nSpw(),0), nbelowsnr(msmc_p->nSpw(),0), nbelowbl(msmc_p->nSpw(),0);
3678 :
3679 : // Start Collecting counts
3680 31 : CalCounts* calCounts = new CalCounts();
3681 31 : calCounts->initCounts(msmc_p->nSpw(),msmc_p->nAnt(), svc_p->nPar());
3682 : // Set up results record
3683 31 : std::map<Int, std::map<String, Vector<Int>>> resultMap;
3684 :
3685 :
3686 : #ifdef _OPENMP
3687 31 : Tsetup+=(omp_get_wtime()-time0);
3688 : #endif
3689 :
3690 31 : Int nGood(0);
3691 31 : vi.originChunks();
3692 31 : Int nGlobalChunks=0; // counts VI chunks globally
3693 256 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
3694 :
3695 : #ifdef _OPENMP
3696 225 : time0=omp_get_wtime();
3697 : #endif
3698 :
3699 : // Data will accumulate here
3700 225 : SDBList sdbs;
3701 :
3702 : // Gather the chunks/VBs for this solution
3703 : // Solution boundaries will ALWAYS occur on chunk boundaries,
3704 : // though some chunk boundaries will be ignored in the
3705 : // combine='spw' or 'field' context
3706 :
3707 225 : for (Int ichunk=0; // count chunks in this solution
3708 462 : ichunk<nChPSol(isol)&&vi.moreChunks(); // while more chunks needed _and_ available
3709 237 : ++ichunk,vi.nextChunk()) { // advance to next chunk
3710 :
3711 : // Global chunk counter
3712 237 : ++nGlobalChunks;
3713 :
3714 : // Loop over VB2s in this chunk
3715 : // (we get more then one when preavg<solint)
3716 237 : Int ivb(0);
3717 237 : for (vi.origin();
3718 27756 : vi.more();
3719 27519 : ++ivb,vi.next()) {
3720 :
3721 : // Add this VB to the SDBList
3722 : #ifdef _OPENMP
3723 27519 : Double Tadd0=omp_get_wtime();
3724 : #endif
3725 :
3726 27519 : sdbs.add(*vb);
3727 :
3728 : #ifdef _OPENMP
3729 27519 : Tadd+=(omp_get_wtime()-Tadd0);
3730 : #endif
3731 :
3732 : // Keep track of spws seen but not included in solving
3733 27519 : Int ispw=vb->spectralWindows()(0);
3734 27519 : if (spwwts(ispw)<0) spwwts(ispw)=0.0f;
3735 27519 : spwwts(ispw)+=sum(vb->weightSpectrum());
3736 :
3737 : /*
3738 : Double modv(86400.0);
3739 : cout.precision(7);
3740 : cout << "isol=" << isol
3741 : << " ichk="<<ichunk
3742 : << " ivb="<<ivb
3743 : << " sc="<<vb->scan()(0)
3744 : << " fld="<<vb->fieldId()(0)
3745 : << " spw="<<vb->spectralWindows()(0)
3746 : << " nr="<<vb->nRows()
3747 : << " nch="<<vb->nChannels() << flush;
3748 : cout << " f="<<mean(vb->getFrequencies(0))/1e9
3749 : << "->"<< mean(sdbs.freqs())/1e9
3750 : << " t="<<fmod(mean(vb->time()),modv)
3751 : << "->"<< fmod(sdbs.aggregateTime(),modv)
3752 : << " tC="<<fmod(mean(vb->timeCentroid()),modv)
3753 : << "->"<< fmod(sdbs.aggregateTimeCentroid(),modv)
3754 : << " wt="<<sum(vb->weightSpectrum())
3755 : << " nSDB=" << sdbs.nSDB()
3756 : << " nGlChks=" << nGlobalChunks
3757 : << " (keych=" << vi.keyChange() << ")"
3758 : << endl;
3759 : //*/
3760 :
3761 : } // VI2 subchunks (VB2s)
3762 : } // VI2 chunks
3763 :
3764 : // Which spw is this?
3765 225 : Int thisSpw(sdbs.aggregateSpw());
3766 :
3767 : // Expecting a solution
3768 225 : nexp(thisSpw)+=1;
3769 : // Initialize the antennaMap_ in svc
3770 225 : svc_p->clearMap();
3771 : // Get expected and data unflagged accumulation
3772 225 : svc_p->expectedUnflagged(sdbs);
3773 :
3774 :
3775 : #ifdef _OPENMP
3776 225 : Tgather+=(omp_get_wtime()-time0);
3777 225 : time0=omp_get_wtime();
3778 : #endif
3779 :
3780 :
3781 225 : if (sdbs.Ok()) {
3782 :
3783 : // Some unflagged data, so Attempting a solution
3784 225 : natt(thisSpw)+=1;
3785 :
3786 : // make phase- or amp-only, if necessary
3787 225 : sdbs.enforceAPonData(svc_p->apmode());
3788 :
3789 : // zero cross-hand weights, if necessary
3790 225 : sdbs.enforceSolveWeights(svc_p->phandonly());
3791 :
3792 : // Synchronize meta-info in SVC
3793 225 : svc_p->syncSolveMeta(sdbs);
3794 :
3795 : // Set or verify freqs in the caltable
3796 : //svc_p->setOrVerifyCTFrequencies(thisSpw);
3797 225 : svc_p->setCTFrequencies(thisSpw);
3798 :
3799 : // Size the solvePar arrays inside SVC
3800 : // (smart: if freqDepPar()=F, uses 1)
3801 : // returns the number of channel solutions to iterate over
3802 : //Int nChanSol=svc_p->sizeSolveParCurrSpw(sdbs.nChannels());
3803 225 : Int nChanSol=svc_p->sizeSolveParCurrSpw((svc_p->freqDepPar() ? sdbs.nChannels() : 1));
3804 :
3805 225 : if (svc_p->useGenericSolveOne()) {
3806 :
3807 : // We'll use the generic solver
3808 0 : VisCalSolver2 vcs(svc_p->solmode(),svc_p->rmsthresh());
3809 :
3810 : // Guess from the data
3811 0 : svc_p->guessPar(sdbs,corrDepFlags_);
3812 :
3813 0 : Bool totalGoodSol(False); // Will be set True if any channel is good
3814 : //for (Int ich=0;ich<nChanSol;++ich) {
3815 0 : for (Int ich=nChanSol-1;ich>-1;--ich) {
3816 0 : svc_p->markTimer();
3817 0 : svc_p->focusChan()=ich;
3818 :
3819 : // Execute the solve
3820 0 : Bool goodsol=vcs.solve(*ve_p,*svc_p,sdbs);
3821 :
3822 0 : if (goodsol) {
3823 0 : totalGoodSol=True;
3824 :
3825 0 : svc_p->formSolveSNR();
3826 0 : svc_p->applySNRThreshold();
3827 : }
3828 : else
3829 0 : svc_p->currMetaNote();
3830 :
3831 : // Record solution in its channel, good or bad
3832 0 : if (svc_p->freqDepPar())
3833 0 : svc_p->keep1(ich);
3834 :
3835 : } // ich
3836 :
3837 0 : if (totalGoodSol) {
3838 : // Keep this good solution, and count it
3839 0 : svc_p->keepNCT();
3840 0 : ++nGood;
3841 0 : nsuc(thisSpw)+=1;
3842 : }
3843 :
3844 0 : } // useGenericSolveOne
3845 : else {
3846 : // Use self-directed individual solve
3847 : // TBD: return T/F for goodness?
3848 : //cout << "Calling selfSolveOne()!!!!!!" << endl;
3849 225 : svc_p->selfSolveOne(sdbs);
3850 :
3851 : // Keep this good solution, and count it
3852 225 : svc_p->keepNCT();
3853 225 : ++nGood;
3854 225 : nsuc(thisSpw)+=1;
3855 : }
3856 :
3857 : } // sdbs.Ok()
3858 : else {
3859 : // Synchronize meta-info in SVC
3860 0 : svc_p->syncSolveMeta(sdbs);
3861 0 : cout << "Found no unflagged data at:";
3862 0 : svc_p->currMetaNote();
3863 : }
3864 : //cout << endl;
3865 : // Get all the antenna value counts
3866 225 : resultMap = svc_p->getAntennaMap();
3867 225 : calCounts->addAntennaCounts(thisSpw,msmc_p->nAnt(), svc_p->nPar(), resultMap);
3868 :
3869 : #ifdef _OPENMP
3870 225 : Tsolve+=(omp_get_wtime()-time0);
3871 : #endif
3872 :
3873 : // throw(AipsError("EARLY ESCAPE!!"));
3874 :
3875 225 : } // isol
3876 :
3877 : // Report nGood to logger
3878 : logSink() << " Found good "
3879 31 : << svc_p->typeName() << " solutions in "
3880 : << nGood << " solution intervals."
3881 31 : << LogIO::POST;
3882 :
3883 : #ifdef _OPENMP
3884 : #ifdef REPORT_CAL_TIMING
3885 : cout << "Calibrater::genericGatherAndSolve Timing: " << endl;
3886 : cout << " setup=" << Tsetup
3887 : << " gather=" << Tgather
3888 : << " (SDBadd=" << Tadd << ")"
3889 : << " solve=" << Tsolve
3890 : << " total=" << Tsetup+Tgather+Tsolve
3891 : // << " tick=" << omp_get_wtick()
3892 : << endl;
3893 : #endif
3894 : #endif
3895 :
3896 : // Report spws that were seen but not solved
3897 31 : Vector<Bool> unsolspw=(spwwts==0.0f);
3898 31 : summarize_uncalspws(unsolspw, "solv");
3899 :
3900 : // throw(AipsError("EARLY ESCAPE!!"));
3901 :
3902 31 : if (nGood>0) {
3903 31 : if (svc_p->typeName()!="BPOLY") { // needed?
3904 : // apply refant, etc.
3905 31 : svc_p->clearRefantMap();
3906 31 : svc_p->globalPostSolveTinker();
3907 :
3908 : // write to disk
3909 31 : svc_p->storeNCT();
3910 : }
3911 : }
3912 : else {
3913 : logSink() << "No output calibration table written."
3914 0 : << LogIO::POST;
3915 : }
3916 :
3917 : // Collect and update the refants
3918 31 : std::map<Int, std::map<Int, Int>> refantMap;
3919 31 : refantMap = svc_p->getRefantMap();
3920 31 : calCounts->updateRefants(msmc_p->nSpw(), msmc_p->nAnt(), refantMap);
3921 : // Compile all accumulated counts into a record
3922 31 : resRec_ = calCounts->makeRecord(msmc_p->nAnt(), vi.nPolarizationIds());
3923 :
3924 : // print a matrix to the logger
3925 31 : logSink() << "----- PER ANTENNA INFO -----" << LogIO::POST;
3926 : // print the antenna list
3927 31 : logSink() << " ";
3928 347 : for (Int ant=0; ant<msmc_p->nAnt(); ant++) {
3929 316 : logSink() << " ANT: " << ant << " ";
3930 : }
3931 31 : logSink() << LogIO::POST;
3932 :
3933 155 : for (Int spw=0; spw < msmc_p->nSpw(); spw++) {
3934 :
3935 124 : logSink() << LogIO::NORMAL << "SPW: " << spw;
3936 :
3937 1388 : for (Int ant=0; ant < msmc_p->nAnt(); ant++) {
3938 1264 : logSink() << " " << calCounts->antMapVal(spw, ant, "above_minsnr") << " ";
3939 : }
3940 124 : logSink() << LogIO::POST;
3941 : }
3942 :
3943 31 : logSink() << "----- PER SPW INFO -----" << LogIO::POST;
3944 : // print the result fields
3945 31 : logSink() << " ";
3946 31 : logSink() << " expected data_unflagged above_minblperant above_minsnr";
3947 31 : logSink() << LogIO::POST;
3948 155 : for (Int spw = 0; spw < msmc_p->nSpw(); spw++) {
3949 124 : logSink() << LogIO::NORMAL << "SPW: " << spw << " ";
3950 124 : logSink() << calCounts->spwMapVal(spw, "expected") << " ";
3951 124 : logSink() << calCounts->spwMapVal(spw, "data_unflagged") << " ";
3952 124 : logSink() << calCounts->spwMapVal(spw, "above_minblperant") << " ";
3953 124 : logSink() << calCounts->spwMapVal(spw, "above_minsnr");
3954 124 : logSink() << LogIO::POST;
3955 : }
3956 :
3957 31 : logSink() << "----- GLOBAL INFO -----" << LogIO::POST;
3958 31 : logSink() << "expected data_unflagged above_minblperant above_minsnr";
3959 31 : logSink() << LogIO::POST;
3960 31 : logSink() << calCounts->totalMapVal("expected") << " ";
3961 31 : logSink() << calCounts->totalMapVal("data_unflagged") << " ";
3962 31 : logSink() << calCounts->totalMapVal("above_minblperant") << " ";
3963 31 : logSink() << calCounts->totalMapVal("above_minsnr");
3964 31 : logSink() << LogIO::POST;
3965 :
3966 : // Fill activity record
3967 :
3968 31 : actRec_ = Record();
3969 31 : actRec_.define("origin","Calibrater::genericGatherAndSolve");
3970 31 : actRec_.define("nExpected",nexp);
3971 31 : actRec_.define("nAttempt",natt);
3972 31 : actRec_.define("nSucceed",nsuc);
3973 :
3974 : //cout << nexp << ", " << natt << ", " << nsuc << endl;
3975 :
3976 : {
3977 31 : Record solveRec=svc_p->solveActionRec();
3978 31 : if (solveRec.nfields()>0)
3979 0 : actRec_.merge(solveRec);
3980 31 : }
3981 :
3982 : // Reach here, all is good
3983 31 : return True;
3984 :
3985 35 : }
3986 :
3987 :
3988 0 : void Calibrater::writeHistory(LogIO& /*os*/, Bool /*cliCommand*/)
3989 : {
3990 : // Disabled 2016/04/19: avoid direct MS.HISTORY updates from
3991 : // below the python level, FOR NOW
3992 :
3993 0 : return;
3994 : /*
3995 : if (!historytab_p.isNull()) {
3996 : if (histLockCounter_p == 0) {
3997 : historytab_p.lock(false);
3998 : }
3999 : ++histLockCounter_p;
4000 :
4001 : os.postLocally();
4002 : if (cliCommand) {
4003 : hist_p->cliCommand(os);
4004 : } else {
4005 : hist_p->addMessage(os);
4006 : }
4007 :
4008 : if (histLockCounter_p == 1) {
4009 : historytab_p.unlock();
4010 : }
4011 : if (histLockCounter_p > 0) {
4012 : --histLockCounter_p;
4013 : }
4014 : } else {
4015 : os << LogIO::SEVERE << "calibrater is not yet initialized" << LogIO::POST;
4016 : }
4017 : */
4018 : }
4019 :
4020 32 : void Calibrater::setCalFilterConfiguration(String const &type,
4021 : Record const &config) {
4022 : // currently only SDDoubleCircleGainCal requires data filtering
4023 32 : if (type.startsWith("SDGAIN_OTFD")) {
4024 0 : calFilterConfig_p.define("mode", "SDGAIN_OTFD");
4025 0 : if (config.isDefined("smooth")) {
4026 0 : calFilterConfig_p.define("smooth", config.asBool("smooth"));
4027 : }
4028 0 : if (config.isDefined("radius")) {
4029 0 : calFilterConfig_p.define("radius", config.asString("radius"));
4030 : }
4031 : }
4032 32 : }
4033 :
4034 : // *********************************************
4035 : // OldCalibrater implementations that use vs_p
4036 :
4037 0 : OldCalibrater::OldCalibrater():
4038 : Calibrater(),
4039 0 : vs_p(0),
4040 0 : rawvs_p(0)
4041 : {
4042 : // cout << "This is the OLD VI2-aware Calibrater" << endl;
4043 0 : }
4044 :
4045 0 : OldCalibrater::OldCalibrater(String msname):
4046 : Calibrater(msname),
4047 0 : vs_p(0),
4048 0 : rawvs_p(0)
4049 : {
4050 0 : }
4051 :
4052 : /*
4053 : OldCalibrater::OldCalibrater(const OldCalibrater & other) :
4054 : Calibrater(other)
4055 : {
4056 : operator=(other);
4057 : }
4058 :
4059 : OldCalibrater& OldCalibrater::operator=(const OldCalibrater & other)
4060 : {
4061 : Calibrater::operator=(other); // copy parental units
4062 : vs_p=other.vs_p;
4063 : rawvs_p=other.rawvs_p;
4064 : return *this;
4065 : }
4066 : */
4067 :
4068 0 : OldCalibrater::~OldCalibrater()
4069 : {
4070 0 : OldCalibrater::cleanupVisSet();
4071 0 : }
4072 :
4073 :
4074 : // Select data (using MSSelection syntax)
4075 0 : void OldCalibrater::selectvis(const String& time,
4076 : const String& spw,
4077 : const String& scan,
4078 : const String& field,
4079 : const String& intent,
4080 : const String& obsIDs,
4081 : const String& baseline,
4082 : const String& uvrange,
4083 : const String& chanmode,
4084 : const Int& nchan,
4085 : const Int& start,
4086 : const Int& step,
4087 : const MRadialVelocity& mStart,
4088 : const MRadialVelocity& mStep,
4089 : const String& msSelect)
4090 : {
4091 : // Define primary measurement set selection criteria
4092 : // Inputs:
4093 : // time
4094 : // spw
4095 : // scan
4096 : // field
4097 : // intent
4098 : // obsIDs
4099 : // baseline
4100 : // uvrange
4101 : // chanmode const String& Frequency/velocity selection mode
4102 : // ("channel", "velocity" or
4103 : // "opticalvelocity")
4104 : // nchan const Int& No of channels to select
4105 : // start const Int& Start channel to select
4106 : // step const Int& Channel increment
4107 : // mStart const MRadialVelocity& Start radial vel. to select
4108 : // mStep const MRadialVelocity& Radial velocity increment
4109 : // msSelect const String& MS selection string (TAQL)
4110 : // Output to private data:
4111 : //
4112 0 : logSink() << LogOrigin("Calibrater","selectvis") << LogIO::NORMAL3;
4113 :
4114 : try {
4115 :
4116 : /*
4117 : cout << "time = " << time << " " << time.length() <<endl;
4118 : cout << "spw = " << spw << " " << spw.length() <<endl;
4119 : cout << "scan = " << scan << " " << scan.length() <<endl;
4120 : cout << "field = " << field << " " << field.length() <<endl;
4121 : cout << "baseline = " << baseline << " " << baseline.length() << endl;
4122 : cout << "uvrange = " << uvrange << " " << uvrange.length() << endl;
4123 : */
4124 :
4125 0 : logSink() << "Selecting data" << LogIO::POST;
4126 :
4127 : // Apply selection to the original MeasurementSet
4128 0 : logSink() << "Performing selection on MeasurementSet" << endl;
4129 :
4130 0 : if (mssel_p) {
4131 0 : delete mssel_p;
4132 0 : mssel_p=0;
4133 : };
4134 :
4135 : // Report non-trivial user selections
4136 0 : if (time!="")
4137 0 : logSink() << " Selecting on time: '" << time << "'" << endl;
4138 0 : if (spw!="")
4139 0 : logSink() << " Selecting on spw: '" << spw << "'" << endl;
4140 0 : if (scan!="")
4141 0 : logSink() << " Selecting on scan: '" << scan << "'" << endl;
4142 0 : if (field!="")
4143 0 : logSink() << " Selecting on field: '" << field << "'" << endl;
4144 0 : if (intent!="")
4145 0 : logSink() << " Selecting on intent: '" << intent << "'" << endl;
4146 0 : if(obsIDs != "")
4147 0 : logSink() << " Selecting by observation IDs: '" << obsIDs << "'" << endl;
4148 0 : if (baseline!="")
4149 0 : logSink() << " Selecting on antenna/baseline: '" << baseline << "'" << endl;
4150 0 : if (uvrange!="")
4151 0 : logSink() << " Selecting on uvrange: '" << uvrange << "'" << endl;
4152 0 : if (msSelect!="")
4153 0 : logSink() << " Selecting with TaQL: '" << msSelect << "'" << endl;
4154 0 : logSink() << LogIO::POST;
4155 :
4156 :
4157 : // Assume no selection, for starters
4158 : // gmoellen 2012/01/30 mssel_p = new MeasurementSet(sorted, ms_p);
4159 0 : mssel_p = new MeasurementSet(*ms_p);
4160 :
4161 : // Apply user-supplied selection
4162 0 : Bool nontrivsel=false;
4163 : // gmoellen 2012/01/30 nontrivsel= mssSetData(MeasurementSet(sorted, ms_p),
4164 :
4165 : // Ensure use of a fresh MSSelection object
4166 0 : if (mss_p) { delete mss_p; mss_p=NULL; }
4167 0 : mss_p=new MSSelection();
4168 0 : nontrivsel= mssSetData(*ms_p,
4169 0 : *mssel_p,"",
4170 : time,baseline,
4171 : field,spw,
4172 : uvrange,msSelect,
4173 : "",scan,"",intent, obsIDs,mss_p);
4174 :
4175 : // Keep any MR status for the MS
4176 0 : mssel_p->setMemoryResidentSubtables(ms_p->getMrsEligibility());
4177 :
4178 : // If non-trivial MSSelection invoked and nrow reduced:
4179 0 : if(nontrivsel && mssel_p->nrow()<ms_p->nrow()) {
4180 :
4181 : // Escape if no rows selected
4182 0 : if (mssel_p->nrow()==0)
4183 0 : throw(AipsError("Specified selection selects zero rows!"));
4184 :
4185 : // ...otherwise report how many rows are selected
4186 0 : logSink() << "By selection " << ms_p->nrow()
4187 0 : << " rows are reduced to " << mssel_p->nrow()
4188 0 : << LogIO::POST;
4189 : }
4190 : else {
4191 : // Selection did nothing:
4192 0 : logSink() << "Selection did not drop any rows" << LogIO::POST;
4193 : }
4194 :
4195 : // Now, re-create the associated VisSet
4196 0 : if(vs_p) delete vs_p; vs_p=0;
4197 0 : Block<int> sort(0);
4198 0 : Matrix<Int> noselection;
4199 : // gmoellen 2012/01/30 vs_p = new VisSet(*mssel_p,sort,noselection);
4200 0 : vs_p = new VisSet(*mssel_p,sort,noselection,false,0.0,false,false);
4201 0 : AlwaysAssert(vs_p, AipsError);
4202 :
4203 : // Attempt to use MSSelection for channel selection
4204 : // if user not using the old way
4205 0 : if (chanmode=="none") {
4206 0 : selectChannel(spw);
4207 : }
4208 : else {
4209 : // Reluctantly use the old-fashioned way
4210 : logSink() << LogIO::WARN
4211 : << "You have used the old-fashioned mode parameter" << endl
4212 : << "for channel selection. It still works, for now," << endl
4213 : << "but this will be eliminated in the near future." << endl
4214 : << "Please begin using the new channel selection" << endl
4215 0 : << "syntax in the spw parameter." << LogIO::POST;
4216 0 : selectChannel(chanmode,nchan,start,step,mStart,mStep);
4217 : }
4218 :
4219 0 : }
4220 0 : catch (MSSelectionError& x) {
4221 : // Re-initialize with the existing MS
4222 0 : logSink() << LogOrigin("Calibrater","selectvis",WHERE)
4223 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4224 0 : << LogIO::POST;
4225 : // jagonzal (CAS-4110): I guess it is not necessary to create these columns when the selection is empty
4226 0 : initialize(*ms_p,false,false,false);
4227 0 : throw(AipsError("Error in data selection specification: " + x.getMesg()));
4228 0 : }
4229 0 : catch (AipsError x) {
4230 : // Re-initialize with the existing MS
4231 0 : logSink() << LogOrigin("Calibrater","selectvis",WHERE)
4232 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4233 0 : << LogIO::POST;
4234 : // jagonzal (CAS-4110): I guess it is not necessary to create these columns when the selection is empty.
4235 0 : initialize(*ms_p,false,false,false);
4236 0 : throw(AipsError("Error in Calibrater::selectvis(): " + x.getMesg()));
4237 0 : }
4238 0 : };
4239 :
4240 :
4241 0 : Bool OldCalibrater::setapply (const String& type,
4242 : const Record& applypar)
4243 : {
4244 0 : logSink() << LogOrigin("Calibrater", "setapply(type, applypar)");
4245 :
4246 : // First try to create the requested VisCal object
4247 0 : VisCal *vc(NULL);
4248 :
4249 : try {
4250 :
4251 0 : if(!ok())
4252 0 : throw(AipsError("Calibrater not prepared for setapply."));
4253 :
4254 0 : String upType=type;
4255 0 : upType.upcase();
4256 :
4257 : logSink() << LogIO::NORMAL
4258 : << "Arranging to APPLY:"
4259 0 : << LogIO::POST;
4260 :
4261 : // Add a new VisCal to the apply list
4262 0 : vc = createVisCal(upType,*vs_p);
4263 :
4264 0 : vc->setApply(applypar);
4265 :
4266 : logSink() << LogIO::NORMAL << ". "
4267 0 : << vc->applyinfo()
4268 0 : << LogIO::POST;
4269 :
4270 0 : } catch (AipsError x) {
4271 : logSink() << LogIO::SEVERE << x.getMesg()
4272 : << " Check inputs and try again."
4273 0 : << LogIO::POST;
4274 0 : if (vc) delete vc;
4275 0 : throw(AipsError("Error in Calibrater::setapply."));
4276 : return false;
4277 0 : }
4278 :
4279 : // Creation apparently successful, so add to the apply list
4280 : // TBD: consolidate with above?
4281 : try {
4282 :
4283 0 : uInt napp=vc_p.nelements();
4284 0 : vc_p.resize(napp+1,false,true);
4285 0 : vc_p[napp] = vc;
4286 0 : vc=NULL;
4287 :
4288 : // Maintain sort of apply list
4289 0 : ve_p->setapply(vc_p);
4290 :
4291 0 : return true;
4292 :
4293 0 : } catch (AipsError x) {
4294 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4295 0 : << LogIO::POST;
4296 0 : if (vc) delete vc;
4297 0 : throw(AipsError("Error in Calibrater::setapply."));
4298 : return false;
4299 0 : }
4300 : return false;
4301 : }
4302 :
4303 : // Set up apply-able calibration via a Cal Library
4304 0 : Bool OldCalibrater::setcallib(Record callib) {
4305 :
4306 0 : logSink() << LogOrigin("Calibrater", "setcallib(callib)");
4307 :
4308 : // cout << "Calibrater::setcallib: callib.isFixed() = " << boolalpha << callib.isFixed() << endl;
4309 :
4310 0 : uInt ntab=callib.nfields();
4311 :
4312 : // cout << "callib.nfields() = " << ntab << endl;
4313 :
4314 : // Do some preliminary per-table verification
4315 0 : for (uInt itab=0;itab<ntab;++itab) {
4316 :
4317 0 : String tabname=callib.name(itab);
4318 :
4319 : // Insist that the table exists on disk
4320 0 : if (!Table::isReadable(tabname))
4321 0 : throw(AipsError("Caltable "+tabname+" does not exist."));
4322 :
4323 0 : }
4324 :
4325 : // Tables exist, so deploy them...
4326 :
4327 0 : for (uInt itab=0;itab<ntab;++itab) {
4328 :
4329 0 : String tabname=callib.name(itab);
4330 :
4331 : // Get the type from the table
4332 0 : String upType=calTableType(tabname);
4333 0 : upType.upcase();
4334 :
4335 : // Add table name to the record
4336 0 : Record thistabrec=callib.asrwRecord(itab);
4337 0 : thistabrec.define("tablename",tabname);
4338 :
4339 : // First try to create the requested VisCal object
4340 0 : VisCal *vc(NULL);
4341 :
4342 : try {
4343 :
4344 0 : if(!ok())
4345 0 : throw(AipsError("Calibrater not prepared for setapply."));
4346 :
4347 : logSink() << LogIO::NORMAL
4348 : << "Arranging to APPLY:"
4349 0 : << LogIO::POST;
4350 :
4351 : // Add a new VisCal to the apply list
4352 0 : vc = createVisCal(upType,*vs_p);
4353 :
4354 :
4355 : // ingest this table according to its callib
4356 0 : vc->setCallib(thistabrec,*mssel_p);
4357 :
4358 0 : } catch (AipsError x) {
4359 : logSink() << LogIO::SEVERE << x.getMesg()
4360 : << " Check inputs and try again."
4361 0 : << LogIO::POST;
4362 0 : if (vc) delete vc;
4363 0 : throw(AipsError("Error in Calibrater::setapply."));
4364 : return false;
4365 0 : }
4366 :
4367 : // Creation apparently successful, so add to the apply list
4368 : // TBD: consolidate with above?
4369 : try {
4370 :
4371 0 : uInt napp=vc_p.nelements();
4372 0 : vc_p.resize(napp+1,false,true);
4373 0 : vc_p[napp] = vc;
4374 0 : vc=NULL;
4375 :
4376 : // Maintain sort of apply list
4377 0 : ve_p->setapply(vc_p);
4378 :
4379 0 : } catch (AipsError x) {
4380 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4381 0 : << LogIO::POST;
4382 0 : if (vc) delete vc;
4383 0 : throw(AipsError("Error in Calibrater::setapply."));
4384 : return false;
4385 0 : }
4386 0 : }
4387 :
4388 : // All ok, if we get this far!
4389 0 : return true;
4390 :
4391 : }
4392 :
4393 :
4394 : // Set up apply-able calibration via a Cal Library
4395 0 : Bool OldCalibrater::setcallib2(Record callib, const casacore::MeasurementSet* ms) {
4396 :
4397 0 : logSink() << LogOrigin("OldCalibrater", "setcallib2(callib)");
4398 :
4399 : // cout << "Calibrater::setcallib2(callib) : " << boolalpha << callib << endl;
4400 :
4401 0 : uInt ntab=callib.nfields();
4402 :
4403 : // cout << "callib.nfields() = " << ntab << endl;
4404 :
4405 : // Do some preliminary per-table verification
4406 0 : for (uInt itab=0;itab<ntab;++itab) {
4407 :
4408 0 : String tabname=callib.name(itab);
4409 :
4410 : // Trap parang
4411 : // TBD...
4412 : // if (tabname=="<parang>")
4413 : // continue;
4414 :
4415 : // Insist that the table exists on disk
4416 0 : if (!Table::isReadable(tabname))
4417 0 : throw(AipsError("Caltable "+tabname+" does not exist."));
4418 :
4419 0 : }
4420 :
4421 : // Tables exist, so deploy them...
4422 :
4423 : // Local MS object for callib parsing (only)
4424 : // MeasurementSet lms(msname_p,Table::Update);
4425 : //cout << "OLD lms" << endl;
4426 :
4427 :
4428 : // Local const MS object for callib parsing (only)
4429 0 : const MeasurementSet *lmsp(0);
4430 0 : if (ms) {
4431 : // Use supplied MS (from outside), if specified...
4432 : // TBD: should we verify same base MS as ms_p/mssel_p?
4433 0 : lmsp=ms;
4434 : }
4435 : else {
4436 : // ...use internal one instead
4437 0 : lmsp=mssel_p;
4438 : }
4439 : // Reference for use below
4440 0 : const MeasurementSet &lms(*lmsp);
4441 :
4442 : // Get some global shape info:
4443 0 : Int MSnAnt = lms.antenna().nrow();
4444 0 : Int MSnSpw = lms.spectralWindow().nrow();
4445 :
4446 0 : for (uInt itab=0;itab<ntab;++itab) {
4447 :
4448 0 : String tabname=callib.name(itab);
4449 :
4450 : // Get the type from the table
4451 0 : String upType=calTableType(tabname);
4452 0 : upType.upcase();
4453 :
4454 : // Add table name to the record
4455 0 : Record thistabrec=callib.asrwRecord(itab);
4456 0 : thistabrec.define("tablename",tabname);
4457 :
4458 : // First try to create the requested VisCal object
4459 0 : VisCal *vc(NULL);
4460 :
4461 : try {
4462 :
4463 : // if(!ok())
4464 : // throw(AipsError("Calibrater not prepared for setapply."));
4465 :
4466 : logSink() << LogIO::NORMAL
4467 : << "Arranging to APPLY:"
4468 0 : << LogIO::POST;
4469 :
4470 : // Add a new VisCal to the apply list
4471 0 : vc = createVisCal(upType,msname_p,MSnAnt,MSnSpw);
4472 :
4473 : // ingest this table according to its callib
4474 0 : vc->setCallib(thistabrec,lms);
4475 :
4476 0 : } catch (AipsError x) {
4477 : logSink() << LogIO::SEVERE << x.getMesg()
4478 : << " Check inputs and try again."
4479 0 : << LogIO::POST;
4480 0 : if (vc) delete vc;
4481 0 : throw(AipsError("Error in Calibrater::callib2."));
4482 : return false;
4483 0 : }
4484 :
4485 : // Creation apparently successful, so add to the apply list
4486 : // TBD: consolidate with above?
4487 : try {
4488 :
4489 0 : uInt napp=vc_p.nelements();
4490 0 : vc_p.resize(napp+1,false,true);
4491 0 : vc_p[napp] = vc;
4492 0 : vc=NULL;
4493 :
4494 : // Maintain sort of apply list
4495 0 : ve_p->setapply(vc_p);
4496 :
4497 0 : } catch (AipsError x) {
4498 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4499 0 : << LogIO::POST;
4500 0 : if (vc) delete vc;
4501 0 : throw(AipsError("Error in Calibrater::setapply."));
4502 : return false;
4503 0 : }
4504 0 : }
4505 : // All ok, if we get this far!
4506 0 : return true;
4507 :
4508 : }
4509 :
4510 :
4511 0 : Bool OldCalibrater::setsolve (const String& type,
4512 : const Record& solvepar) {
4513 :
4514 : // Attempt to create the solvable object
4515 0 : SolvableVisCal *svc(NULL);
4516 : try {
4517 :
4518 0 : if(!ok())
4519 0 : throw(AipsError("Calibrater not prepared for setsolve."));
4520 :
4521 0 : String upType = type;
4522 0 : upType.upcase();
4523 :
4524 : // Clean out any old solve that was lying around
4525 0 : unsetsolve();
4526 :
4527 : logSink() << LogIO::NORMAL
4528 : << "Arranging to SOLVE:"
4529 0 : << LogIO::POST;
4530 :
4531 : // Create the new SolvableVisCal
4532 0 : svc = createSolvableVisCal(upType,*vs_p);
4533 0 : svc->setSolve(solvepar);
4534 :
4535 : logSink() << LogIO::NORMAL << ". "
4536 0 : << svc->solveinfo()
4537 0 : << LogIO::POST;
4538 :
4539 : // Creation apparently successful, keep it
4540 0 : svc_p=svc;
4541 0 : svc=NULL;
4542 :
4543 0 : return true;
4544 :
4545 0 : } catch (AipsError x) {
4546 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4547 0 : << LogIO::POST;
4548 0 : unsetsolve();
4549 0 : if (svc) delete svc;
4550 0 : throw(AipsError("Error in Calibrater::setsolve."));
4551 : return false;
4552 0 : }
4553 : return false;
4554 : }
4555 :
4556 0 : Bool OldCalibrater::correct(String mode)
4557 : {
4558 0 : logSink() << LogOrigin("Calibrater","correct") << LogIO::NORMAL;
4559 :
4560 0 : Bool retval = true;
4561 :
4562 : try {
4563 :
4564 : // make mode all-caps
4565 0 : String upmode=mode;
4566 0 : upmode.upcase();
4567 :
4568 : // If trialmode=T, only the flags will be set
4569 : // (and only written if not TRIAL)
4570 0 : Bool trialmode=(upmode.contains("TRIAL") ||
4571 0 : upmode.contains("FLAGONLY"));
4572 :
4573 : // Set up VisSet and its VisibilityIterator.
4574 :
4575 0 : VisibilityIterator::DataColumn whichOutCol = configureForCorrection ();
4576 :
4577 0 : VisIter& vi(vs_p->iter());
4578 0 : VisBufferAutoPtr vb (vi);
4579 0 : vi.origin();
4580 :
4581 : // Pass each timestamp (VisBuffer) to VisEquation for correction
4582 :
4583 0 : Vector<Bool> uncalspw(vi.numberSpw()); // Used to accumulate error messages
4584 0 : uncalspw.set(false); // instead of bombing the user
4585 0 : uncalspw.set(False); // instead of bombing the user
4586 : // in a loop.
4587 :
4588 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
4589 :
4590 0 : for (vi.origin(); vi.more(); vi++) {
4591 :
4592 0 : uInt spw = vb->spectralWindow();
4593 :
4594 : // Re-initialize weights from sigma column
4595 0 : vb->resetWeightMat();
4596 :
4597 : // If we can calibrate this vb, do it...
4598 0 : if (ve_p->spwOK(spw)){
4599 :
4600 : // throws exception if nothing to apply
4601 0 : ve_p->correct(*vb,trialmode);
4602 :
4603 : }
4604 : // ...else don't, prepare warning, and possibly set flags
4605 : else{
4606 :
4607 : // set uncalspw for warning message
4608 0 : uncalspw[spw] = true;
4609 : // set the flags, if we are being strict
4610 0 : if (upmode.contains("STRICT"))
4611 : // set the flags
4612 : // (don't touch the data/weights, which are initialized)
4613 0 : vb->flag().set(true);
4614 : }
4615 :
4616 : // Only if not a trial run, trigger write to disk
4617 0 : if (!upmode.contains("TRIAL")) {
4618 :
4619 0 : if (upmode.contains("CAL")) {
4620 0 : vi.setVis (vb->visCube(), whichOutCol);
4621 0 : vi.setWeightMat(vb->weightMat());
4622 : }
4623 :
4624 0 : if (upmode.contains("FLAG"))
4625 0 : vi.setFlag (vb->flag());
4626 :
4627 : }
4628 :
4629 : }
4630 : }
4631 :
4632 0 : vs_p->flush (); // Flush to disk
4633 :
4634 : // Now that we're out of the loop, summarize any errors.
4635 :
4636 0 : retval = summarize_uncalspws(uncalspw, "correct",
4637 0 : upmode.contains("STRICT"));
4638 :
4639 0 : actRec_=Record();
4640 0 : actRec_.define("origin","Calibrater::correct");
4641 0 : actRec_.defineRecord("VisEquation",ve_p->actionRec());
4642 :
4643 0 : }
4644 0 : catch (AipsError x) {
4645 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4646 0 : << LogIO::POST;
4647 :
4648 0 : logSink() << "Resetting all calibration application settings." << LogIO::POST;
4649 0 : unsetapply();
4650 :
4651 0 : throw(AipsError("Error in Calibrater::correct."));
4652 : retval = false; // Not that it ever gets here...
4653 0 : }
4654 0 : return retval;
4655 : }
4656 :
4657 0 : Bool OldCalibrater::corrupt() {
4658 :
4659 0 : logSink() << LogOrigin("Calibrater","corrupt") << LogIO::NORMAL;
4660 0 : Bool retval = true;
4661 :
4662 : try {
4663 :
4664 0 : if (!ok())
4665 0 : throw(AipsError("Calibrater not prepared for corrupt!"));
4666 :
4667 : // Nominally, we write out to the MODEL_DATA, unless absent
4668 0 : VisibilityIterator::DataColumn whichOutCol(VisibilityIterator::Model);
4669 :
4670 0 : if (!ms_p->tableDesc().isColumn("MODEL_DATA"))
4671 0 : throw(AipsError("MODEL_DATA column unexpectedly absent. Cannot corrupt."));
4672 :
4673 : // Ensure apply list non-zero and properly sorted
4674 0 : ve_p->setapply(vc_p);
4675 :
4676 : // Report the types that will be applied
4677 0 : applystate();
4678 :
4679 : // Arrange for iteration over data
4680 0 : Block<Int> columns;
4681 : // include scan iteration
4682 0 : columns.resize(5);
4683 0 : columns[0]=MS::ARRAY_ID;
4684 0 : columns[1]=MS::SCAN_NUMBER;
4685 0 : columns[2]=MS::FIELD_ID;
4686 0 : columns[3]=MS::DATA_DESC_ID;
4687 0 : columns[4]=MS::TIME;
4688 0 : vs_p->resetVisIter(columns,0.0);
4689 0 : VisIter& vi(vs_p->iter());
4690 0 : VisBuffer vb(vi);
4691 :
4692 : // Pass each timestamp (VisBuffer) to VisEquation for corruption.
4693 0 : Vector<Bool> uncalspw(vi.numberSpw()); // Used to accumulate error messages
4694 0 : uncalspw.set(false); // instead of bombing the user
4695 : // in a loop.
4696 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
4697 0 : Int spw = vi.spectralWindow();
4698 :
4699 : // Only proceed if spw can be calibrated
4700 0 : if (ve_p->spwOK(spw)) {
4701 :
4702 0 : for (vi.origin(); vi.more(); vi++) {
4703 :
4704 : // Corrupt the MODEL_DATA
4705 : // (note we are not treating weights and flags)
4706 0 : ve_p->corrupt(vb); // throws exception if nothing to apply
4707 0 : vi.setVis(vb.modelVisCube(),whichOutCol);
4708 :
4709 : }
4710 : }
4711 : else
4712 0 : uncalspw[spw] = true;
4713 : }
4714 : // Flush to disk
4715 0 : vs_p->flush();
4716 :
4717 : // Now that we're out of the loop, summarize any errors.
4718 0 : retval = summarize_uncalspws(uncalspw, "corrupt");
4719 0 : }
4720 0 : catch (AipsError x) {
4721 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4722 0 : << LogIO::POST;
4723 :
4724 0 : logSink() << "Resetting all calibration application settings." << LogIO::POST;
4725 0 : unsetapply();
4726 :
4727 0 : throw(AipsError("Error in Calibrater::corrupt."));
4728 : retval = false; // Not that it ever gets here...
4729 0 : }
4730 0 : return retval;
4731 : }
4732 :
4733 0 : Bool OldCalibrater::initWeightsWithTsys(String wtmode, Bool dowtsp,
4734 : String tsystable, String gainfield, String interp, Vector<Int> spwmap) {
4735 :
4736 0 : logSink() << LogOrigin("Calibrater", "initWeightsWithTsys")
4737 0 : << LogIO::NORMAL;
4738 0 : Bool retval = true;
4739 :
4740 : try {
4741 :
4742 0 : if (!ok())
4743 0 : throw(AipsError("Calibrater not prepared for initWeights!"));
4744 :
4745 0 : String uptype = calTableType(tsystable);
4746 0 : if (!uptype.contains("TSYS")) {
4747 0 : throw(AipsError(
4748 0 : "Invalid calibration table type for Tsys weighting."));
4749 : }
4750 : // Set record format for calibration table application information
4751 0 : RecordDesc applyparDesc;
4752 0 : applyparDesc.addField("t", TpDouble);
4753 0 : applyparDesc.addField("table", TpString);
4754 0 : applyparDesc.addField("interp", TpString);
4755 0 : applyparDesc.addField("spw", TpArrayInt);
4756 0 : applyparDesc.addField("fieldstr", TpString);
4757 0 : applyparDesc.addField("calwt", TpBool);
4758 0 : applyparDesc.addField("spwmap", TpArrayInt);
4759 0 : applyparDesc.addField("opacity", TpArrayDouble);
4760 :
4761 : // Create record with the requisite field values
4762 0 : Record applypar(applyparDesc);
4763 0 : applypar.define("t", 0.0);
4764 0 : applypar.define("table", tsystable);
4765 0 : applypar.define("interp", interp);
4766 0 : applypar.define("spw", getSpwIdx(""));
4767 0 : applypar.define("fieldstr", gainfield);
4768 0 : applypar.define("calwt", true);
4769 0 : applypar.define("spwmap", spwmap);
4770 0 : applypar.define("opacity", Vector<Double>(1, 0.0));
4771 :
4772 0 : if (vc_p.nelements() > 0) {
4773 0 : logSink() << LogIO::WARN << "Resetting all calibration application settings." << LogIO::POST;
4774 0 : unsetapply();
4775 : }
4776 0 : logSink() << LogIO::NORMAL << "Weight initialization does not support selection. Resetting MS selection." << LogIO::POST;
4777 0 : selectvis();
4778 0 : StandardTsys vc = StandardTsys(*vs_p);
4779 0 : vc.setApply(applypar);
4780 :
4781 0 : logSink() << LogIO::NORMAL << ". " << vc.applyinfo() << LogIO::POST;
4782 0 : PtrBlock<VisCal*> vcb(1, &vc);
4783 : // Maintain sort of apply list
4784 0 : ve_p->setapply(vcb);
4785 :
4786 : // Detect WEIGHT_SPECTRUM and SIGMA_SPECTRUM
4787 0 : TableDesc mstd = ms_p->actualTableDesc();
4788 0 : String colWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
4789 0 : Bool wtspexists = mstd.isColumn(colWtSp);
4790 0 : String colSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
4791 0 : Bool sigspexists = mstd.isColumn(colSigSp);
4792 0 : Bool addsigsp = (dowtsp && !sigspexists);
4793 :
4794 : // Some log info
4795 0 : bool use_exposure = false;
4796 0 : if (wtmode == "tsys") {
4797 : logSink()
4798 : << "Initializing SIGMA and WEIGHT according to channel bandwidth and Tsys. NOTE this is an expert mode."
4799 0 : << LogIO::WARN << LogIO::POST;
4800 0 : } else if (wtmode == "tinttsys") {
4801 : logSink()
4802 : << "Initializing SIGMA and WEIGHT according to channel bandwidth, integration time, and Tsys. NOTE this is an expert mode."
4803 0 : << LogIO::WARN << LogIO::POST;
4804 0 : use_exposure = true;
4805 : } else {
4806 0 : throw(AipsError("Unrecognized wtmode specified: " + wtmode));
4807 : }
4808 :
4809 : // Force dowtsp if the column already exists
4810 0 : if (wtspexists && !dowtsp) {
4811 : logSink() << "Found WEIGHT_SPECTRUM; will force its initialization."
4812 0 : << LogIO::POST;
4813 0 : dowtsp = true;
4814 : }
4815 :
4816 : // Report that we are initializing the WEIGHT_SPECTRUM, and prepare to do so.
4817 0 : if (dowtsp) {
4818 :
4819 : // Ensure WEIGHT_SPECTRUM really exists at all
4820 : // (often it exists but is empty)
4821 0 : if (!wtspexists) {
4822 0 : logSink() << "Creating WEIGHT_SPECTRUM." << LogIO::POST;
4823 :
4824 : // Nominal defaulttileshape
4825 0 : IPosition dts(3, 4, 32, 1024);
4826 :
4827 : // Discern DATA's default tile shape and use it
4828 0 : const Record dminfo = ms_p->dataManagerInfo();
4829 0 : for (uInt i = 0; i < dminfo.nfields(); ++i) {
4830 0 : Record col = dminfo.asRecord(i);
4831 : //if (upcase(col.asString("NAME"))=="TILEDDATA") {
4832 0 : if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
4833 0 : dts = IPosition(
4834 0 : col.asRecord("SPEC").asArrayInt(
4835 0 : "DEFAULTTILESHAPE"));
4836 : //cout << "Found DATA's default tile: " << dts << endl;
4837 0 : break;
4838 : }
4839 0 : }
4840 :
4841 : // Add the column
4842 0 : String colWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
4843 0 : TableDesc tdWtSp;
4844 0 : tdWtSp.addColumn(
4845 0 : ArrayColumnDesc<Float>(colWtSp, "weight spectrum", 2));
4846 0 : TiledShapeStMan wtSpStMan("TiledWgtSpectrum", dts);
4847 0 : ms_p->addColumn(tdWtSp, wtSpStMan);
4848 0 : } else
4849 0 : logSink() << "Found WEIGHT_SPECTRUM." << LogIO::POST;
4850 : // Ensure WEIGHT_SPECTRUM really exists at all
4851 : // (often it exists but is empty)
4852 0 : if (!sigspexists) {
4853 0 : logSink() << "Creating SIGMA_SPECTRUM." << LogIO::POST;
4854 :
4855 : // Nominal defaulttileshape
4856 0 : IPosition dts(3, 4, 32, 1024);
4857 :
4858 : // Discern DATA's default tile shape and use it
4859 0 : const Record dminfo = ms_p->dataManagerInfo();
4860 0 : for (uInt i = 0; i < dminfo.nfields(); ++i) {
4861 0 : Record col = dminfo.asRecord(i);
4862 : //if (upcase(col.asString("NAME"))=="TILEDDATA") {
4863 0 : if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
4864 0 : dts = IPosition(
4865 0 : col.asRecord("SPEC").asArrayInt(
4866 0 : "DEFAULTTILESHAPE"));
4867 : //cout << "Found DATA's default tile: " << dts << endl;
4868 0 : break;
4869 : }
4870 0 : }
4871 :
4872 : // Add the column
4873 0 : String colSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
4874 0 : TableDesc tdSigSp;
4875 0 : tdSigSp.addColumn(
4876 0 : ArrayColumnDesc<Float>(colSigSp, "sigma spectrum", 2));
4877 0 : TiledShapeStMan sigSpStMan("TiledSigtSpectrum", dts);
4878 0 : ms_p->addColumn(tdSigSp, sigSpStMan);
4879 : {
4880 0 : TableDesc loctd = ms_p->actualTableDesc();
4881 0 : String loccolSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
4882 0 : AlwaysAssert(loctd.isColumn(loccolSigSp),AipsError);
4883 0 : }
4884 0 : }
4885 : }
4886 : else {
4887 0 : if (sigspexists) {
4888 0 : logSink() << "Removing SIGMA_SPECTRUM for non-channelized weight." << LogIO::POST;
4889 : if (true || ms_p->canRemoveColumn(colSigSp)) {
4890 0 : ms_p->removeColumn(colSigSp);
4891 : }
4892 : else
4893 : logSink() << LogIO::WARN << "Failed to remove SIGMA_SPECTRUM column. Values in SIGMA and SIGMA_SPECTRUM columns may be inconsistent after the operation." << LogIO::POST;
4894 : }
4895 : }
4896 :
4897 : // Arrange for iteration over data
4898 : // TBD: Be sure this sort is optimal for creating WS?
4899 0 : Block<Int> columns;
4900 : // include scan iteration
4901 0 : columns.resize(5);
4902 0 : columns[0] = MS::ARRAY_ID;
4903 0 : columns[1] = MS::SCAN_NUMBER;
4904 0 : columns[2] = MS::FIELD_ID;
4905 0 : columns[3] = MS::DATA_DESC_ID;
4906 0 : columns[4] = MS::TIME;
4907 :
4908 0 : vi::SortColumns sc(columns);
4909 0 : vi::VisibilityIterator2 vi2(*ms_p, sc, true);
4910 0 : vi::VisBuffer2 *vb = vi2.getVisBuffer();
4911 :
4912 0 : MSColumns mscol(*ms_p);
4913 0 : const MSSpWindowColumns& msspw(mscol.spectralWindow());
4914 0 : uInt nSpw = msspw.nrow();
4915 0 : Vector<Double> effChBw(nSpw, 0.0);
4916 0 : for (uInt ispw = 0; ispw < nSpw; ++ispw) {
4917 0 : effChBw[ispw] = mean(msspw.effectiveBW()(ispw));
4918 : }
4919 :
4920 0 : Int ivb(0);
4921 0 : for (vi2.originChunks(); vi2.moreChunks(); vi2.nextChunk()) {
4922 :
4923 0 : for (vi2.origin(); vi2.more(); vi2.next(), ++ivb) {
4924 :
4925 0 : Int spw = vb->spectralWindows()(0);
4926 :
4927 0 : auto nrow = vb->nRows();
4928 0 : Int nchan = vb->nChannels();
4929 0 : Int ncor = vb->nCorrelations();
4930 :
4931 : // Prepare for WEIGHT_SPECTRUM and SIGMA_SPECTRUM, if nec.
4932 0 : Cube<Float> newwtsp(0, 0, 0), newsigsp(0, 0, 0);
4933 0 : if (dowtsp) {
4934 0 : newwtsp.resize(ncor, nchan, nrow);
4935 0 : newwtsp.set(1.0);
4936 0 : newsigsp.resize(ncor, nchan, nrow);
4937 0 : newsigsp.set(1.0);
4938 : }
4939 :
4940 0 : if (ve_p->spwOK(spw)) {
4941 :
4942 : // Re-initialize weight info from sigma info
4943 : // This is smart wrt spectral weights, etc.
4944 : // (this makes W and WS, if present, "dirty" in the vb)
4945 : // TBD: only do this if !trial (else: avoid the I/O)
4946 : // vb->resetWeightsUsingSigma();
4947 : // Handle non-trivial modes
4948 : // Init WEIGHT, SIGMA from bandwidth & time
4949 0 : Matrix<Float> newwt(ncor, nrow), newsig(ncor, nrow);
4950 0 : newwt.set(1.0);
4951 0 : newsig.set(1.0);
4952 :
4953 : // Detect ACs
4954 0 : const Vector<Int> a1(vb->antenna1());
4955 0 : const Vector<Int> a2(vb->antenna2());
4956 0 : Vector<Bool> ac(a1 == a2);
4957 :
4958 : // XCs need an extra factor of 2
4959 0 : Vector<Float> xcfactor(nrow, 2.0);
4960 0 : xcfactor(ac) = 1.0; // (but not ACs)
4961 :
4962 : // The row-wise integration time
4963 0 : Vector<Float> expo(nrow);
4964 0 : convertArray(expo, vb->exposure());
4965 :
4966 : // Set weights to channel bandwidth first.
4967 0 : newwt.set(Float(effChBw(spw)));
4968 :
4969 : // For each correlation, apply exposure and xcfactor
4970 0 : for (Int icor = 0; icor < ncor; ++icor) {
4971 :
4972 0 : Vector<Float> wt(newwt.row(icor));
4973 0 : if (use_exposure) {
4974 0 : wt *= expo;
4975 : }
4976 0 : wt *= xcfactor;
4977 0 : if (dowtsp) {
4978 0 : for (Int ich = 0; ich < nchan; ++ich) {
4979 : Vector<Float> wtspi(
4980 0 : newwtsp(Slice(icor, 1, 1),
4981 0 : Slice(ich, 1, 1), Slice()).nonDegenerate(
4982 0 : IPosition(1, 2)));
4983 0 : wtspi = wt;
4984 0 : }
4985 : }
4986 0 : }
4987 : // Handle SIGMA_SPECTRUM
4988 0 : if (dowtsp) {
4989 0 : newsigsp = 1.0f / sqrt(newwtsp);
4990 : }
4991 : // sig from wt is inverse sqrt
4992 0 : newsig = 1.0f / sqrt(newwt);
4993 :
4994 : // Arrange write-back of both SIGMA and WEIGHT
4995 0 : vb->setSigma(newsig);
4996 0 : vb->setWeight(newwt);
4997 0 : if (dowtsp) {
4998 0 : vb->initWeightSpectrum(newwtsp);
4999 0 : vb->initSigmaSpectrum(newsigsp);
5000 : }
5001 : // Force writeback to disk (need to initialize weight/sigma before applying cal table)
5002 0 : vb->writeChangesBack();
5003 :
5004 : // Arrange for _in-place_ apply on CORRECTED_DATA (init from DATA)
5005 : // (this makes CD "dirty" in the vb)
5006 : // TBD: only do this if !trial (else: avoid the I/O)
5007 0 : vb->setVisCubeCorrected(vb->visCube());
5008 :
5009 : // Make flagcube dirty in the vb
5010 : // NB: we must _always_ do this I/O (even trial mode)
5011 0 : vb->setFlagCube(vb->flagCube());
5012 :
5013 : // Make all vb "not dirty"; we'll carefully arrange the writeback below
5014 0 : vb->dirtyComponentsClear();
5015 :
5016 : // throws exception if nothing to apply
5017 0 : ve_p->correct2(*vb, false, dowtsp);
5018 :
5019 0 : if (dowtsp) {
5020 0 : vb->setWeightSpectrum(vb->weightSpectrum());
5021 : // If WS was calibrated, set W to its channel-axis median
5022 0 : vb->setWeight( partialMedians(vb->weightSpectrum(), IPosition(1, 1)) );
5023 0 : newsigsp = 1.0f / sqrt(vb->weightSpectrum());
5024 0 : vb->initSigmaSpectrum(newsigsp);
5025 0 : vb->setSigma( partialMedians(newsigsp, IPosition(1, 1)) );
5026 : } else {
5027 0 : vb->setWeight(vb->weight());
5028 0 : newsig = 1.0f / sqrt(vb->weight());
5029 0 : vb->setSigma(newsig);
5030 : }
5031 : // Force writeback to disk
5032 0 : vb->writeChangesBack();
5033 :
5034 0 : } else {//Not calibrating the spw
5035 0 : if (dowtsp && !wtspexists) {
5036 : // newly created WS Need to initialize
5037 0 : vb->initWeightSpectrum(newwtsp);
5038 : }
5039 0 : if (addsigsp) {
5040 : // newly created SS Need to initialize
5041 0 : vb->initSigmaSpectrum(newsigsp);
5042 0 : vb->writeChangesBack();
5043 : }
5044 : }
5045 0 : }
5046 : }
5047 : // clear-up Tsys caltable from list of apply
5048 0 : unsetapply();
5049 :
5050 0 : } catch (AipsError x) {
5051 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
5052 0 : << LogIO::POST;
5053 :
5054 : logSink() << "Resetting all calibration application settings."
5055 0 : << LogIO::POST;
5056 0 : unsetapply();
5057 :
5058 0 : throw(AipsError("Error in Calibrater::initWeights."));
5059 : retval = false; // Not that it ever gets here...
5060 0 : }
5061 0 : return retval;
5062 : }
5063 :
5064 0 : Bool OldCalibrater::solve() {
5065 :
5066 0 : logSink() << LogOrigin("Calibrater","solve") << LogIO::NORMAL3;
5067 :
5068 : try {
5069 :
5070 0 : if (!ok())
5071 0 : throw(AipsError("Calibrater not prepared for solve."));
5072 :
5073 : // Handle nothing-to-solve-for case
5074 0 : if (!svc_p)
5075 0 : throw(AipsError("Please run setsolve before attempting to solve."));
5076 :
5077 : // Handle specified caltable
5078 0 : if (true && svc_p) {
5079 :
5080 : /*
5081 : cout << "name: " << svc_p->calTableName() << endl;
5082 : cout << boolalpha;
5083 : cout << "append? " << svc_p->append() << endl;
5084 : cout << "opened? " << Table::isOpened(svc_p->calTableName()) << endl;
5085 : cout << "readable? " << Table::isReadable(svc_p->calTableName()) << endl;
5086 : // cout << "writable? " << Table::isWritable(svc_p->calTableName()) << endl;
5087 : cout << "canDelete? " << Table::canDeleteTable(svc_p->calTableName(),true) << endl;
5088 : */
5089 :
5090 : // If table exists (readable) and not deletable
5091 : // we have to abort (append=T requires deletable)
5092 0 : if ( Table::isReadable(svc_p->calTableName()) &&
5093 0 : !TableUtil::canDeleteTable(svc_p->calTableName()) ) {
5094 0 : throw(AipsError("Specified caltable ("+svc_p->calTableName()+") exists and\n cannot be replaced (or appended to) because it appears to be open somewhere (Quit plotcal?)."));
5095 : }
5096 : }
5097 :
5098 : // Arrange VisEquation for solve
5099 0 : ve_p->setsolve(*svc_p);
5100 :
5101 : // Ensure apply list properly sorted w.r.t. solvable term
5102 0 : ve_p->setapply(vc_p);
5103 :
5104 : // Report what is being applied and solved-for
5105 0 : applystate();
5106 0 : solvestate();
5107 :
5108 :
5109 : // Report correct/corrupt apply order
5110 : // ve_p->state();
5111 :
5112 : // Set the channel mask
5113 0 : svc_p->setChanMask(chanmask_);
5114 :
5115 : // Generally use standard solver
5116 0 : if (svc_p->useGenericGatherForSolve())
5117 0 : genericGatherAndSolve(); // using VisBuffGroupAcc
5118 : else {
5119 : //cout << "Fully self-directed data gather and solve" << endl;
5120 : // Fully self-directed data gather and solve
5121 0 : svc_p->selfGatherAndSolve(*vs_p,*ve_p);
5122 : }
5123 :
5124 0 : svc_p->clearChanMask();
5125 :
5126 0 : } catch (AipsError x) {
5127 0 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
5128 :
5129 0 : logSink() << "Reseting entire solve/apply state." << LogIO::POST;
5130 0 : reset();
5131 :
5132 0 : throw(AipsError("Error in Calibrater::solve."));
5133 : return false;
5134 0 : }
5135 :
5136 0 : return true;
5137 :
5138 : }
5139 :
5140 0 : Vector<Double> OldCalibrater::modelfit(const Int& niter,
5141 : const String& stype,
5142 : const Vector<Double>& par,
5143 : const Vector<Bool>& vary,
5144 : const String& file) {
5145 :
5146 : /*
5147 : cout << "Calibrater::modelfit" << endl;
5148 : cout << " niter = " << niter << endl;
5149 : cout << " stype = " << stype << endl;
5150 : cout << " par = " << par << endl;
5151 : cout << " vary = " << vary << endl;
5152 : cout << " file = " << file << endl;
5153 : */
5154 : // logSink() << LogOrigin("Calibrater","modelfit") << LogIO::NORMAL;
5155 :
5156 : try {
5157 0 : if(!ok()) throw(AipsError("Calibrater not ok()"));
5158 :
5159 : // Construct UVMod with the VisSet
5160 0 : UVMod uvmod(*vs_p);
5161 :
5162 0 : if (stype=="P")
5163 0 : uvmod.setModel(ComponentType::POINT, par, vary);
5164 0 : else if (stype=="G")
5165 0 : uvmod.setModel(ComponentType::GAUSSIAN, par, vary);
5166 0 : else if (stype=="D")
5167 0 : uvmod.setModel(ComponentType::DISK, par, vary);
5168 : else
5169 0 : throw(AipsError("Unrecognized component type in Calibrater::modelfit."));
5170 :
5171 : // Run the fit
5172 0 : uvmod.modelfit(niter,file);
5173 :
5174 : // Return the parameter vector
5175 0 : return uvmod.par();
5176 :
5177 0 : } catch (AipsError x) {
5178 0 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
5179 0 : throw(AipsError("Error in Calibrater::modelfit."));
5180 :
5181 : return Vector<Double>();
5182 0 : }
5183 :
5184 : }
5185 :
5186 0 : void OldCalibrater::fluxscale(const String& infile,
5187 : const String& outfile,
5188 : const Vector<Int>& refField,
5189 : const Vector<Int>& refSpwMap,
5190 : const Vector<Int>& tranField,
5191 : const Bool& append,
5192 : const Float& inGainThres,
5193 : const String& antSel,
5194 : const String& timerangeSel,
5195 : const String& scanSel,
5196 : SolvableVisCal::fluxScaleStruct& oFluxScaleFactor,
5197 : const String& oListFile,
5198 : const Bool& incremental,
5199 : const Int& fitorder,
5200 : const Bool& display) {
5201 :
5202 : // throw(AipsError("Method 'fluxscale' is temporarily disabled."));
5203 :
5204 : // TBD: write inputs to MSHistory
5205 0 : logSink() << LogOrigin("Calibrater","fluxscale") << LogIO::NORMAL3;
5206 :
5207 0 : SolvableVisCal *fsvj_(NULL);
5208 : try {
5209 : // If infile is Calibration table
5210 0 : if (Table::isReadable(infile) &&
5211 0 : TableUtil::tableInfo(infile).type()=="Calibration") {
5212 :
5213 : // get calibration type
5214 0 : String caltype;
5215 0 : caltype = TableUtil::tableInfo(infile).subType();
5216 : logSink() << "Table " << infile
5217 : << " is of type: "<< caltype
5218 0 : << LogIO::POST;
5219 0 : String message="Table "+infile+" is of type: "+caltype;
5220 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
5221 :
5222 : // form selection
5223 0 : String select="";
5224 : // Selection is empty for case of no tran specification
5225 0 : if (tranField.nelements()>0) {
5226 :
5227 : // All selected fields
5228 0 : Vector<Int> allflds = concatenateArray(refField,tranField);
5229 :
5230 : // Assemble TaQL
5231 0 : ostringstream selectstr;
5232 0 : selectstr << "FIELD_ID IN [";
5233 0 : for (Int iFld=0; iFld<allflds.shape(); iFld++) {
5234 0 : if (iFld>0) selectstr << ", ";
5235 0 : selectstr << allflds(iFld);
5236 : }
5237 0 : selectstr << "]";
5238 0 : select=selectstr.str();
5239 0 : }
5240 :
5241 : // Construct proper SVC object
5242 0 : if (caltype == "G Jones") {
5243 0 : fsvj_ = createSolvableVisCal("G",*vs_p);
5244 0 : } else if (caltype == "T Jones") {
5245 0 : fsvj_ = createSolvableVisCal("T",*vs_p);
5246 : } else {
5247 : // Can't process other than G and T (add B?)
5248 0 : ostringstream typeErr;
5249 0 : typeErr << "Type " << caltype
5250 0 : << " not supported in fluxscale.";
5251 :
5252 0 : throw(AipsError(typeErr.str()));
5253 0 : }
5254 :
5255 : // fill table with selection
5256 0 : RecordDesc applyparDesc;
5257 0 : applyparDesc.addField ("table", TpString);
5258 0 : applyparDesc.addField ("select", TpString);
5259 0 : Record applypar(applyparDesc);
5260 0 : applypar.define ("table", infile);
5261 0 : applypar.define ("select", select);
5262 0 : fsvj_->setApply(applypar);
5263 :
5264 : //Bool incremental=false;
5265 : // Make fluxscale calculation
5266 0 : Vector<String> fldnames(MSFieldColumns(ms_p->field()).name().getColumn());
5267 : //fsvj_->fluxscale(refField,tranField,refSpwMap,fldnames,oFluxScaleFactor,
5268 0 : fsvj_->fluxscale(outfile,refField,tranField,refSpwMap,fldnames,inGainThres,antSel,
5269 : timerangeSel,scanSel,oFluxScaleFactor, oListFile,incremental,fitorder,display);
5270 : // oListFile);
5271 :
5272 : // If no outfile specified, use infile (overwrite!)
5273 0 : String out(outfile);
5274 0 : if (out.length()==0)
5275 0 : out = infile;
5276 :
5277 : // Store result
5278 0 : if (append) {
5279 0 : logSink() << "Appending result to " << out << LogIO::POST;
5280 0 : String message="Appending result to "+out;
5281 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
5282 0 : } else {
5283 0 : logSink() << "Storing result in " << out << LogIO::POST;
5284 0 : String message="Storing result in "+out;
5285 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
5286 0 : }
5287 0 : fsvj_->storeNCT(out,append);
5288 :
5289 : // Clean up
5290 0 : delete fsvj_;
5291 :
5292 0 : } else {
5293 : // Table not found/unreadable, or not Calibration
5294 0 : ostringstream tabErr;
5295 0 : tabErr << "File " << infile
5296 0 : << " does not exist or is not a Calibration Table.";
5297 :
5298 0 : throw(AipsError(tabErr.str()));
5299 :
5300 0 : }
5301 0 : } catch (AipsError x) {
5302 :
5303 : logSink() << LogIO::SEVERE
5304 : << "Caught Exception: "
5305 : << x.getMesg()
5306 0 : << LogIO::POST;
5307 :
5308 : // Clean up
5309 0 : if (fsvj_) delete fsvj_;
5310 :
5311 : // Write to MS History table
5312 : // String message="Caught Exception: "+x.getMesg();
5313 : // MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
5314 :
5315 0 : throw(AipsError("Error in Calibrater::fluxscale."));
5316 :
5317 : return;
5318 :
5319 0 : }
5320 0 : return;
5321 :
5322 :
5323 : }
5324 :
5325 0 : void OldCalibrater::accumulate(const String& intab,
5326 : const String& incrtab,
5327 : const String& outtab,
5328 : const String& fields,
5329 : const String& calFields,
5330 : const String& interp,
5331 : const Double& t,
5332 : const Vector<Int>& spwmap) {
5333 :
5334 : // logSink() << LogOrigin("Calibrater","accumulate") << LogIO::NORMAL;
5335 :
5336 0 : logSink() << "Beginning accumulate." << LogIO::POST;
5337 :
5338 : // SVJ objects:
5339 0 : SolvableVisCal *incal_(NULL), *incrcal_(NULL);
5340 :
5341 : try {
5342 :
5343 : /*
5344 : cout << "intab = " << intab << endl;
5345 : cout << "incrtab = " << incrtab << endl;
5346 : cout << "outtab = " << outtab << endl;
5347 : cout << "fields = " << fields << endl;
5348 : cout << "calFields = " << calFields << endl;
5349 : cout << "interp = " << interp << endl;
5350 : cout << "t = " << t << endl;
5351 : */
5352 :
5353 : // Incremental table's type sets the type we are dealing with
5354 0 : String caltype=calTableType(incrtab);
5355 :
5356 : // If no input cumulative timescale specified, then
5357 : // a valid input cumulative table must be specified
5358 0 : if (t < 0.0) {
5359 :
5360 0 : String intype=calTableType(intab);
5361 :
5362 0 : if (intype!=caltype) {
5363 :
5364 0 : ostringstream typeErr;
5365 0 : typeErr << "Table " << intab
5366 0 : << " is not the same type as "
5367 0 : << incrtab << " (" << caltype << ")";
5368 0 : throw(AipsError(typeErr.str()));
5369 0 : }
5370 0 : }
5371 :
5372 : // At this point all is ok; we will:
5373 : // o fill from intab and accumulate to it (t<0), OR
5374 : // o create a new cumulative table from scratch (t>0)
5375 :
5376 : // If creating a new cumulative table, it must span the whole dataset,
5377 : // so reset data selection to whole MS, and setup iterator
5378 0 : if (t>0.0) {
5379 0 : selectvis();
5380 0 : Block<Int> columns;
5381 0 : columns.resize(4);
5382 0 : columns[0]=MS::ARRAY_ID;
5383 0 : columns[1]=MS::TIME;
5384 0 : columns[2]=MS::FIELD_ID;
5385 0 : columns[3]=MS::DATA_DESC_ID;
5386 0 : vs_p->resetVisIter(columns,t);
5387 0 : }
5388 :
5389 : // logSink() << "Table " << infile
5390 : // << " is of type: "<< caltype
5391 : // << LogIO::POST;
5392 :
5393 0 : incal_ = createSolvableVisCal(caltype,*vs_p);
5394 0 : incrcal_ = createSolvableVisCal(caltype,*vs_p);
5395 :
5396 : // TBD: move to svj.setAccumulate?
5397 0 : if ( !(incal_->accumulatable()) ) {
5398 0 : ostringstream typeErr;
5399 0 : typeErr << "Type " << caltype
5400 0 : << " not yet supported in accumulate.";
5401 0 : throw(AipsError(typeErr.str()));
5402 0 : }
5403 :
5404 : // At this point, accumulation apparently supported,
5405 : // so continue...
5406 :
5407 : // initialize the cumulative solutions
5408 0 : incal_->setAccumulate(*vs_p,intab,"",t,-1);
5409 :
5410 :
5411 : /*
5412 : // form selection on incr table
5413 : String incrSel="";
5414 : if (calFields.shape()>0) {
5415 :
5416 : // Assemble TaQL
5417 : ostringstream selectstr;
5418 : selectstr << "FIELD_ID IN [";
5419 : for (Int iFld=0; iFld<calFields.shape(); iFld++) {
5420 : if (iFld>0) selectstr << ", ";
5421 : selectstr << calFields(iFld);
5422 : }
5423 : selectstr << "]";
5424 : incrSel=selectstr.str();
5425 : }
5426 : */
5427 :
5428 : // fill incr table with selection
5429 : logSink() << "Preparing to accumulate calibration from table: "
5430 : << incrtab
5431 0 : << LogIO::POST;
5432 :
5433 : // Set record format for calibration table application information
5434 0 : RecordDesc applyparDesc;
5435 0 : applyparDesc.addField ("t", TpDouble);
5436 0 : applyparDesc.addField ("table", TpString);
5437 : // applyparDesc.addField ("select", TpString);
5438 0 : applyparDesc.addField ("fieldstr", TpString);
5439 0 : applyparDesc.addField ("interp", TpString);
5440 0 : applyparDesc.addField ("spwmap",TpArrayInt);
5441 :
5442 : // Create record with the requisite field values
5443 0 : Record applypar(applyparDesc);
5444 0 : applypar.define ("t", t);
5445 0 : applypar.define ("table", incrtab);
5446 : // applypar.define ("select", incrSel);
5447 0 : applypar.define ("fieldstr", calFields);
5448 0 : applypar.define ("interp", interp);
5449 0 : applypar.define ("spwmap",spwmap);
5450 :
5451 0 : incrcal_->setApply(applypar);
5452 :
5453 0 : Vector<Int> fldidx(0);
5454 0 : if (fields.length()>0)
5455 0 : fldidx=getFieldIdx(fields);
5456 :
5457 : // All ready, now do the accumulation
5458 0 : incal_->accumulate(incrcal_,fldidx);
5459 :
5460 : // ...and store the result
5461 : logSink() << "Storing accumulated calibration in table: "
5462 : << outtab
5463 0 : << LogIO::POST;
5464 :
5465 0 : if (outtab != "")
5466 0 : incal_->calTableName()=outtab;
5467 :
5468 0 : incal_->storeNCT();
5469 :
5470 0 : delete incal_;
5471 0 : delete incrcal_;
5472 :
5473 : logSink() << "Finished accumulation."
5474 0 : << LogIO::POST;
5475 :
5476 0 : } catch (AipsError x) {
5477 : logSink() << LogIO::SEVERE
5478 : << "Caught Exception: "
5479 : << x.getMesg()
5480 0 : << LogIO::POST;
5481 :
5482 0 : if (incal_) delete incal_;
5483 0 : if (incrcal_) delete incrcal_;
5484 :
5485 0 : throw(AipsError("Error in Calibrater::accumulate."));
5486 : return;
5487 0 : }
5488 0 : return;
5489 :
5490 : }
5491 :
5492 0 : void OldCalibrater::specifycal(const String& type,
5493 : const String& caltable,
5494 : const String& time,
5495 : const String& spw,
5496 : const String& antenna,
5497 : const String& pol,
5498 : const Vector<Double>& parameter,
5499 : const String& infile,
5500 : const Bool& uniform) {
5501 :
5502 0 : logSink() << LogOrigin("Calibrater","specifycal") << LogIO::NORMAL;
5503 :
5504 : // SVJ objects:
5505 0 : SolvableVisCal *cal_(NULL);
5506 :
5507 : try {
5508 :
5509 : // Set record format for calibration table application information
5510 0 : RecordDesc specifyDesc;
5511 0 : specifyDesc.addField ("caltable", TpString);
5512 0 : specifyDesc.addField ("time", TpString);
5513 0 : specifyDesc.addField ("spw", TpArrayInt);
5514 0 : specifyDesc.addField ("antenna", TpArrayInt);
5515 0 : specifyDesc.addField ("pol", TpString);
5516 0 : specifyDesc.addField ("parameter", TpArrayDouble);
5517 0 : specifyDesc.addField ("caltype",TpString);
5518 0 : specifyDesc.addField ("infile",TpString);
5519 0 : specifyDesc.addField ("uniform",TpBool);
5520 :
5521 : // Create record with the requisite field values
5522 0 : Record specify(specifyDesc);
5523 0 : specify.define ("caltable", caltable);
5524 0 : specify.define ("time", time);
5525 0 : if (spw=="*")
5526 0 : specify.define ("spw",Vector<Int>(1,-1));
5527 : else
5528 0 : specify.define ("spw",getSpwIdx(spw));
5529 0 : if (antenna=="*")
5530 0 : specify.define ("antenna",Vector<Int>(1,-1) );
5531 : else
5532 0 : specify.define ("antenna",getAntIdx(antenna));
5533 0 : specify.define ("pol",pol);
5534 0 : specify.define ("parameter",parameter);
5535 0 : specify.define ("caltype",type);
5536 0 : specify.define ("infile",infile);
5537 0 : specify.define ("uniform",uniform);
5538 :
5539 : // Now do it
5540 0 : String utype=upcase(type);
5541 0 : if (utype=="G" || utype.contains("AMP") || utype.contains("PH"))
5542 0 : cal_ = createSolvableVisCal("G",*vs_p);
5543 0 : else if (utype=='K' || utype.contains("SBD") || utype.contains("DELAY"))
5544 0 : cal_ = createSolvableVisCal("K",*vs_p);
5545 0 : else if (utype.contains("MBD"))
5546 0 : cal_ = createSolvableVisCal("K",*vs_p); // as of 5.3, KMBD is just K
5547 0 : else if (utype.contains("ANTPOS"))
5548 0 : cal_ = createSolvableVisCal("KANTPOS",*vs_p);
5549 0 : else if (utype.contains("TSYS"))
5550 0 : cal_ = createSolvableVisCal("TSYS",*vs_p);
5551 0 : else if (utype.contains("EVLAGAIN") ||
5552 0 : utype.contains("SWP") ||
5553 0 : utype.contains("RQ"))
5554 0 : cal_ = createSolvableVisCal("EVLASWP",*vs_p);
5555 0 : else if (utype.contains("OPAC"))
5556 0 : cal_ = createSolvableVisCal("TOPAC",*vs_p);
5557 0 : else if (utype.contains("GC") || utype.contains("EFF"))
5558 0 : cal_ = createSolvableVisCal("GAINCURVE",*vs_p);
5559 0 : else if (utype.contains("TEC"))
5560 0 : cal_ = createSolvableVisCal("TEC",*vs_p);
5561 0 : else if (utype.contains("SDSKY_PS"))
5562 0 : cal_ = createSolvableVisCal("SDSKY_PS",*vs_p);
5563 0 : else if (utype.contains("SDSKY_RASTER"))
5564 0 : cal_ = createSolvableVisCal("SDSKY_RASTER",*vs_p);
5565 0 : else if (utype.contains("SDSKY_OTF"))
5566 0 : cal_ = createSolvableVisCal("SDSKY_OTF",*vs_p);
5567 : else
5568 0 : throw(AipsError("Unrecognized caltype."));
5569 :
5570 : // set up for specification (set up the CalSet)
5571 0 : cal_->setSpecify(specify);
5572 :
5573 : // fill with specified values
5574 0 : cal_->specify(specify);
5575 :
5576 : // Store result
5577 0 : cal_->storeNCT();
5578 :
5579 0 : delete cal_;
5580 :
5581 0 : } catch (AipsError x) {
5582 : logSink() << LogIO::SEVERE
5583 : << "Caught Exception: "
5584 : << x.getMesg()
5585 0 : << LogIO::POST;
5586 :
5587 0 : if (cal_) delete cal_;
5588 :
5589 0 : throw(AipsError("Error in Calibrater::specifycal."));
5590 : return;
5591 0 : }
5592 0 : return;
5593 :
5594 : }
5595 :
5596 0 : Bool OldCalibrater::smooth(const String& infile,
5597 : String& outfile, // const Bool& append,
5598 : const String& smoothtype,
5599 : const Double& smoothtime,
5600 : const String& fields,
5601 : const bool ratesmooth)
5602 : {
5603 :
5604 : // TBD: support append?
5605 : // TBD: spw selection?
5606 :
5607 0 : logSink() << LogOrigin("Calibrater","smooth") << LogIO::NORMAL;
5608 :
5609 0 : logSink() << "Beginning smoothing/interpolating method." << LogIO::POST;
5610 :
5611 :
5612 : // A pointer to an SVC
5613 0 : SolvableVisCal *svc(NULL);
5614 :
5615 : try {
5616 :
5617 : // Handle no in file
5618 0 : if (infile=="")
5619 0 : throw(AipsError("Please specify an input calibration table."));
5620 :
5621 : // Handle bad smoothtype
5622 0 : if (smoothtype!="mean" && smoothtype!="median")
5623 0 : throw(AipsError("Unrecognized smooth type!"));
5624 :
5625 : // Handle bad smoothtime
5626 0 : if (smoothtime<=0)
5627 0 : throw(AipsError("Please specify a strictly positive smoothtime."));
5628 :
5629 : // Handle no outfile
5630 0 : if (outfile=="") {
5631 0 : outfile=infile;
5632 : logSink() << "Will overwrite input file with smoothing result."
5633 0 : << LogIO::POST;
5634 : }
5635 :
5636 :
5637 0 : svc = createSolvableVisCal(calTableType(infile),*vs_p);
5638 :
5639 0 : if (svc->smoothable()) {
5640 :
5641 : // Fill calibration table using setApply
5642 0 : RecordDesc applyparDesc;
5643 0 : applyparDesc.addField ("table", TpString);
5644 0 : Record applypar(applyparDesc);
5645 0 : applypar.define ("table", infile);
5646 0 : svc->setApply(applypar);
5647 :
5648 : // Convert refFields/transFields to index lists
5649 0 : Vector<Int> fldidx(0);
5650 0 : if (fields.length()>0)
5651 0 : fldidx=getFieldIdx(fields);
5652 :
5653 : // Delegate to SVC
5654 0 : svc->smooth(fldidx,smoothtype,smoothtime,ratesmooth);
5655 :
5656 : // Store the result on disk
5657 : // if (append) logSink() << "Appending result to " << outfile << LogIO::POST;
5658 : //else
5659 0 : logSink() << "Storing result in " << outfile << LogIO::POST;
5660 :
5661 :
5662 0 : if (outfile != "")
5663 0 : svc->calTableName()=outfile;
5664 0 : svc->storeNCT();
5665 :
5666 : // Clean up
5667 0 : if (svc) delete svc; svc=NULL;
5668 :
5669 : // Apparently, it worked
5670 0 : return true;
5671 :
5672 0 : }
5673 : else
5674 0 : throw(AipsError("This type ("+svc->typeName()+") does not support smoothing."));
5675 :
5676 0 : } catch (AipsError x) {
5677 :
5678 : logSink() << LogIO::SEVERE
5679 : << "Caught Exception: "
5680 : << x.getMesg()
5681 0 : << LogIO::POST;
5682 : // Clean up
5683 0 : if (svc) delete svc; svc=NULL;
5684 :
5685 0 : throw(AipsError("Error in Calibrater::smooth."));
5686 :
5687 : return false;
5688 0 : }
5689 : return false;
5690 : }
5691 :
5692 : // List a calibration table
5693 0 : Bool OldCalibrater::listCal(const String& infile,
5694 : const String& field,
5695 : const String& antenna,
5696 : const String& spw,
5697 : const String& listfile,
5698 : const Int& pagerows) {
5699 :
5700 0 : SolvableVisCal *svc(NULL);
5701 0 : logSink() << LogOrigin("Calibrater","listCal");
5702 :
5703 : try {
5704 :
5705 : // Trap (currently) unsupported types
5706 0 : if (upcase(calTableType(infile))=="GSPLINE" ||
5707 0 : upcase(calTableType(infile))=="BPOLY")
5708 0 : throw(AipsError("GSPLINE and BPOLY tables cannot currently be listed."));
5709 :
5710 : // Get user's selected fields, ants
5711 0 : Vector<Int> ufldids=getFieldIdx(field);
5712 0 : Vector<Int> uantids=getAntIdx(antenna);
5713 :
5714 0 : String newSpw = spw;
5715 0 : Bool defaultSelect = false;
5716 0 : if (spw.empty()) { // list all channels (default)
5717 0 : defaultSelect = true;
5718 0 : newSpw = "*";
5719 : logSink() << LogIO::NORMAL1 << "Spws selected: ALL" << endl
5720 0 : << "Channels selected: ALL" << LogIO::POST;
5721 : }
5722 : // Get user's selected spw and channels
5723 0 : Vector<Int> uspwids=getSpwIdx(newSpw);
5724 0 : Matrix<Int> uchanids=getChanIdx(newSpw);
5725 0 : if (!defaultSelect) {
5726 : logSink() << LogIO::NORMAL1 << "Spw and Channel selection matrix: "
5727 : << endl << "Each rows shows: [ Spw , Start Chan , Stop Chan , Chan Step ]"
5728 0 : << endl << uchanids << LogIO::POST;
5729 : }
5730 : logSink() << LogIO::DEBUG2
5731 : << "uspwids = " << uspwids << endl
5732 0 : << "uchanids = " << uchanids << LogIO::POST;
5733 :
5734 : // By default, do first spw, first chan
5735 0 : if (uspwids.nelements()==0) {
5736 0 : uchanids.resize(1,4);
5737 0 : uchanids=0;
5738 : }
5739 :
5740 : // Set record format for calibration table application information
5741 0 : RecordDesc applyparDesc;
5742 0 : applyparDesc.addField ("table", TpString);
5743 :
5744 : // Create record with the requisite field values
5745 0 : Record applypar(applyparDesc);
5746 0 : applypar.define ("table", infile);
5747 :
5748 : // Generate the VisCal to be listed
5749 0 : svc = createSolvableVisCal(calTableType(infile),*vs_p);
5750 0 : svc->setApply(applypar);
5751 :
5752 : // list it
5753 0 : svc->listCal(ufldids,uantids,uchanids, //uchanids(0,0),uchanids(0,1),
5754 : listfile,pagerows);
5755 :
5756 0 : if (svc) delete svc; svc=NULL;
5757 :
5758 0 : return true;
5759 :
5760 0 : } catch (AipsError x) {
5761 :
5762 : logSink() << LogIO::SEVERE
5763 : << "Caught Exception: "
5764 : << x.getMesg()
5765 0 : << LogIO::POST;
5766 : // Clean up
5767 0 : if (svc) delete svc; svc=NULL;
5768 :
5769 0 : throw(AipsError("Error in Calibrater::listCal."));
5770 :
5771 : return false;
5772 0 : }
5773 : return false;
5774 :
5775 : }
5776 :
5777 0 : Bool OldCalibrater::initialize(MeasurementSet& inputMS,
5778 : Bool compress,
5779 : Bool addScratch, Bool addModel) {
5780 :
5781 0 : logSink() << LogOrigin("Calibrater","") << LogIO::NORMAL3;
5782 :
5783 : try {
5784 0 : timer_p.mark();
5785 :
5786 : // Set pointer ms_p from input MeasurementSet
5787 0 : if (ms_p) {
5788 0 : *ms_p=inputMS;
5789 : } else {
5790 0 : ms_p = new MeasurementSet(inputMS);
5791 0 : AlwaysAssert(ms_p,AipsError);
5792 : };
5793 :
5794 : // Disabled 2016/04/19 (gmoellen): avoid direct MS.HISTORY
5795 : // updates from below the python level, FOR NOW
5796 :
5797 : /*
5798 :
5799 : // Setup to write LogIO to HISTORY Table in MS
5800 : if(!(Table::isReadable(ms_p->historyTableName()))){
5801 : // create a new HISTORY table if its not there
5802 : TableRecord &kws = ms_p->rwKeywordSet();
5803 : SetupNewTable historySetup(ms_p->historyTableName(),
5804 : MSHistory::requiredTableDesc(),Table::New);
5805 : kws.defineTable(MS::keywordName(MS::HISTORY), Table(historySetup));
5806 : MSHistoryHandler::addMessage(*ms_p, "HISTORY Table added by Calibrater",
5807 : "Calibrater","","Calibrater::initialize()");
5808 : }
5809 : historytab_p=Table(ms_p->historyTableName(),
5810 : TableLock(TableLock::UserNoReadLocking), Table::Update);
5811 : // jagonzal (CAS-4110): When the selectvis method throws an exception the initialize method
5812 : // is called again to leave the calibrater in a proper state, and since there was a previous
5813 : // initialization the history handler was already created, and has to be destroyed before
5814 : // creating a new one to avoid leaveing the HISTORY table opened.
5815 : if (hist_p) delete hist_p;
5816 : hist_p= new MSHistoryHandler(*ms_p, "calibrater");
5817 :
5818 : // (2016/04/19) */
5819 :
5820 : // Remember the ms's name
5821 0 : msname_p=ms_p->tableName();
5822 :
5823 :
5824 : // Add/init scr cols, if requested (init is hard-wired)
5825 0 : if (addScratch || addModel) {
5826 0 : Bool alsoinit=true;
5827 0 : VisSetUtil::addScrCols(*ms_p,addModel,addScratch,alsoinit,compress);
5828 : }
5829 :
5830 : // Set the selected MeasurementSet to be the same initially
5831 : // as the input MeasurementSet
5832 0 : if (mssel_p) delete mssel_p;
5833 0 : mssel_p=new MeasurementSet(*ms_p);
5834 :
5835 : logSink() << LogIO::NORMAL
5836 : << "Initializing nominal selection to the whole MS."
5837 0 : << LogIO::POST;
5838 :
5839 :
5840 : // Create a VisSet with no selection
5841 : // (gmoellen 2012/02/06: this merely makes a VisIter now)
5842 0 : if (vs_p) {
5843 0 : delete vs_p;
5844 0 : vs_p=0;
5845 : };
5846 0 : Block<Int> nosort(0);
5847 0 : Matrix<Int> noselection;
5848 0 : Double timeInterval=0;
5849 : // gmoellen 2012/02/06 vs_p=new VisSet(*ms_p,nosort,noselection,addScratch,timeInterval,compress, addModel);
5850 0 : vs_p=new VisSet(*ms_p,nosort,noselection,false,timeInterval,false,false);
5851 :
5852 : // Size-up the chanmask PB
5853 0 : initChanMask();
5854 :
5855 : // Create the associated VisEquation
5856 : // TBD: move to ctor and make it non-pointer
5857 0 : if (ve_p) {
5858 0 : delete ve_p;
5859 0 : ve_p=0;
5860 : };
5861 0 : ve_p=new VisEquation();
5862 :
5863 : // Reset the apply/solve VisCals
5864 0 : reset(true,true);
5865 :
5866 0 : return true;
5867 :
5868 0 : } catch (AipsError x) {
5869 0 : logSink() << LogOrigin("Calibrater","initialize",WHERE)
5870 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
5871 0 : << LogIO::POST;
5872 0 : cleanup();
5873 0 : cleanupVisSet();
5874 0 : if (ms_p) delete ms_p; ms_p=NULL;
5875 0 : if (hist_p) delete hist_p; hist_p=NULL;
5876 :
5877 0 : throw(AipsError("Error in Calibrater::initialize()"));
5878 : return false;
5879 0 : }
5880 : return false;
5881 : }
5882 :
5883 0 : Bool OldCalibrater::initCalSet(const Int& calSet)
5884 : {
5885 :
5886 : // logSink() << LogOrigin("Calibrater","initCalSet") << LogIO::NORMAL3;
5887 :
5888 0 : if (vs_p) {
5889 :
5890 0 : Block<Int> columns;
5891 : // include scan iteration, for more optimal iteration
5892 0 : columns.resize(5);
5893 0 : columns[0]=MS::ARRAY_ID;
5894 0 : columns[1]=MS::SCAN_NUMBER;
5895 0 : columns[2]=MS::FIELD_ID;
5896 0 : columns[3]=MS::DATA_DESC_ID;
5897 0 : columns[4]=MS::TIME;
5898 0 : vs_p->resetVisIter(columns,0.0);
5899 :
5900 0 : vs_p->initCalSet(calSet);
5901 0 : return true;
5902 0 : }
5903 : else {
5904 0 : throw(AipsError("Calibrater cannot initCalSet"));
5905 : return false;
5906 : }
5907 : }
5908 :
5909 0 : Bool OldCalibrater::cleanupVisSet() {
5910 :
5911 : // logSink() << LogOrigin("OldCalibrater","cleanupVisSet") << LogIO::NORMAL;
5912 :
5913 0 : if(vs_p) delete vs_p; vs_p=0;
5914 :
5915 : // Delete chanmask
5916 0 : initChanMask();
5917 :
5918 0 : return true;
5919 :
5920 : }
5921 :
5922 0 : VisibilityIterator::DataColumn OldCalibrater::configureForCorrection ()
5923 : {
5924 0 : if (!ok())
5925 0 : throw(AipsError("Calibrater not prepared for correct!"));
5926 :
5927 : // Nominally, we write out to the CORRECTED_DATA, unless absent
5928 0 : VisibilityIterator::DataColumn whichOutCol(VisibilityIterator::Corrected);
5929 :
5930 0 : if (!ms_p->tableDesc().isColumn("CORRECTED_DATA"))
5931 0 : throw(AipsError("CORRECTED_DATA column unexpectedly absent. Cannot correct."));
5932 :
5933 : // Ensure apply list non-zero and properly sorted
5934 0 : ve_p->setapply(vc_p);
5935 :
5936 : // Report the types that will be applied
5937 0 : applystate();
5938 :
5939 : // Arrange for iteration over data
5940 0 : Block<Int> columns;
5941 : // include scan iteration
5942 0 : columns.resize(5);
5943 0 : columns[0]=MS::ARRAY_ID;
5944 0 : columns[1]=MS::SCAN_NUMBER;
5945 0 : columns[2]=MS::FIELD_ID;
5946 0 : columns[3]=MS::DATA_DESC_ID;
5947 0 : columns[4]=MS::TIME;
5948 :
5949 : // Reset the VisibilityIterator in the VisSet.
5950 0 : vs_p->resetVisIter (columns, 0.0);
5951 :
5952 0 : return whichOutCol;
5953 0 : }
5954 :
5955 0 : void OldCalibrater::selectChannel(const String& spw) {
5956 :
5957 : // Initialize the chanmask_
5958 0 : initChanMask();
5959 :
5960 0 : if (mss_p && mssel_p) {
5961 :
5962 : // Refresh the frequencySelections object to feed to VI2, if relevant
5963 0 : frequencySelections_p.reset(new vi::FrequencySelections());
5964 :
5965 0 : vi::FrequencySelectionUsingChannels usingChannels;
5966 0 : usingChannels.add(*mss_p,mssel_p);
5967 0 : frequencySelections_p->add(usingChannels);
5968 :
5969 : // cout << usingChannels.toString() << endl;
5970 : // cout << "FS.size() = " << frequencySelections_p->size() << endl;
5971 :
5972 0 : }
5973 :
5974 0 : Matrix<Int> chansel = getChanIdx(spw);
5975 0 : uInt nselspw=chansel.nrow();
5976 :
5977 0 : if (nselspw==0)
5978 : logSink() << "Frequency selection: Selecting all channels in all spws."
5979 0 : << LogIO::POST;
5980 : else {
5981 :
5982 0 : logSink() << "Frequency selection: " << LogIO::POST;
5983 :
5984 : // Trap non-unit step (for now)
5985 0 : if (ntrue(chansel.column(3)==1)!=nselspw) {
5986 : logSink() << LogIO::WARN
5987 : << "Calibration does not support non-unit channel stepping; "
5988 : << "using step=1."
5989 0 : << LogIO::POST;
5990 0 : chansel.column(3)=1;
5991 : }
5992 :
5993 0 : Int nspw=vs_p->numberSpw();
5994 0 : Vector<Int> nChan0;
5995 0 : nChan0 = vs_p->numberChan();
5996 :
5997 0 : Vector<Int> uspw(chansel.column(0));
5998 0 : Vector<Int> ustart(chansel.column(1));
5999 0 : Vector<Int> uend(chansel.column(2));
6000 :
6001 0 : Vector<Int> start(nspw,INT_MAX);
6002 0 : Vector<Int> end(nspw,-INT_MAX);
6003 0 : logSink() << LogIO::NORMAL;
6004 0 : for (uInt i=0;i<nselspw;++i) {
6005 :
6006 0 : Int& spw=uspw(i);
6007 :
6008 : // Initialize this spw mask, if necessary (def = masked)
6009 0 : if (!chanmask_[spw])
6010 0 : chanmask_[spw]=new Vector<Bool>(nChan0(spw),true);
6011 :
6012 : // revise net start/end/nchan
6013 0 : start(spw)=min(start(spw),ustart(i));
6014 0 : end(spw)=max(end(spw),uend(i));
6015 0 : Int nchan=end(spw)-start(spw)+1; // net inclusive nchan
6016 :
6017 : // User's
6018 0 : Int step=chansel(i,3);
6019 0 : Int unchan=uend(i)-ustart(i)+1;
6020 :
6021 : // Update the mask (false = valid)
6022 0 : (*chanmask_[spw])(Slice(ustart(i),unchan))=false;
6023 :
6024 :
6025 : logSink() << ". Spw " << spw << ":"
6026 0 : << ustart(i) << "~" << uend(i)
6027 0 : << " (" << uend(i)-ustart(i)+1 << " channels,"
6028 : << " step by " << step << ")"
6029 0 : << endl;
6030 :
6031 : /*
6032 : cout << i << " " << spw << " {"
6033 : << start(spw) << " [" << ustart(i) << " "
6034 : << uend(i) << "] " << end(spw) << "}" << endl;
6035 : cout << "chanmask = ";
6036 : for (Int j=0;j<nChan0(spw);++j) cout << (*chanmask_[spw])(j);
6037 : cout << endl << endl;
6038 : */
6039 :
6040 : // Call via VisSet (avoid call to VisIter::origin)
6041 0 : vs_p->selectChannel(1,start(spw),nchan,step,spw,false);
6042 :
6043 : } // i
6044 0 : logSink() << LogIO::POST;
6045 :
6046 0 : } // non-triv spw selection
6047 :
6048 : // For testing:
6049 : if (false) {
6050 :
6051 : VisIter& vi(vs_p->iter());
6052 : VisBuffer vb(vi);
6053 :
6054 : // Pass each timestamp (VisBuffer) to VisEquation for correction
6055 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
6056 : vi.origin();
6057 : // for (vi.origin(); vi.more(); vi++)
6058 : cout << vb.spectralWindow() << " "
6059 : << vb.nChannel() << " "
6060 : << vb.channel() << " "
6061 : << vb.visCube().shape()
6062 : << endl;
6063 : }
6064 : }
6065 :
6066 0 : }
6067 :
6068 0 : void OldCalibrater::initChanMask() {
6069 :
6070 0 : for (uInt i=0;i<chanmask_.nelements();++i)
6071 0 : if (chanmask_[i])
6072 0 : delete chanmask_[i];
6073 0 : if (vs_p) {
6074 0 : chanmask_.resize(vs_p->numberSpw(),true);
6075 0 : chanmask_=NULL;
6076 : }
6077 : else {
6078 : // throw(AipsError("Trouble sizing chanmask!"));
6079 : // just don't support channel masking:
6080 0 : chanmask_.resize(0,true);
6081 : }
6082 :
6083 0 : }
6084 :
6085 : // Select on channel in the VisSet
6086 0 : void OldCalibrater::selectChannel(const String& mode,
6087 : const Int& nchan,
6088 : const Int& start, const Int& step,
6089 : const MRadialVelocity& mStart,
6090 : const MRadialVelocity& mStep) {
6091 :
6092 : // Set data selection variables
6093 0 : dataMode_p=mode;
6094 0 : dataNchan_p=nchan;
6095 0 : if (dataNchan_p<0) dataNchan_p=0;
6096 0 : dataStart_p=start;
6097 0 : if (dataStart_p<0) dataNchan_p=0;
6098 0 : dataStep_p=step;
6099 0 : if (dataStep_p<1) dataNchan_p=1;
6100 :
6101 0 : mDataStart_p=mStart;
6102 0 : mDataStep_p=mStep;
6103 :
6104 : // Select on frequency channel
6105 0 : if(dataMode_p=="channel") {
6106 : // *** this bit here is temporary till we unifomize data selection
6107 : //Getting the selected SPWs
6108 0 : MSMainColumns msc(*mssel_p);
6109 0 : Vector<Int> dataDescID = msc.dataDescId().getColumn();
6110 : Bool dum;
6111 0 : Sort sort( dataDescID.getStorage(dum),sizeof(Int) );
6112 0 : sort.sortKey((uInt)0,TpInt);
6113 0 : Vector<uInt> index,uniq;
6114 0 : sort.sort(index,dataDescID.nelements());
6115 0 : uInt nSpw = sort.unique(uniq,index);
6116 :
6117 0 : Vector<Int> selectedSpw(nSpw);
6118 0 : Vector<Int> nChan(nSpw);
6119 0 : for (uInt k=0; k < nSpw; ++k) {
6120 0 : selectedSpw[k]=dataDescID[index[uniq[k]]];
6121 0 : nChan[k]=vs_p->numberChan()(selectedSpw[k]);
6122 :
6123 : }
6124 0 : if(dataNchan_p==0) dataNchan_p=vs_p->numberChan()(selectedSpw[0]);
6125 0 : if(dataStart_p<0) {
6126 : logSink() << LogIO::SEVERE << "Illegal start pixel = "
6127 0 : << dataStart_p << LogIO::POST;
6128 : }
6129 0 : Int end = Int(dataStart_p) + Int(dataNchan_p) * Int(dataStep_p);
6130 0 : for (uInt k=0; k < selectedSpw.nelements() ; ++k){
6131 0 : if(end < 1 || end > nChan[k]) {
6132 : logSink() << LogIO::SEVERE << "Illegal step pixel = " << dataStep_p
6133 0 : << " in Spw " << selectedSpw[k]
6134 0 : << LogIO::POST;
6135 : }
6136 : logSink() << "Selecting "<< dataNchan_p
6137 : << " channels, starting at visibility channel "
6138 : << dataStart_p << " stepped by "
6139 0 : << dataStep_p << " in Spw " << selectedSpw[k] << LogIO::POST;
6140 :
6141 : // Set frequency channel selection for all spectral window id's
6142 : Int nch;
6143 : //Vector<Int> nChan=vs_p->numberChan();
6144 : //Int nSpw=vs_p->numberSpw();
6145 0 : if (dataNchan_p==0) {
6146 0 : nch=nChan(k);
6147 : }else {
6148 0 : nch=dataNchan_p;
6149 : };
6150 0 : vs_p->selectChannel(1,dataStart_p,nch,dataStep_p,selectedSpw[k]);
6151 :
6152 : }
6153 0 : }
6154 : // Select on velocity
6155 0 : else if (dataMode_p=="velocity") {
6156 0 : MVRadialVelocity mvStart(mDataStart_p.get("m/s"));
6157 0 : MVRadialVelocity mvStep(mDataStep_p.get("m/s"));
6158 : MRadialVelocity::Types
6159 0 : vType((MRadialVelocity::Types)mDataStart_p.getRefPtr()->getType());
6160 : logSink() << "Selecting "<< dataNchan_p
6161 : << " channels, starting at radio velocity " << mvStart
6162 : << " stepped by " << mvStep << ", reference frame is "
6163 0 : << MRadialVelocity::showType(vType) << LogIO::POST;
6164 0 : vs_p->iter().selectVelocity(Int(dataNchan_p), mvStart, mvStep,
6165 : vType, MDoppler::RADIO);
6166 0 : }
6167 :
6168 : // Select on optical velocity
6169 0 : else if (dataMode_p=="opticalvelocity") {
6170 0 : MVRadialVelocity mvStart(mDataStart_p.get("m/s"));
6171 0 : MVRadialVelocity mvStep(mDataStep_p.get("m/s"));
6172 : MRadialVelocity::Types
6173 0 : vType((MRadialVelocity::Types)mDataStart_p.getRefPtr()->getType());
6174 : logSink() << "Selecting "<< dataNchan_p
6175 : << " channels, starting at optical velocity " << mvStart
6176 : << " stepped by " << mvStep << ", reference frame is "
6177 0 : << MRadialVelocity::showType(vType) << LogIO::POST;
6178 0 : vs_p->iter().selectVelocity(Int(dataNchan_p), mvStart, mvStep,
6179 : vType, MDoppler::OPTICAL);
6180 0 : }
6181 :
6182 :
6183 0 : }
6184 :
6185 0 : Bool OldCalibrater::ok() {
6186 :
6187 0 : if(vs_p && ms_p && mssel_p && ve_p) {
6188 0 : return true;
6189 : }
6190 : else {
6191 0 : logSink() << "Calibrater is not yet initialized" << LogIO::POST;
6192 0 : return false;
6193 : }
6194 : }
6195 :
6196 0 : Bool OldCalibrater::genericGatherAndSolve() {
6197 :
6198 : #ifdef _OPENMP
6199 0 : Double Tsetup(0.0),Tgather(0.0),Tsolve(0.0);
6200 0 : Double time0=omp_get_wtime();
6201 : #endif
6202 :
6203 : //cout << "Generic gather and solve." << endl;
6204 :
6205 : // Create the solver
6206 0 : VisCalSolver vcs;
6207 :
6208 : // Inform logger/history
6209 0 : logSink() << "Solving for " << svc_p->typeName()
6210 0 : << LogIO::POST;
6211 :
6212 : // Initialize the svc according to current VisSet
6213 : // (this counts intervals, sizes CalSet)
6214 0 : Vector<Int> nChunkPerSol;
6215 0 : Int nSol = svc_p->sizeUpSolve(*vs_p,nChunkPerSol);
6216 :
6217 : // Create the in-memory (New)CalTable
6218 0 : svc_p->createMemCalTable();
6219 :
6220 : // The iterator, VisBuffer
6221 0 : VisIter& vi(vs_p->iter());
6222 0 : VisBuffer vb(vi);
6223 :
6224 0 : Vector<Int> slotidx(vs_p->numberSpw(),-1);
6225 :
6226 : // We will remember which spws couldn't be processed
6227 0 : Vector<Bool> unsolspw(vi.numberSpw());
6228 0 : unsolspw.set(false);
6229 :
6230 : // Manage verbosity of partial channel averaging
6231 0 : Vector<Bool> verb(vi.numberSpw(),true);
6232 :
6233 0 : Vector<Int64> nexp(vi.numberSpw(),0), natt(vi.numberSpw(),0),nsuc(vi.numberSpw(),0);
6234 :
6235 : #ifdef _OPENMP
6236 0 : Tsetup+=(omp_get_wtime()-time0);
6237 : #endif
6238 :
6239 0 : Int nGood(0);
6240 0 : vi.originChunks();
6241 0 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
6242 :
6243 : #ifdef _OPENMP
6244 0 : time0=omp_get_wtime();
6245 : #endif
6246 :
6247 0 : nexp(vi.spectralWindow())+=1;
6248 :
6249 : // capture obs, scan info so we can set it later
6250 : // (and not rely on what the VB averaging code can't properly do)
6251 0 : Vector<Int> scv,obsv;
6252 0 : Int solscan=vi.scan(scv)(0),solobs=vi.observationId(obsv)(0);
6253 :
6254 : // Arrange to accumulate
6255 : // VisBuffAccumulator vba(vs_p->numberAnt(),svc_p->preavg(),false);
6256 0 : VisBuffGroupAcc vbga(vs_p->numberAnt(),vs_p->numberSpw(),vs_p->numberFld(),svc_p->preavg());
6257 :
6258 0 : for (Int ichunk=0;ichunk<nChunkPerSol(isol);++ichunk) {
6259 :
6260 : // Current _chunk_'s spw
6261 0 : Int spw(vi.spectralWindow());
6262 :
6263 : // Only accumulate for solve if we can pre-calibrate
6264 0 : if (ve_p->spwOK(spw)) {
6265 :
6266 : // Collapse each timestamp in this chunk according to VisEq
6267 : // with calibration and averaging
6268 0 : for (vi.origin(); vi.more(); vi++) {
6269 :
6270 : // Force read of the field Id
6271 0 : vb.fieldId();
6272 :
6273 : // Apply the channel mask (~no-op, if unnecessary)
6274 0 : svc_p->applyChanMask(vb);
6275 :
6276 : // This forces the data/model/wt I/O, and applies
6277 : // any prior calibrations
6278 0 : ve_p->collapse(vb);
6279 :
6280 : // If permitted/required by solvable component, normalize
6281 0 : if (svc_p->normalizable())
6282 0 : vb.normalize();
6283 :
6284 : // If this solve not freqdep, and channels not averaged yet, do so
6285 0 : if (!svc_p->freqDepMat() && vb.nChannel()>1)
6286 0 : vb.freqAveCubes();
6287 :
6288 0 : if (svc_p->freqDepPar() &&
6289 0 : svc_p->fsolint()!="none" &&
6290 0 : svc_p->fintervalCh()>0.0) {
6291 : // cout << "svc_p->currSpw() = " << svc_p->currSpw() << endl;
6292 0 : if (verb(spw))
6293 : logSink() << " Reducing nchan in spw "
6294 : << spw
6295 0 : << " from " << vb.nChannel();
6296 0 : vb.channelAve(svc_p->chanAveBounds(spw));
6297 :
6298 : // Kludge for 3.4 to reset corr-indep flag to correct channel axis shape
6299 : // (because we use vb.flag() below, rather than vb.flagCube())
6300 0 : vb.flag().assign(operator>(partialNTrue(vb.flagCube(),IPosition(1,0)),0UL));
6301 :
6302 0 : if (verb(spw)) {
6303 : logSink() << " to "
6304 0 : << vb.nChannel() << LogIO::POST;
6305 0 : verb(spw)=false; // suppress future verbosity in this spw
6306 : }
6307 : }
6308 :
6309 : // Accumulate collapsed vb in a time average
6310 : // (only if the vb contains any unflagged data)
6311 0 : if (nfalse(vb.flag())>0)
6312 0 : vbga.accumulate(vb);
6313 :
6314 : }
6315 : }
6316 : else
6317 : // This spw not accumulated for solve
6318 0 : unsolspw(spw)=true;
6319 :
6320 : // Advance the VisIter, if possible
6321 0 : if (vi.moreChunks()) vi.nextChunk();
6322 :
6323 : }
6324 :
6325 : // Finalize the averged VisBuffer
6326 0 : vbga.finalizeAverage();
6327 :
6328 : // Establish meta-data for this interval
6329 : // (some of this may be used _during_ solve)
6330 : // (this sets currSpw() in the SVC)
6331 0 : Bool vbOk=(vbga.nBuf()>0 && svc_p->syncSolveMeta(vbga));
6332 :
6333 0 : svc_p->overrideObsScan(solobs,solscan);
6334 :
6335 : #ifdef _OPENMP
6336 0 : Tgather+=(omp_get_wtime()-time0);
6337 0 : time0=omp_get_wtime();
6338 : #endif
6339 :
6340 0 : if (vbOk) {
6341 :
6342 : // Use spw of first VB in vbga
6343 : // TBD: (currSpw==thisSpw) here?? (I.e., use svc_p->currSpw()? currSpw is prot!)
6344 0 : Int thisSpw=svc_p->spwMap()(vbga(0).spectralWindow());
6345 :
6346 0 : natt(thisSpw)+=1;
6347 :
6348 0 : slotidx(thisSpw)++;
6349 :
6350 : // Make data amp- or phase-only, if needed
6351 0 : vbga.enforceAPonData(svc_p->apmode());
6352 :
6353 : // Select on correlation via weights, according to the svc
6354 0 : vbga.enforceSolveCorrWeights(svc_p->phandonly());
6355 :
6356 0 : if (svc_p->useGenericSolveOne()) {
6357 : // generic individual solve
6358 :
6359 : //cout << "Generic individual solve: isol=" << isol << endl;
6360 :
6361 : // First guess
6362 0 : svc_p->guessPar(vbga(0));
6363 :
6364 : // Solve for each parameter channel (in curr Spw)
6365 :
6366 : // (NB: force const version of nChanPar() [why?])
6367 : // for (Int ich=0;ich<((const SolvableVisCal*)svc_p)->nChanPar();++ich) {
6368 0 : Bool totalGoodSol(false);
6369 0 : for (Int ich=((const SolvableVisCal*)svc_p)->nChanPar()-1;ich>-1;--ich) {
6370 : // for (Int ich=0;ich<((const SolvableVisCal*)svc_p)->nChanPar();++ich) {
6371 :
6372 : // If pars chan-dep, SVC mechanisms for only one channel at a time
6373 0 : svc_p->markTimer();
6374 0 : svc_p->focusChan()=ich;
6375 :
6376 : // svc_p->state();
6377 :
6378 : // cout << "Starting solution..." << endl;
6379 :
6380 : // Pass VE, SVC, VB to solver
6381 0 : Bool goodSoln=vcs.solve(*ve_p,*svc_p,vbga);
6382 :
6383 : // cout << "goodSoln= " << boolalpha << goodSoln << endl;
6384 :
6385 : // If good...
6386 0 : if (goodSoln) {
6387 0 : totalGoodSol=true;
6388 :
6389 0 : svc_p->formSolveSNR();
6390 0 : svc_p->applySNRThreshold();
6391 :
6392 : // ..and file this solution in the correct slot
6393 0 : if (svc_p->freqDepPar())
6394 0 : svc_p->keep1(ich);
6395 : // svc_p->keep(slotidx(thisSpw));
6396 : // Int n=svc_p->nSlots(thisSpw);
6397 : // svc_p->printActivity(n,slotidx(thisSpw),vi.fieldId(),thisSpw,nGood);
6398 :
6399 : }
6400 : else {
6401 : // report where this failure occured
6402 0 : svc_p->currMetaNote();
6403 0 : if (svc_p->freqDepPar())
6404 : // We must record a flagged solution for this channel
6405 0 : svc_p->keep1(ich);
6406 : }
6407 :
6408 : } // parameter channels
6409 :
6410 0 : if (totalGoodSol) {
6411 0 : svc_p->keepNCT();
6412 0 : nsuc(thisSpw)+=1;
6413 : }
6414 :
6415 : // Count good solutions.
6416 0 : if (totalGoodSol) nGood++;
6417 :
6418 : }
6419 : else {
6420 : //cout << "Self-directed individual solve: isol=" << isol << endl;
6421 : // self-directed individual solve
6422 : // TBD: selfSolveOne should return T/F for "good"
6423 0 : svc_p->selfSolveOne(vbga);
6424 :
6425 : // File this solution in the correct slot of the CalSet
6426 0 : svc_p->keepNCT();
6427 :
6428 0 : nGood++;
6429 0 : nsuc(thisSpw)+=1;
6430 : }
6431 :
6432 : } // vbOK
6433 :
6434 : #ifdef _OPENMP
6435 0 : Tsolve+=(omp_get_wtime()-time0);
6436 : #endif
6437 :
6438 0 : } // isol
6439 :
6440 : logSink() << " Found good "
6441 0 : << svc_p->typeName() << " solutions in "
6442 : << nGood << " slots."
6443 0 : << LogIO::POST;
6444 : #ifdef _OPENMP
6445 : #ifdef REPORT_CAL_TIMING
6446 : cout << "OldCalibrater::genericGatherAndSolve Timing: " << endl;
6447 : cout << " setup=" << Tsetup
6448 : << " gather=" << Tgather
6449 : << " solve=" << Tsolve
6450 : << " total=" << Tsetup+Tgather+Tsolve
6451 : << endl;
6452 : #endif
6453 : #endif
6454 :
6455 0 : summarize_uncalspws(unsolspw, "solv");
6456 :
6457 : // Store whole of result in a caltable
6458 0 : if (nGood==0) {
6459 : logSink() << "No output calibration table written."
6460 0 : << LogIO::POST;
6461 : }
6462 : else {
6463 :
6464 : // TBD: Remove BPOLY specificity here
6465 0 : if (svc_p->typeName()!="BPOLY") {
6466 : // Do global post-solve tinkering (e.g., phase-only, normalization, etc.)
6467 0 : svc_p->globalPostSolveTinker();
6468 :
6469 : // write the table
6470 0 : svc_p->storeNCT();
6471 :
6472 : }
6473 : }
6474 :
6475 : // Fill activity record
6476 0 : actRec_=Record();
6477 0 : actRec_.define("origin","Calibrater::genericGatherAndSolve");
6478 0 : actRec_.define("nExpected",nexp);
6479 0 : actRec_.define("nAttempt",natt);
6480 0 : actRec_.define("nSucceed",nsuc);
6481 :
6482 : {
6483 0 : Record solveRec=svc_p->solveActionRec();
6484 0 : if (solveRec.nfields()>0)
6485 0 : actRec_.merge(solveRec);
6486 0 : }
6487 :
6488 :
6489 :
6490 0 : return true;
6491 :
6492 0 : }
6493 :
6494 :
6495 : } //# NAMESPACE CASA - END
|
