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 0 : CalCounts::CalCounts():
91 0 : antennaMap_(),
92 0 : spwMap_(),
93 0 : totalMap_()
94 0 : {}
95 :
96 0 : void CalCounts::initCounts(Int NSpw, Int NAnt, Int NPol) {
97 0 : Vector<Int> polShape(NPol, 0);
98 0 : Vector<Int> singleVector(1, 0);
99 :
100 0 : totalMap_["expected"] = polShape;
101 0 : totalMap_["data_unflagged"] = polShape;
102 0 : totalMap_["above_minblperant"] = polShape;
103 0 : totalMap_["above_minsnr"] = polShape;
104 :
105 0 : for (Int spw=0; spw<NSpw; spw++) {
106 : // Init spw map keys
107 0 : spwMap_[spw]["expected"] = polShape;
108 0 : spwMap_[spw]["data_unflagged"] = polShape;
109 0 : spwMap_[spw]["above_minblperant"] = polShape;
110 0 : spwMap_[spw]["above_minsnr"] = polShape;
111 0 : for (Int ant=0; ant<NAnt; ant++) {
112 : // Init antenna map keys
113 0 : antennaMap_[spw][ant]["expected"] = polShape;
114 0 : antennaMap_[spw][ant]["data_unflagged"] = polShape;
115 0 : antennaMap_[spw][ant]["above_minblperant"] = polShape;
116 0 : antennaMap_[spw][ant]["above_minsnr"] = polShape;
117 0 : antennaMap_[spw][ant]["used_as_refant"] = singleVector;
118 : }
119 : }
120 0 : }
121 :
122 0 : 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 0 : Vector<Int> spwExp(NPol, 0), spwUnflag(NPol, 0), spwMinsnr(NPol, 0);
125 :
126 0 : for (Int ant=0; ant<NAnt; ant++) {
127 0 : if (resultMap.find(ant) == resultMap.end()) {continue;};
128 : // add if used as refant
129 0 : antennaMap_[spw][ant]["used_as_refant"][0] += resultMap[ant]["used_as_refant"][0];
130 0 : for (Int pol=0; pol<NPol; pol++) {
131 0 : antennaMap_[spw][ant]["expected"][pol] += resultMap[ant]["expected"][pol];
132 0 : antennaMap_[spw][ant]["data_unflagged"][pol] += resultMap[ant]["data_unflagged"][pol];
133 0 : antennaMap_[spw][ant]["above_minblperant"][pol] += resultMap[ant]["above_minblperant"][pol];
134 0 : antennaMap_[spw][ant]["above_minsnr"][pol] += resultMap[ant]["above_minsnr"][pol];
135 : // if any antenna value is good add to spw total
136 0 : if (resultMap[ant]["expected"][pol] == 1) {spwExp[pol] = 1;};
137 0 : if (resultMap[ant]["data_unflagged"][pol] == 1) {spwUnflag[pol] = 1;};
138 0 : if (resultMap[ant]["above_minblperant"][pol] == 1) {spwMinsnr[pol] = 1;};
139 0 : if (resultMap[ant]["above_minsnr"][pol] == 1) {spwMinsnr[pol] = 1;};
140 : }
141 : }
142 : // Add the totals to the spw Map
143 0 : for (Int pol=0; pol<NPol; pol++) {
144 0 : spwMap_[spw]["expected"][pol] += spwExp[pol];
145 0 : spwMap_[spw]["data_unflagged"][pol] += spwUnflag[pol];
146 0 : spwMap_[spw]["above_minblperant"][pol] += spwMinsnr[pol];
147 0 : spwMap_[spw]["above_minsnr"][pol] += spwMinsnr[pol];
148 : // Acumulate total overall values
149 0 : totalMap_["expected"][pol] += spwExp[pol];
150 0 : totalMap_["data_unflagged"][pol] += spwUnflag[pol];
151 0 : totalMap_["above_minblperant"][pol] += spwMinsnr[pol];
152 0 : totalMap_["above_minsnr"][pol] += spwMinsnr[pol];
153 : }
154 :
155 0 : }
156 :
157 0 : void CalCounts::updateRefants(Int NSpw, Int NAnt, std::map<Int, std::map<Int, Int>> refantMap) {
158 0 : for (Int spw=0; spw<NSpw; spw++) {
159 0 : for (Int ant=0; ant<NAnt; ant++) {
160 0 : antennaMap_[spw][ant]["used_as_refant"][0] += refantMap[spw][ant];
161 : }
162 : }
163 0 : }
164 :
165 0 : Vector<Int> CalCounts::antMapVal(Int spw, Int ant, String gate) {
166 0 : return antennaMap_[spw][ant][gate];
167 : }
168 :
169 0 : Vector<Int> CalCounts::spwMapVal(Int spw, String gate) {
170 0 : return spwMap_[spw][gate];
171 : }
172 :
173 0 : Vector<Int> CalCounts::totalMapVal(String gate) {
174 0 : return totalMap_[gate];
175 : }
176 :
177 0 : Record CalCounts::makeRecord(Int NAnt, Int NPol) {
178 :
179 0 : Vector<Int> totExp(NPol, 0), totUnflag(NPol, 0), totMinsnr(NPol, 0);
180 :
181 0 : nSpw = spwMap_.size();
182 0 : Record containerRec = Record();
183 0 : Record resultRec = Record();
184 :
185 0 : resultRec.define("expected", totalMap_["expected"]);
186 0 : resultRec.define("data_unflagged", totalMap_["data_unflagged"]);
187 0 : resultRec.define("above_minblperant", totalMap_["above_minblperant"]);
188 0 : resultRec.define("above_minsnr", totalMap_["above_minsnr"]);
189 :
190 0 : for (Int spw=0; spw<nSpw; spw++) {
191 0 : Record spwRec = Record();
192 0 : spwRec.define("expected", spwMap_[spw]["expected"]);
193 0 : spwRec.define("data_unflagged", spwMap_[spw]["data_unflagged"]);
194 0 : spwRec.define("above_minblperant", spwMap_[spw]["above_minblperant"]);
195 0 : spwRec.define("above_minsnr", spwMap_[spw]["above_minsnr"]);
196 0 : for (Int ant=0; ant<NAnt; ant++) {
197 0 : Record subRec = Record();
198 0 : subRec.define("expected", antennaMap_[spw][ant]["expected"]);
199 0 : subRec.define("data_unflagged", antennaMap_[spw][ant]["data_unflagged"]);
200 0 : subRec.define("above_minblperant", antennaMap_[spw][ant]["above_minblperant"]);
201 0 : subRec.define("above_minsnr", antennaMap_[spw][ant]["above_minsnr"]);
202 0 : subRec.define("used_as_refant", antennaMap_[spw][ant]["used_as_refant"]);
203 0 : spwRec.defineRecord(RecordFieldId("ant"+to_string(ant)), subRec);
204 0 : }
205 0 : resultRec.defineRecord(RecordFieldId("spw"+to_string(spw)), spwRec);
206 0 : }
207 0 : containerRec.defineRecord(RecordFieldId("solvestats"), resultRec);
208 :
209 0 : return containerRec;
210 0 : }
211 :
212 0 : CalCounts::~CalCounts() {}
213 :
214 11 : Calibrater::Calibrater():
215 11 : ms_p(0),
216 11 : mssel_p(0),
217 11 : mss_p(0),
218 11 : frequencySelections_p(nullptr),
219 11 : msmc_p(0),
220 11 : ve_p(0),
221 11 : vc_p(),
222 11 : svc_p(0),
223 11 : histLockCounter_p(),
224 11 : hist_p(0),
225 11 : usingCalLibrary_(false),
226 11 : corrDepFlags_(false), // default (==traditional behavior)
227 11 : actRec_(),
228 11 : resRec_(),
229 11 : simdata_p(false),
230 33 : ssvp_p()
231 : {
232 : // cout << "This is the NEW VI2-aware Calibrater" << endl;
233 11 : }
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 33 : Calibrater::~Calibrater()
332 : {
333 11 : cleanup();
334 11 : if (msmc_p) delete msmc_p; msmc_p=0;
335 11 : if (ms_p) delete ms_p; ms_p=0;
336 11 : if (hist_p) delete hist_p; hist_p=0;
337 :
338 22 : }
339 :
340 11 : Calibrater* Calibrater::factory(Bool old) {
341 :
342 11 : Calibrater* cal(NULL);
343 :
344 11 : if (old)
345 0 : cal = new OldCalibrater();
346 : else
347 : //throw(AipsError("VI2-aware Calibrater not yet available..."));
348 11 : cal = new Calibrater();
349 :
350 11 : 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 0 : 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 0 : Bool Calibrater::initialize(MeasurementSet& inputMS,
381 : Bool compress,
382 : Bool addScratch, Bool addModel) {
383 :
384 0 : logSink() << LogOrigin("Calibrater","") << LogIO::NORMAL3;
385 :
386 : try {
387 0 : timer_p.mark();
388 :
389 : // Set pointer ms_p from input MeasurementSet
390 0 : if (ms_p) {
391 0 : *ms_p=inputMS;
392 : } else {
393 0 : ms_p = new MeasurementSet(inputMS);
394 0 : AlwaysAssert(ms_p,AipsError);
395 : };
396 :
397 : // Remember the ms's name
398 0 : msname_p=ms_p->tableName();
399 :
400 : // Initialize the meta-info server that will be shared with VisCals
401 0 : if (msmc_p) delete msmc_p;
402 0 : msmc_p = new MSMetaInfoForCal(*ms_p);
403 :
404 : // Add/init scr cols, if requested (init is hard-wired)
405 0 : 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 0 : << LogIO::POST;
415 :
416 0 : if (mssel_p) delete mssel_p;
417 0 : mssel_p=new MeasurementSet(*ms_p);
418 :
419 : // Create the associated VisEquation
420 : // TBD: move to ctor and make it non-pointer
421 0 : if (ve_p) {
422 0 : delete ve_p;
423 0 : ve_p=0;
424 : };
425 0 : ve_p=new VisEquation();
426 :
427 : // Reset the apply/solve VisCals
428 0 : reset(true,true);
429 :
430 0 : 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 0 : 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 0 : 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 0 : logSink() << "Selecting data" << LogIO::POST;
500 :
501 : // Apply selection to the original MeasurementSet
502 0 : logSink() << "Performing selection on MeasurementSet" << endl;
503 :
504 0 : if (mssel_p) {
505 0 : delete mssel_p;
506 0 : mssel_p=0;
507 : };
508 :
509 : // Report non-trivial user selections
510 0 : if (time!="")
511 0 : logSink() << " Selecting on time: '" << time << "'" << endl;
512 0 : if (spw!="")
513 0 : logSink() << " Selecting on spw: '" << spw << "'" << endl;
514 0 : if (scan!="")
515 0 : logSink() << " Selecting on scan: '" << scan << "'" << endl;
516 0 : if (field!="")
517 0 : logSink() << " Selecting on field: '" << field << "'" << endl;
518 0 : if (intent!="")
519 0 : logSink() << " Selecting on intent: '" << intent << "'" << endl;
520 0 : if(obsIDs != "")
521 0 : logSink() << " Selecting by observation IDs: '" << obsIDs << "'" << endl;
522 0 : if (baseline!="")
523 0 : logSink() << " Selecting on antenna/baseline: '" << baseline << "'" << endl;
524 0 : if (uvrange!="")
525 0 : logSink() << " Selecting on uvrange: '" << uvrange << "'" << endl;
526 0 : if (msSelect!="")
527 0 : logSink() << " Selecting with TaQL: '" << msSelect << "'" << endl;
528 0 : logSink() << LogIO::POST;
529 :
530 :
531 : // Assume no selection, for starters
532 0 : mssel_p = new MeasurementSet(*ms_p);
533 :
534 : // Apply user-supplied selection
535 0 : Bool nontrivsel=false;
536 :
537 : // Ensure use of a fresh MSSelection object
538 0 : if (mss_p) { delete mss_p; mss_p=NULL; }
539 0 : mss_p=new MSSelection();
540 0 : nontrivsel= mssSetData(*ms_p,
541 0 : *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 0 : mssel_p->setMemoryResidentSubtables(ms_p->getMrsEligibility());
549 :
550 : // If non-trivial MSSelection invoked and nrow reduced:
551 0 : if(nontrivsel && mssel_p->nrow()<ms_p->nrow()) {
552 :
553 : // Escape if no rows selected
554 0 : if (mssel_p->nrow()==0)
555 0 : throw(AipsError("Specified selection selects zero rows!"));
556 :
557 : // ...otherwise report how many rows are selected
558 0 : logSink() << "By selection " << ms_p->nrow()
559 0 : << " rows are reduced to " << mssel_p->nrow()
560 0 : << LogIO::POST;
561 :
562 : // Revise the msmc_p so it contains only the _selected_ MS
563 0 : if (msmc_p) delete msmc_p;
564 0 : msmc_p = new MSMetaInfoForCal(*mssel_p, ms_p->tableName());
565 :
566 : }
567 : else {
568 : // Selection did nothing:
569 0 : 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 0 : if (chanmode=="none") {
575 0 : 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 0 : };
607 :
608 :
609 0 : 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 0 : logSink() << LogOrigin("Calibrater",
621 : "setapply(type, t, table, spw, field, interp, calwt, spwmap, opacity)")
622 0 : << LogIO::NORMAL;
623 :
624 :
625 : // cout << "Calibrater::setapply: field="<< field << endl;
626 :
627 : // Set record format for calibration table application information
628 0 : RecordDesc applyparDesc;
629 0 : applyparDesc.addField ("t", TpDouble);
630 0 : applyparDesc.addField ("table", TpString);
631 0 : applyparDesc.addField ("interp", TpString);
632 0 : applyparDesc.addField ("spw", TpArrayInt);
633 : // applyparDesc.addField ("field", TpArrayInt);
634 0 : applyparDesc.addField ("fieldstr", TpString);
635 0 : applyparDesc.addField ("calwt",TpBool);
636 0 : applyparDesc.addField ("spwmap",TpArrayInt);
637 0 : applyparDesc.addField ("opacity",TpArrayDouble);
638 :
639 : // Create record with the requisite field values
640 0 : Record applypar(applyparDesc);
641 0 : applypar.define ("t", t);
642 0 : 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 0 : applypar.define ("interp", interp);
651 0 : applypar.define ("spw",getSpwIdx(spw));
652 : // applypar.define ("field",getFieldIdx(field));
653 0 : applypar.define ("fieldstr",field);
654 0 : applypar.define ("calwt",calwt);
655 0 : applypar.define ("spwmap",spwmap);
656 0 : applypar.define ("opacity", opacity);
657 :
658 0 : String upType=type;
659 0 : upType.upcase();
660 0 : if (upType=="")
661 : // Get type from table
662 0 : upType = calTableType(table);
663 :
664 0 : return setapply(upType,applypar);
665 :
666 0 : }
667 :
668 0 : Bool Calibrater::setapply (const String& type,
669 : const Record& applypar)
670 : {
671 0 : logSink() << LogOrigin("Calibrater", "setapply(type, applypar)");
672 :
673 : // First try to create the requested VisCal object
674 0 : VisCal *vc(NULL);
675 :
676 : try {
677 :
678 0 : if(!ok())
679 0 : throw(AipsError("Calibrater not prepared for setapply."));
680 :
681 0 : String upType=type;
682 0 : upType.upcase();
683 :
684 : logSink() << LogIO::NORMAL
685 : << "Arranging to APPLY:"
686 0 : << LogIO::POST;
687 :
688 : // Add a new VisCal to the apply list
689 0 : vc = createVisCal(upType,*msmc_p);
690 :
691 0 : vc->setApply(applypar);
692 :
693 : logSink() << LogIO::NORMAL << ". "
694 0 : << vc->applyinfo()
695 0 : << LogIO::POST;
696 :
697 0 : } 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 0 : uInt napp=vc_p.nelements();
711 0 : vc_p.resize(napp+1,false,true);
712 0 : vc_p[napp] = vc;
713 0 : vc=NULL;
714 :
715 : // Maintain sort of apply list
716 0 : ve_p->setapply(vc_p);
717 :
718 0 : 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 0 : Bool Calibrater::setcallib2(Record callib, const MeasurementSet* ms) {
766 :
767 0 : logSink() << LogOrigin("Calibrater", "setcallib2(callib)");
768 :
769 : // cout << "Calibrater::setcallib2(callib) : " << boolalpha << callib << endl;
770 :
771 0 : uInt ntab=callib.nfields();
772 :
773 : // cout << "callib.nfields() = " << ntab << endl;
774 :
775 : // Do some preliminary per-table verification
776 0 : for (uInt itab=0;itab<ntab;++itab) {
777 :
778 0 : 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 0 : if (!Table::isReadable(tabname))
787 0 : throw(AipsError("Caltable "+tabname+" does not exist."));
788 :
789 0 : }
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 0 : const MeasurementSet *lmsp(0);
801 0 : 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 0 : lmsp=mssel_p;
811 : }
812 : // Reference for use below
813 0 : const MeasurementSet &lms(*lmsp);
814 :
815 :
816 0 : for (uInt itab=0;itab<ntab;++itab) {
817 :
818 0 : String tabname=callib.name(itab);
819 :
820 : // Get the type from the table
821 0 : String upType=calTableType(tabname);
822 0 : upType.upcase();
823 :
824 : // Add table name to the record
825 0 : Record thistabrec=callib.asrwRecord(itab);
826 0 : thistabrec.define("tablename",tabname);
827 :
828 : // First try to create the requested VisCal object
829 0 : 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 0 : << LogIO::POST;
839 :
840 : // Add a new VisCal to the apply list
841 0 : vc = createVisCal(upType,*msmc_p);
842 :
843 : // ingest this table according to its callib
844 0 : 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 0 : uInt napp=vc_p.nelements();
860 0 : vc_p.resize(napp+1,false,true);
861 0 : vc_p[napp] = vc;
862 0 : vc=NULL;
863 :
864 : // Maintain sort of apply list
865 0 : 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 0 : }
875 :
876 : // Signal use of CalLibrary
877 0 : usingCalLibrary_=true;
878 :
879 : // All ok, if we get this far!
880 0 : 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 0 : 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 0 : logSink() << LogOrigin("Calibrater","setsolve") << LogIO::NORMAL3;
941 :
942 : // Create a record description containing the solver parameters
943 0 : RecordDesc solveparDesc;
944 0 : solveparDesc.addField ("solint", TpString);
945 0 : solveparDesc.addField ("preavg", TpDouble);
946 0 : solveparDesc.addField ("apmode", TpString);
947 0 : solveparDesc.addField ("refant", TpArrayInt);
948 0 : solveparDesc.addField ("refantmode", TpString);
949 0 : solveparDesc.addField ("minblperant", TpInt);
950 0 : solveparDesc.addField ("table", TpString);
951 0 : solveparDesc.addField ("append", TpBool);
952 0 : solveparDesc.addField ("solnorm", TpBool);
953 0 : solveparDesc.addField ("normtype", TpString);
954 0 : solveparDesc.addField ("type", TpString);
955 0 : solveparDesc.addField ("combine", TpString);
956 0 : solveparDesc.addField ("maxgap", TpInt);
957 0 : solveparDesc.addField ("cfcache", TpString);
958 0 : solveparDesc.addField ("painc", TpDouble);
959 0 : solveparDesc.addField ("fitorder", TpInt);
960 0 : solveparDesc.addField ("solmode", TpString);
961 0 : solveparDesc.addField ("rmsthresh", TpArrayDouble);
962 :
963 : // fringe-fit specific fields
964 0 : solveparDesc.addField ("zerorates", TpBool);
965 0 : solveparDesc.addField ("minsnr", TpFloat);
966 0 : solveparDesc.addField ("globalsolve", TpBool);
967 0 : solveparDesc.addField ("delaywindow", TpArrayDouble);
968 0 : solveparDesc.addField ("ratewindow", TpArrayDouble);
969 0 : solveparDesc.addField ("niter", TpInt);
970 0 : solveparDesc.addField ("corrcomb", TpString);
971 0 : solveparDesc.addField ("paramactive", TpArrayBool);
972 0 : solveparDesc.addField ("concatspws", TpBool);
973 :
974 : // single dish specific fields
975 0 : solveparDesc.addField ("fraction", TpFloat);
976 0 : solveparDesc.addField ("numedge", TpInt);
977 0 : solveparDesc.addField ("radius", TpString);
978 0 : solveparDesc.addField ("smooth", TpBool);
979 :
980 :
981 : // Create a solver record with the requisite field values
982 0 : Record solvepar(solveparDesc);
983 0 : solvepar.define ("solint", solint);
984 0 : solvepar.define ("preavg", preavg);
985 0 : String upmode=apmode;
986 0 : upmode.upcase();
987 0 : solvepar.define ("apmode", upmode);
988 0 : solvepar.define ("refant", getRefantIdxList(refant));
989 0 : solvepar.define ("refantmode", refantmode);
990 0 : solvepar.define ("minblperant", minblperant);
991 0 : solvepar.define ("table", table);
992 0 : solvepar.define ("append", append);
993 0 : solvepar.define ("solnorm", solnorm);
994 0 : solvepar.define ("normtype", normtype);
995 : // Fringe-fit specific
996 0 : solvepar.define ("minsnr", minsnr);
997 0 : solvepar.define ("zerorates", zerorates);
998 0 : solvepar.define ("globalsolve", globalsolve);
999 0 : solvepar.define ("niter", niter);
1000 0 : solvepar.define ("corrcomb", corrcomb);
1001 0 : solvepar.define ("delaywindow", delaywindow);
1002 0 : solvepar.define ("ratewindow", ratewindow);
1003 0 : solvepar.define ("solmode", solmode);
1004 0 : solvepar.define ("rmsthresh", rmsthresh);
1005 0 : solvepar.define ("paramactive", paramactive);
1006 0 : solvepar.define ("concatspws", concatspws);
1007 :
1008 0 : String uptype=type;
1009 0 : uptype.upcase();
1010 0 : solvepar.define ("type", uptype);
1011 :
1012 0 : String upcomb=combine;
1013 0 : upcomb.upcase();
1014 0 : solvepar.define("combine",upcomb);
1015 0 : solvepar.define("maxgap",fillgaps);
1016 0 : solvepar.define ("cfcache", cfcache);
1017 0 : solvepar.define ("painc", painc);
1018 0 : solvepar.define("fitorder", fitorder);
1019 :
1020 : // single dish specific
1021 0 : solvepar.define("fraction", fraction);
1022 0 : solvepar.define("numedge", numedge);
1023 0 : solvepar.define("radius", radius);
1024 0 : solvepar.define("smooth", smooth);
1025 :
1026 0 : return setsolve(type,solvepar);
1027 0 : }
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 0 : Bool Calibrater::setsolve (const String& type,
1125 : const Record& solvepar) {
1126 :
1127 : // Attempt to create the solvable object
1128 0 : SolvableVisCal *svc(NULL);
1129 : try {
1130 :
1131 0 : if(!ok())
1132 0 : throw(AipsError("Calibrater not prepared for setsolve."));
1133 :
1134 0 : String upType = type;
1135 0 : upType.upcase();
1136 :
1137 : // Clean out any old solve that was lying around
1138 0 : unsetsolve();
1139 :
1140 : logSink() << LogIO::NORMAL
1141 : << "Arranging to SOLVE:"
1142 0 : << LogIO::POST;
1143 :
1144 : // Create the new SolvableVisCal
1145 0 : svc = createSolvableVisCal(upType,*msmc_p);
1146 0 : svc->setSolve(solvepar);
1147 :
1148 : logSink() << LogIO::NORMAL << ". "
1149 0 : << svc->solveinfo()
1150 0 : << LogIO::POST;
1151 :
1152 : // Creation apparently successful, keep it
1153 0 : svc_p=svc;
1154 0 : svc=NULL;
1155 :
1156 : // if calibration specific data filter is necessary
1157 : // keep configuration parameter as a record
1158 0 : setCalFilterConfiguration(upType, solvepar);
1159 :
1160 0 : return true;
1161 :
1162 0 : } 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 0 : 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 0 : << LogIO::POST;
1338 :
1339 0 : Int napp(vc_p.nelements());
1340 0 : if (napp>0)
1341 0 : for (Int iapp=0;iapp<napp;++iapp)
1342 : logSink() << LogIO::NORMAL << ". "
1343 0 : << vc_p[iapp]->applyinfo()
1344 0 : << LogIO::POST;
1345 : else
1346 : logSink() << LogIO::NORMAL << ". "
1347 : << "(None)"
1348 0 : << LogIO::POST;
1349 :
1350 0 : return true;
1351 :
1352 : }
1353 :
1354 :
1355 0 : 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 0 : << LogIO::POST;
1362 :
1363 0 : if (svc_p)
1364 : logSink() << LogIO::NORMAL << ". "
1365 0 : << svc_p->solveinfo()
1366 0 : << LogIO::POST;
1367 : else
1368 : logSink() << LogIO::NORMAL << ". "
1369 : << "(None)"
1370 0 : << LogIO::POST;
1371 :
1372 0 : return true;
1373 : }
1374 :
1375 11 : 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 11 : if (apply)
1381 11 : unsetapply();
1382 :
1383 : // Delete the VisCal solve object
1384 11 : if (solve)
1385 11 : unsetsolve();
1386 :
1387 11 : return true;
1388 : }
1389 :
1390 : // Delete all (default) or one VisCal in apply list
1391 11 : Bool Calibrater::unsetapply(const Int& which) {
1392 :
1393 : // logSink() << LogOrigin("Calibrater","unsetapply") << LogIO::NORMAL;
1394 :
1395 : try {
1396 11 : if (which<0) {
1397 11 : for (uInt i=0;i<vc_p.nelements();i++)
1398 0 : if (vc_p[i]) delete vc_p[i];
1399 11 : 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 11 : if(ve_p) ve_p->setapply(vc_p);
1407 :
1408 11 : 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 11 : Bool Calibrater::unsetsolve() {
1421 :
1422 : // logSink() << LogOrigin("Calibrater","unsetsolve") << LogIO::NORMAL;
1423 :
1424 : try {
1425 11 : if (svc_p) delete svc_p;
1426 11 : svc_p=NULL;
1427 :
1428 11 : if(ve_p) ve_p->setsolve(*svc_p);
1429 :
1430 11 : 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 0 : Calibrater::setCorrDepFlags(const Bool& corrDepFlags)
1443 : {
1444 :
1445 0 : logSink() << LogOrigin("Calibrater","setCorrDepFlags") << LogIO::NORMAL;
1446 :
1447 : // Set it
1448 0 : corrDepFlags_=corrDepFlags;
1449 :
1450 0 : logSink() << "Setting correlation dependent flags = " << (corrDepFlags_ ? "True" : "False") << LogIO::POST;
1451 :
1452 0 : 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 0 : Bool Calibrater::solve() {
2606 :
2607 0 : logSink() << LogOrigin("Calibrater","solve") << LogIO::NORMAL3;
2608 :
2609 : try {
2610 :
2611 0 : if (!ok())
2612 0 : throw(AipsError("Calibrater not prepared for solve."));
2613 :
2614 : // Handle nothing-to-solve-for case
2615 0 : if (!svc_p)
2616 0 : throw(AipsError("Please run setsolve before attempting to solve."));
2617 :
2618 : // Handle specified caltable
2619 0 : if (svc_p) {
2620 :
2621 : // If table exists (readable) and not deletable
2622 : // we have to abort (append=T requires deletable)
2623 0 : if ( Table::isReadable(svc_p->calTableName()) &&
2624 0 : !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 0 : ve_p->setsolve(*svc_p);
2631 :
2632 : // Ensure apply list properly sorted w.r.t. solvable term
2633 0 : ve_p->setapply(vc_p);
2634 :
2635 : // Report what is being applied and solved-for
2636 0 : applystate();
2637 0 : solvestate();
2638 :
2639 : // Report correct/corrupt apply order
2640 : // ve_p->state();
2641 :
2642 : // Generally use standard solver
2643 0 : if (svc_p->useGenericGatherForSolve())
2644 0 : 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 0 : } catch (AipsError x) {
2653 0 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
2654 :
2655 0 : logSink() << "Reseting entire solve/apply state." << LogIO::POST;
2656 0 : reset();
2657 :
2658 0 : throw(AipsError("Error in Calibrater::solve."));
2659 : return false;
2660 0 : }
2661 :
2662 0 : return true;
2663 :
2664 : }
2665 :
2666 :
2667 0 : Bool Calibrater::summarize_uncalspws(const Vector<Bool>& uncalspw,
2668 : const String& origin,
2669 : Bool strictflag)
2670 : {
2671 0 : Bool hadprob = false;
2672 0 : uInt totNspw = uncalspw.nelements();
2673 :
2674 0 : for(uInt i = 0; i < totNspw; ++i){
2675 0 : if(uncalspw[i]){
2676 0 : hadprob = true;
2677 0 : break;
2678 : }
2679 : }
2680 0 : 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 0 : return !hadprob;
2698 : }
2699 :
2700 0 : 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 0 : logSink() << LogOrigin("Calibrater","fluxscale") << LogIO::NORMAL3;
2721 :
2722 : //outfile check
2723 0 : if (outfile=="")
2724 0 : throw(AipsError("output fluxscaled caltable name must be specified!"));
2725 : else {
2726 0 : 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 0 : Vector<Int> refidx(0);
2732 :
2733 0 : if (refFields.length()>0)
2734 0 : refidx=getFieldIdx(refFields);
2735 : else
2736 0 : throw(AipsError("A reference field must be specified!"));
2737 :
2738 0 : if (tranFields.length()>0)
2739 0 : tranidx=getFieldIdx(tranFields);
2740 :
2741 : // Call Vector<Int> version:
2742 0 : fluxscale(infile,outfile,refidx,refSpwMap,tranidx,append,inGainThres,antSel,timerangeSel,
2743 : scanSel,oFluxScaleFactor,oListFile,incremental,fitorder,display);
2744 :
2745 0 : }
2746 :
2747 0 : 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 0 : logSink() << LogOrigin("Calibrater","fluxscale") << LogIO::NORMAL3;
2767 :
2768 0 : SolvableVisCal *fsvj_(NULL);
2769 : try {
2770 : // If infile is Calibration table
2771 0 : if (Table::isReadable(infile) &&
2772 0 : TableUtil::tableInfo(infile).type()=="Calibration") {
2773 :
2774 : // get calibration type
2775 0 : String caltype;
2776 0 : caltype = TableUtil::tableInfo(infile).subType();
2777 : logSink() << "Table " << infile
2778 : << " is of type: "<< caltype
2779 0 : << LogIO::POST;
2780 0 : String message="Table "+infile+" is of type: "+caltype;
2781 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
2782 :
2783 : // form selection
2784 0 : String select="";
2785 : // Selection is empty for case of no tran specification
2786 0 : if (tranField.nelements()>0) {
2787 :
2788 : // All selected fields
2789 0 : Vector<Int> allflds = concatenateArray(refField,tranField);
2790 :
2791 : // Assemble TaQL
2792 0 : ostringstream selectstr;
2793 0 : selectstr << "FIELD_ID IN [";
2794 0 : for (Int iFld=0; iFld<allflds.shape(); iFld++) {
2795 0 : if (iFld>0) selectstr << ", ";
2796 0 : selectstr << allflds(iFld);
2797 : }
2798 0 : selectstr << "]";
2799 0 : select=selectstr.str();
2800 0 : }
2801 :
2802 : // Construct proper SVC object
2803 0 : if (caltype == "G Jones") {
2804 0 : 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 0 : RecordDesc applyparDesc;
2818 0 : applyparDesc.addField ("table", TpString);
2819 0 : applyparDesc.addField ("select", TpString);
2820 0 : Record applypar(applyparDesc);
2821 0 : applypar.define ("table", infile);
2822 0 : applypar.define ("select", select);
2823 0 : fsvj_->setApply(applypar);
2824 :
2825 : //Bool incremental=false;
2826 : // Make fluxscale calculation
2827 0 : Vector<String> fldnames(MSFieldColumns(ms_p->field()).name().getColumn());
2828 : //fsvj_->fluxscale(refField,tranField,refSpwMap,fldnames,oFluxScaleFactor,
2829 0 : 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 0 : String out(outfile);
2835 0 : if (out.length()==0)
2836 0 : out = infile;
2837 :
2838 : // Store result
2839 0 : if (append) {
2840 0 : logSink() << "Appending result to " << out << LogIO::POST;
2841 0 : String message="Appending result to "+out;
2842 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
2843 0 : } else {
2844 0 : logSink() << "Storing result in " << out << LogIO::POST;
2845 0 : String message="Storing result in "+out;
2846 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
2847 0 : }
2848 0 : fsvj_->storeNCT(out,append);
2849 :
2850 : // Clean up
2851 0 : delete fsvj_;
2852 :
2853 0 : } 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 0 : 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 : {
3000 :
3001 : // TBD: support append?
3002 : // TBD: spw selection?
3003 :
3004 0 : logSink() << LogOrigin("Calibrater","smooth") << LogIO::NORMAL;
3005 :
3006 0 : logSink() << "Beginning smoothing/interpolating method." << LogIO::POST;
3007 :
3008 :
3009 : // A pointer to an SVC
3010 0 : SolvableVisCal *svc(NULL);
3011 :
3012 : try {
3013 :
3014 : // Handle no in file
3015 0 : if (infile=="")
3016 0 : throw(AipsError("Please specify an input calibration table."));
3017 :
3018 : // Handle bad smoothtype
3019 0 : if (smoothtype!="mean" && smoothtype!="median")
3020 0 : throw(AipsError("Unrecognized smooth type!"));
3021 :
3022 : // Handle bad smoothtime
3023 0 : if (smoothtime<=0)
3024 0 : throw(AipsError("Please specify a strictly positive smoothtime."));
3025 :
3026 : // Handle no outfile
3027 0 : if (outfile=="") {
3028 0 : outfile=infile;
3029 : logSink() << "Will overwrite input file with smoothing result."
3030 0 : << LogIO::POST;
3031 : }
3032 :
3033 :
3034 0 : svc = createSolvableVisCal(calTableType(infile),*msmc_p);
3035 :
3036 0 : if (svc->smoothable()) {
3037 :
3038 : // Fill calibration table using setApply
3039 0 : RecordDesc applyparDesc;
3040 0 : applyparDesc.addField ("table", TpString);
3041 0 : Record applypar(applyparDesc);
3042 0 : applypar.define ("table", infile);
3043 0 : svc->setApply(applypar);
3044 :
3045 : // Convert refFields/transFields to index lists
3046 0 : Vector<Int> fldidx(0);
3047 0 : if (fields.length()>0)
3048 0 : fldidx=getFieldIdx(fields);
3049 :
3050 : // Delegate to SVC
3051 0 : svc->smooth(fldidx,smoothtype,smoothtime);
3052 :
3053 : // Store the result on disk
3054 : // if (append) logSink() << "Appending result to " << outfile << LogIO::POST;
3055 : //else
3056 0 : logSink() << "Storing result in " << outfile << LogIO::POST;
3057 :
3058 :
3059 0 : if (outfile != "")
3060 0 : svc->calTableName()=outfile;
3061 0 : svc->storeNCT();
3062 :
3063 : // Clean up
3064 0 : if (svc) delete svc; svc=NULL;
3065 :
3066 : // Apparently, it worked
3067 0 : return true;
3068 :
3069 0 : }
3070 : else
3071 0 : throw(AipsError("This type ("+svc->typeName()+") does not support smoothing."));
3072 :
3073 0 : } catch (AipsError x) {
3074 :
3075 : logSink() << LogIO::SEVERE
3076 : << "Caught Exception: "
3077 : << x.getMesg()
3078 0 : << LogIO::POST;
3079 : // Clean up
3080 0 : if (svc) delete svc; svc=NULL;
3081 :
3082 0 : throw(AipsError("Error in Calibrater::smooth."));
3083 :
3084 : return false;
3085 0 : }
3086 : return false;
3087 : }
3088 :
3089 : // Apply new reference antenna to calibration
3090 0 : Bool Calibrater::reRefant(const casacore::String& infile,
3091 : casacore::String& outfile,
3092 : const casacore::String& refantmode,
3093 : const casacore::String& refant)
3094 : {
3095 :
3096 0 : logSink() << LogOrigin("Calibrater","reRefant") << LogIO::NORMAL;
3097 :
3098 : // logSink() << "Beginning smoothing/interpolating method." << LogIO::POST;
3099 :
3100 :
3101 : // A pointer to an SVC
3102 0 : SolvableVisJones *svj(NULL);
3103 :
3104 : try {
3105 :
3106 : // Handle no in file
3107 0 : if (infile=="")
3108 0 : throw(AipsError("Please specify an input calibration table."));
3109 :
3110 : // Handle bad refantmode
3111 0 : if (refantmode!="strict" &&
3112 0 : refantmode!="flex")
3113 0 : throw(AipsError("Unrecognized refantmode!"));
3114 :
3115 : // Handle no outfile
3116 0 : if (outfile=="") {
3117 0 : outfile=infile;
3118 : logSink() << "Will overwrite input file with rerefant result."
3119 0 : << LogIO::POST;
3120 : }
3121 :
3122 :
3123 0 : svj = (SolvableVisJones*) createSolvableVisCal(calTableType(infile),*msmc_p);
3124 :
3125 : // Fill calibration table using setApply
3126 0 : RecordDesc applyparDesc;
3127 0 : applyparDesc.addField ("table", TpString);
3128 0 : Record applypar(applyparDesc);
3129 0 : applypar.define ("table", infile);
3130 0 : svj->setApply(applypar);
3131 :
3132 : // Do the work
3133 0 : svj->refantmode() = refantmode;
3134 0 : svj->refantlist().reference(getRefantIdxList(refant)); // replaces the default list
3135 0 : svj->applyRefAnt();
3136 :
3137 : // Store the result on disk
3138 0 : logSink() << "Storing result in " << outfile << LogIO::POST;
3139 :
3140 0 : if (outfile != "")
3141 0 : svj->calTableName()=outfile;
3142 0 : svj->storeNCT();
3143 :
3144 : // Clean up
3145 0 : if (svj) delete svj; svj=NULL;
3146 :
3147 : // Apparently, it worked
3148 0 : return true;
3149 :
3150 0 : } catch (AipsError x) {
3151 :
3152 : logSink() << LogIO::SEVERE
3153 : << "Caught Exception: "
3154 : << x.getMesg()
3155 0 : << LogIO::POST;
3156 : // Clean up
3157 0 : if (svj) delete svj; svj=NULL;
3158 :
3159 0 : throw(AipsError("Error in Calibrater::reRefant."));
3160 :
3161 : return false;
3162 0 : }
3163 : return false;
3164 : }
3165 :
3166 :
3167 : // List a calibration table
3168 0 : Bool Calibrater::listCal(const String& infile,
3169 : const String& field,
3170 : const String& antenna,
3171 : const String& spw,
3172 : const String& listfile,
3173 : const Int& pagerows) {
3174 :
3175 0 : SolvableVisCal *svc(NULL);
3176 0 : logSink() << LogOrigin("Calibrater","listCal");
3177 :
3178 : try {
3179 :
3180 : // Trap (currently) unsupported types
3181 0 : if (upcase(calTableType(infile))=="GSPLINE" ||
3182 0 : upcase(calTableType(infile))=="BPOLY")
3183 0 : throw(AipsError("GSPLINE and BPOLY tables cannot currently be listed."));
3184 :
3185 : // Get user's selected fields, ants
3186 0 : Vector<Int> ufldids=getFieldIdx(field);
3187 0 : Vector<Int> uantids=getAntIdx(antenna);
3188 :
3189 0 : String newSpw = spw;
3190 0 : Bool defaultSelect = false;
3191 0 : if (spw.empty()) { // list all channels (default)
3192 0 : defaultSelect = true;
3193 0 : newSpw = "*";
3194 : logSink() << LogIO::NORMAL1 << "Spws selected: ALL" << endl
3195 0 : << "Channels selected: ALL" << LogIO::POST;
3196 : }
3197 : // Get user's selected spw and channels
3198 0 : Vector<Int> uspwids=getSpwIdx(newSpw);
3199 0 : Matrix<Int> uchanids=getChanIdx(newSpw);
3200 0 : if (!defaultSelect) {
3201 : logSink() << LogIO::NORMAL1 << "Spw and Channel selection matrix: "
3202 : << endl << "Each rows shows: [ Spw , Start Chan , Stop Chan , Chan Step ]"
3203 0 : << endl << uchanids << LogIO::POST;
3204 : }
3205 : logSink() << LogIO::DEBUG2
3206 : << "uspwids = " << uspwids << endl
3207 0 : << "uchanids = " << uchanids << LogIO::POST;
3208 :
3209 : // By default, do first spw, first chan
3210 0 : if (uspwids.nelements()==0) {
3211 0 : uchanids.resize(1,4);
3212 0 : uchanids=0;
3213 : }
3214 :
3215 : // Set record format for calibration table application information
3216 0 : RecordDesc applyparDesc;
3217 0 : applyparDesc.addField ("table", TpString);
3218 :
3219 : // Create record with the requisite field values
3220 0 : Record applypar(applyparDesc);
3221 0 : applypar.define ("table", infile);
3222 :
3223 : // Generate the VisCal to be listed
3224 0 : svc = createSolvableVisCal(calTableType(infile),*msmc_p);
3225 0 : svc->setApply(applypar);
3226 :
3227 : // list it
3228 0 : svc->listCal(ufldids,uantids,uchanids, //uchanids(0,0),uchanids(0,1),
3229 : listfile,pagerows);
3230 :
3231 0 : if (svc) delete svc; svc=NULL;
3232 :
3233 0 : return true;
3234 :
3235 0 : } catch (AipsError x) {
3236 :
3237 : logSink() << LogIO::SEVERE
3238 : << "Caught Exception: "
3239 : << x.getMesg()
3240 0 : << LogIO::POST;
3241 : // Clean up
3242 0 : if (svc) delete svc; svc=NULL;
3243 :
3244 0 : throw(AipsError("Error in Calibrater::listCal."));
3245 :
3246 : return false;
3247 0 : }
3248 : return false;
3249 :
3250 : }
3251 :
3252 :
3253 :
3254 0 : Bool Calibrater::updateCalTable(const String& caltable) {
3255 :
3256 : // Call the SVC method that knows how
3257 0 : return NewCalTable::CTBackCompat(caltable);
3258 :
3259 : }
3260 :
3261 0 : void Calibrater::selectChannel(const String& spw) {
3262 :
3263 0 : if (mss_p && mssel_p) {
3264 :
3265 : // Refresh the frequencySelections object to feed to VI2, if relevant
3266 0 : frequencySelections_p.reset(new vi::FrequencySelections());
3267 :
3268 0 : vi::FrequencySelectionUsingChannels usingChannels;
3269 0 : usingChannels.add(*mss_p,mssel_p);
3270 0 : frequencySelections_p->add(usingChannels);
3271 :
3272 : // cout << usingChannels.toString() << endl;
3273 : // cout << "FS.size() = " << frequencySelections_p->size() << endl;
3274 :
3275 0 : }
3276 :
3277 : // TBD: Does frequencySelections_p support string info for reporting to logger?
3278 :
3279 0 : Matrix<Int> chansel = getChanIdx(spw);
3280 0 : uInt nselspw=chansel.nrow();
3281 :
3282 0 : if (nselspw==0)
3283 : logSink() << "Frequency selection: Selecting all channels in all spws."
3284 0 : << LogIO::POST;
3285 : else {
3286 :
3287 0 : logSink() << "Frequency selection: " << LogIO::POST;
3288 :
3289 : // Trap non-unit step (for now)
3290 0 : if (ntrue(chansel.column(3)==1)!=nselspw) {
3291 : logSink() << LogIO::WARN
3292 : << "Calibration does not support non-unit channel stepping; "
3293 : << "using step=1."
3294 0 : << LogIO::POST;
3295 0 : chansel.column(3)=1;
3296 : }
3297 :
3298 0 : Vector<Int> uspw(chansel.column(0));
3299 0 : Vector<Int> ustart(chansel.column(1));
3300 0 : Vector<Int> uend(chansel.column(2));
3301 0 : Vector<Int> ustep(chansel.column(3));
3302 :
3303 0 : logSink() << LogIO::NORMAL;
3304 0 : for (uInt i=0;i<nselspw;++i) {
3305 :
3306 0 : logSink() << ". Spw " << uspw(i) << ":"
3307 0 : << ustart(i) << "~" << uend(i)
3308 0 : << " (" << uend(i)-ustart(i)+1 << " channels,"
3309 0 : << " step by " << ustep(i) << ")"
3310 0 : << endl;
3311 :
3312 : } // i
3313 0 : logSink() << LogIO::POST;
3314 :
3315 0 : } // non-triv spw selection
3316 :
3317 0 : }
3318 :
3319 :
3320 11 : Bool Calibrater::cleanup() {
3321 :
3322 : // logSink() << LogOrigin("Calibrater","cleanup") << LogIO::NORMAL;
3323 :
3324 : // Delete the VisCals
3325 11 : reset();
3326 :
3327 : // Delete derived dataset stuff
3328 11 : if(mssel_p) delete mssel_p; mssel_p=0;
3329 11 : if(mss_p) delete mss_p; mss_p=0;
3330 11 : frequencySelections_p.reset();
3331 :
3332 : // Delete the current VisEquation
3333 11 : if(ve_p) delete ve_p; ve_p=0;
3334 :
3335 11 : return true;
3336 :
3337 : }
3338 :
3339 : // Parse refant specification
3340 0 : Vector<Int> Calibrater::getRefantIdxList(const String& refant) {
3341 :
3342 0 : Vector<Int> irefant;
3343 0 : if (refant.length()==0) {
3344 : // Nothing specified, return -1 in single-element vector
3345 0 : irefant.resize(1);
3346 0 : irefant(0)=-1;
3347 : }
3348 : else {
3349 : // parse the specification
3350 0 : MSSelection msselect;
3351 0 : msselect.setAntennaExpr(refant);
3352 0 : Vector<Int> iref=msselect.getAntenna1List(mssel_p);
3353 0 : if (anyLT(iref,0))
3354 0 : cout << "Negated selection (via '!') not yet implemented for refant," << endl << " and will be ignored." << endl;
3355 0 : irefant=iref(iref>-1).getCompressedArray();
3356 0 : if (irefant.nelements()==0) {
3357 0 : irefant.resize(1);
3358 0 : irefant(0)=-1;
3359 : }
3360 0 : }
3361 :
3362 0 : return irefant;
3363 0 : }
3364 :
3365 :
3366 : // Interpret refant *index*
3367 0 : Vector<Int> Calibrater::getAntIdx(const String& antenna) {
3368 :
3369 0 : MSSelection msselect;
3370 0 : msselect.setAntennaExpr(antenna);
3371 0 : return msselect.getAntenna1List(mssel_p);
3372 :
3373 0 : }
3374 :
3375 : // Interpret field indices (MSSelection)
3376 0 : Vector<Int> Calibrater::getFieldIdx(const String& fields) {
3377 :
3378 0 : MSSelection mssel;
3379 0 : mssel.setFieldExpr(fields);
3380 0 : return mssel.getFieldList(mssel_p);
3381 :
3382 0 : }
3383 :
3384 : // Interpret spw indices (MSSelection)
3385 0 : Vector<Int> Calibrater::getSpwIdx(const String& spws) {
3386 :
3387 0 : MSSelection mssel;
3388 0 : mssel.setSpwExpr(spws);
3389 :
3390 : /*
3391 : cout << "mssel.getSpwList(mssel_p) = " << mssel.getSpwList(mssel_p) << endl;
3392 : cout << "mssel.getChanList(mssel_p) = " << mssel.getChanList(mssel_p) << endl; cout << "Vector<Int>() = " << Vector<Int>() << endl;
3393 : */
3394 :
3395 : // Use getChanList (column 0 for spw) because it is
3396 : // more reliable about the number and order of specified spws.
3397 : // return mssel.getSpwList(mssel_p);
3398 0 : Matrix<Int> chanmat=mssel.getChanList(mssel_p);
3399 0 : if (chanmat.nelements()>0)
3400 0 : return chanmat.column(0);
3401 : else
3402 0 : return Vector<Int>();
3403 :
3404 0 : }
3405 :
3406 0 : Matrix<Int> Calibrater::getChanIdx(const String& spw) {
3407 :
3408 0 : MSSelection mssel;
3409 0 : mssel.setSpwExpr(spw);
3410 :
3411 0 : return mssel.getChanList(mssel_p);
3412 :
3413 0 : }
3414 :
3415 :
3416 : // Query apply types if we are calibrating the weights
3417 0 : Bool Calibrater::calWt() {
3418 :
3419 0 : Int napp(vc_p.nelements());
3420 : // Return true as soon as we find a type which is cal'ing wts
3421 0 : for (Int iapp=0;iapp<napp;++iapp)
3422 0 : if (vc_p[iapp] && vc_p[iapp]->calWt())
3423 0 : return true;
3424 :
3425 : // None cal'd weights, so return false
3426 0 : return false;
3427 :
3428 : }
3429 :
3430 0 : Bool Calibrater::ok() {
3431 :
3432 0 : if ((simdata_p ||
3433 0 : (ms_p && mssel_p))
3434 0 : && ve_p) {
3435 0 : return true;
3436 : }
3437 : else {
3438 0 : logSink() << "Calibrater is not yet initialized" << LogIO::POST;
3439 0 : return false;
3440 : }
3441 : }
3442 : // The standard solving mechanism
3443 0 : casacore::Bool Calibrater::genericGatherAndSolve()
3444 : {
3445 :
3446 : #ifdef _OPENMP
3447 0 : Double Tsetup(0.0),Tgather(0.0),Tsolve(0.0),Tadd(0.0);
3448 0 : Double time0=omp_get_wtime();
3449 : #endif
3450 :
3451 : // Create record to store results
3452 0 : resRec_ = Record();
3453 0 : Record attemptRec;
3454 0 : attemptRec = Record();
3455 :
3456 : // Condition solint values
3457 0 : svc_p->reParseSolintForVI2();
3458 :
3459 : // Organize VI2 layers for solving:
3460 0 : CalSolVi2Organizer vi2org;
3461 :
3462 : //-------------------------------------------------
3463 : // NB: Populating the vi2org should probably be delegated to the svc_p
3464 : // since that is where most of the required (non-data) info is,
3465 : // and then some calibration types may introduce layers that are
3466 : // very specific to them (e.g., SD data selection, etc.)
3467 : // e.g., replace the following with something like:
3468 : //
3469 : // svc_p->formVi2LayerFactories(vi2org,mssel_p,ve_p);
3470 : // ssvp_p // a simdata version
3471 : //
3472 :
3473 : // Add the (bottom) data access layer
3474 0 : if (simdata_p)
3475 : // Simulated data (used for testing)
3476 0 : vi2org.addSimIO(ssvp_p);
3477 : else
3478 : {
3479 : // Real (selected) data in an MS (channel selection handled by addDiskIO method)
3480 : // The iteration time-interval is the solution interval
3481 0 : vi2org.addDiskIO(mssel_p,svc_p->solTimeInterval(),
3482 0 : svc_p->combobs(),svc_p->combscan(),
3483 0 : svc_p->combfld(),svc_p->combspw(),
3484 : true, // use MSIter2
3485 0 : frequencySelections_p); // Tell VI2 factory about the freq selection!
3486 : }
3487 :
3488 : // Add ad hoc SD section layer (e.g., OTF select of raster boundaries, etc.)
3489 : // only double circle gain calibration is implemented
3490 0 : bool SD = svc_p->longTypeName().startsWith("SDGAIN_OTFD");
3491 0 : if (SD) {
3492 0 : vi2org.addCalFilter(calFilterConfig_p);
3493 : }
3494 :
3495 : // Add pre-cal layer, using the VisEquation
3496 : // (include control for corr-dep flags)
3497 0 : vi2org.addCalForSolving(*ve_p,corrDepFlags_);
3498 :
3499 :
3500 : // Add the freq-averaging layer, if needed
3501 : //cout << "svc_p->fintervalChV() = " << svc_p->fintervalChV() << endl;
3502 : //cout << "svc_p->fsolint() = " << svc_p->fsolint() << endl;
3503 : // TBD: improve the following logic with new method in SVC...
3504 0 : if (!svc_p->freqDepMat() || // entirely unchannelized cal OR
3505 0 : (svc_p->freqDepPar() && // channelized par and
3506 0 : svc_p->fsolint()!="none" && // some partial channel averaging
3507 0 : anyGT(svc_p->fintervalChV(),1.0)) // explicity specified (non-trivially)
3508 : ) {
3509 0 : Vector<Int> chanbin(svc_p->fintervalChV().nelements());
3510 0 : convertArray(chanbin,svc_p->fintervalChV());
3511 0 : vi2org.addChanAve(chanbin);
3512 0 : }
3513 :
3514 : // Add the time-averaging layer, if needed
3515 : // NB: There is some evidence that doing time-averaging _after_
3516 : // channel averaging may be more expensive then doing it _before_,
3517 : // but this ensures that the meta-info (e.g. time, timeCentroid, etc.)
3518 : // appear correctly in the VB2 accessed at this scope.
3519 : // The problem is that AveragingTvi2 has not explicitly
3520 : // implemented all of the relevant data accessors; it just
3521 : // assumes---incorrectly---that you are going to use its VB2.
3522 : // Making AveragingTvi2 the top layer ensures you do use it.
3523 : // (Hard-wired use of AveragingTvi2 in MSTransform set it up
3524 : // as the top [last] layer, and has been well-tested...)
3525 :
3526 : //cout << "svc_p->solint() = " << svc_p->solint() << endl;
3527 : //cout << "svc_p->solTimeInterval() = " << svc_p->solTimeInterval() << endl;
3528 : //cout << "svc_p->preavg() = " << svc_p->preavg() << endl;
3529 0 : if (!svc_p->solint().contains("int") && // i.e., not "int")
3530 0 : svc_p->preavg() != -DBL_MAX) {
3531 0 : Float avetime(svc_p->solTimeInterval()); // avetime is solint, nominally
3532 : // Use preavg instead, if...
3533 0 : if (svc_p->preavg()>FLT_EPSILON && // ...set meaningfully
3534 0 : svc_p->preavg()<svc_p->solTimeInterval()) // ...and less than solint
3535 0 : avetime=svc_p->preavg();
3536 0 : vi2org.addTimeAve(avetime); // use min of solint and preavg here!
3537 : }
3538 : // small@jive.eu (2021-12-17): Currently we only handle the corrcomb
3539 : // cases of "all" and "none"; there will be a need to extend this, but the
3540 : // PolAverageTVILayerFactory that underlies this feature will also need to be extended
3541 : // to make that possible
3542 0 : if (svc_p->corrcomb().contains("all")) {
3543 : //cerr << "Calibrater::genericGatherAndSolve(): Combining correlations!" << endl;
3544 0 : vi2org.addCorrCombine();
3545 : }
3546 : //else {
3547 : // cerr << "Calibrater::genericGatherAndSolve(): Not combining correlations!" << endl;
3548 : //}
3549 :
3550 : // vi2org should be fully configured at this point
3551 : //-------------------------------------------------
3552 :
3553 :
3554 : // Form the VI2 to drive data iteration below
3555 0 : vi::VisibilityIterator2& vi(vi2org.makeFullVI());
3556 : // cout << "VI Layers: " << vi.ViiType() << endl;
3557 :
3558 : // Establish output freq meta data prior to solving (CAS-12735 and related)
3559 :
3560 : // Discern (even implicitly-)selected spw subset
3561 0 : set<uInt> selspwset;
3562 0 : if (!simdata_p) {
3563 :
3564 : // We can't proceed rationally if we have more than one freqSelection (multiple MSs)
3565 : // (this is a sort-of out-of-the-blue test here....)
3566 0 : AlwaysAssert(frequencySelections_p->size()<2,AipsError);
3567 :
3568 : // TBD can we use the msmc_p? (or maybe it is for the whole MS, not just the selected one?)
3569 0 : MSMetaData msmdtmp(mssel_p,50.0);
3570 0 : selspwset=msmdtmp.getSpwIDs();
3571 0 : }
3572 :
3573 : // Extract selected spw list as a Vector<uInt>
3574 0 : Vector<uInt> selspwlist;
3575 0 : if (!simdata_p && selspwset.size()>0) {
3576 0 : Vector<uInt> vtmp(selspwset.begin(),selspwset.size(),0);
3577 0 : selspwlist.reference(vtmp);
3578 0 : } else {
3579 : // All spws in the MS (this is a last resort!)
3580 : // NB: This is used in simdata_p=True context!
3581 0 : selspwlist.resize(msmc_p->nSpw());
3582 0 : indgen(selspwlist);
3583 : }
3584 : // cout << "selspwlist = " << selspwlist << endl;
3585 :
3586 : // Delegate freq meta calculation to the SVC (wherein this info is stored)
3587 0 : svc_p->discernAndSetSolnFrequencies(vi,selspwlist);
3588 :
3589 :
3590 : // Access to the net VB2 for forming SolveDataBuffers in the SDBList below
3591 : // NB: this assumes that the internal address of the VI2's VB2 will never change!
3592 0 : vi::VisBuffer2 *vb = vi.getImpl()->getVisBuffer();
3593 :
3594 : // VI2 should now be ready....
3595 :
3596 : // Discern number of Solutions and how many VI2 chunks per each
3597 0 : Int nSol(0);
3598 0 : Vector<Int> nChPSol(1000,0); // nominal size; will enlarge as needed, and trim at end
3599 :
3600 : {
3601 : //cout << "Counting chunks and solutions..." << endl;
3602 0 : uInt isol(0);
3603 : // uInt ichk(0);
3604 0 : for (vi.originChunks();vi.moreChunks();vi.nextChunk()) {
3605 : // NB: this loop NEVER touches the the subChunks, since
3606 : // this would cause data I/O, cal, averaging, etc.
3607 :
3608 : // int frame=vi.getImpl()->getReportingFrameOfReference();
3609 : // cout << "ichk=" << ichk << " freqs= " << vi.getImpl()->getFrequencies(-1.0e0,frame,selspwlist(ispw),0);
3610 :
3611 :
3612 : // Enlarge nChPSol Vector, if necessary
3613 0 : if (isol>nChPSol.nelements()-1) {
3614 : //cout << " ...At isol=" << isol << ", resizing nchPSol: " << nChPSol.nelements() << "-->";
3615 0 : nChPSol.resize(isol*2,True); // double length, copy existing elements
3616 0 : nChPSol(Slice(isol,isol,1)).set(0); // init new elements
3617 : //cout << nChPSol.nelements() << endl;
3618 : //cout << "sol counting: " << isol << " " << nChPSol << " (" << ntrue(nChPSol>0) << "/" << nChPSol.nelements() << ")" << endl;
3619 : }
3620 :
3621 : // incr chunk counter for current solution
3622 0 : ++nChPSol(isol);
3623 : //cout << "sol counting: " << isol << " " << nChPSol(Slice(0,min(isol+2,nChPSol.nelements()),1))
3624 : // << " (" << ntrue(nChPSol>0) << "/" << nChPSol.nelements() << ")"
3625 : // << " keyCh=" << vi.keyChange()
3626 : // << endl;
3627 :
3628 : // next chunk is nominally next solution
3629 0 : ++isol;
3630 :
3631 : // If combining spw or field, and the next chunk changes by
3632 : // DATA_DESC_ID or FIELD_ID (and _not_ TIME, which is
3633 : // changing slower than these when combine is on),
3634 : // then the next chunk should be counted in the same soln
3635 0 : if ( (svc_p->combspw() && vi.keyChange()=="DATA_DESC_ID") ||
3636 0 : (svc_p->combfld() && vi.keyChange()=="FIELD_ID") )
3637 0 : --isol; // next one is the same solution interval
3638 :
3639 : }
3640 : // Trim list to realized length, which is the number of solutions we'll
3641 : // attempt below
3642 0 : nSol=ntrue(nChPSol>0);
3643 0 : nChPSol.resize(nSol,true);
3644 : //cout << "nSol =" << nSol << endl;
3645 : //cout << "nChPSol=" << nChPSol << endl;
3646 :
3647 : /* from SVC...
3648 : if (svc_p->combobs())
3649 : logSink() << "Combining observation Ids." << LogIO::POST;
3650 : if (svc_p->combscan())
3651 : logSink() << "Combining scans." << LogIO::POST;
3652 : if (svc_p->combspw())
3653 : logSink() << "Combining spws: " << spwlist << " -> " << spwlab << LogIO::POST;
3654 : if (svc_p->combfld())
3655 : logSink() << "Combining fields." << LogIO::POST;
3656 : */
3657 :
3658 0 : logSink() << "For solint = " << svc_p->solint() << ", found "
3659 : << nSol << " solution intervals."
3660 0 : << LogIO::POST;
3661 : }
3662 :
3663 : // throw(AipsError("EARLY ESCAPE!!"));
3664 :
3665 : // Create the output caltable
3666 : // (this version doesn't revise frequencies)
3667 0 : svc_p->createMemCalTable2();
3668 :
3669 0 : Vector<Float> spwwts(msmc_p->nSpw(),-1.0);
3670 0 : 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);
3671 :
3672 : // Start Collecting counts
3673 0 : CalCounts* calCounts = new CalCounts();
3674 0 : calCounts->initCounts(msmc_p->nSpw(),msmc_p->nAnt(), svc_p->nPar());
3675 : // Set up results record
3676 0 : std::map<Int, std::map<String, Vector<Int>>> resultMap;
3677 :
3678 :
3679 : #ifdef _OPENMP
3680 0 : Tsetup+=(omp_get_wtime()-time0);
3681 : #endif
3682 :
3683 0 : Int nGood(0);
3684 0 : vi.originChunks();
3685 0 : Int nGlobalChunks=0; // counts VI chunks globally
3686 0 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
3687 :
3688 : #ifdef _OPENMP
3689 0 : time0=omp_get_wtime();
3690 : #endif
3691 :
3692 : // Data will accumulate here
3693 0 : SDBList sdbs;
3694 :
3695 : // Gather the chunks/VBs for this solution
3696 : // Solution boundaries will ALWAYS occur on chunk boundaries,
3697 : // though some chunk boundaries will be ignored in the
3698 : // combine='spw' or 'field' context
3699 :
3700 0 : for (Int ichunk=0; // count chunks in this solution
3701 0 : ichunk<nChPSol(isol)&&vi.moreChunks(); // while more chunks needed _and_ available
3702 0 : ++ichunk,vi.nextChunk()) { // advance to next chunk
3703 :
3704 : // Global chunk counter
3705 0 : ++nGlobalChunks;
3706 :
3707 : // Loop over VB2s in this chunk
3708 : // (we get more then one when preavg<solint)
3709 0 : Int ivb(0);
3710 0 : for (vi.origin();
3711 0 : vi.more();
3712 0 : ++ivb,vi.next()) {
3713 :
3714 : // Add this VB to the SDBList
3715 : #ifdef _OPENMP
3716 0 : Double Tadd0=omp_get_wtime();
3717 : #endif
3718 :
3719 0 : sdbs.add(*vb);
3720 :
3721 : #ifdef _OPENMP
3722 0 : Tadd+=(omp_get_wtime()-Tadd0);
3723 : #endif
3724 :
3725 : // Keep track of spws seen but not included in solving
3726 0 : Int ispw=vb->spectralWindows()(0);
3727 0 : if (spwwts(ispw)<0) spwwts(ispw)=0.0f;
3728 0 : spwwts(ispw)+=sum(vb->weightSpectrum());
3729 :
3730 : /*
3731 : Double modv(86400.0);
3732 : cout.precision(7);
3733 : cout << "isol=" << isol
3734 : << " ichk="<<ichunk
3735 : << " ivb="<<ivb
3736 : << " sc="<<vb->scan()(0)
3737 : << " fld="<<vb->fieldId()(0)
3738 : << " spw="<<vb->spectralWindows()(0)
3739 : << " nr="<<vb->nRows()
3740 : << " nch="<<vb->nChannels() << flush;
3741 : cout << " f="<<mean(vb->getFrequencies(0))/1e9
3742 : << "->"<< mean(sdbs.freqs())/1e9
3743 : << " t="<<fmod(mean(vb->time()),modv)
3744 : << "->"<< fmod(sdbs.aggregateTime(),modv)
3745 : << " tC="<<fmod(mean(vb->timeCentroid()),modv)
3746 : << "->"<< fmod(sdbs.aggregateTimeCentroid(),modv)
3747 : << " wt="<<sum(vb->weightSpectrum())
3748 : << " nSDB=" << sdbs.nSDB()
3749 : << " nGlChks=" << nGlobalChunks
3750 : << " (keych=" << vi.keyChange() << ")"
3751 : << endl;
3752 : //*/
3753 :
3754 : } // VI2 subchunks (VB2s)
3755 : } // VI2 chunks
3756 :
3757 : // Which spw is this?
3758 0 : Int thisSpw(sdbs.aggregateSpw());
3759 :
3760 : // Expecting a solution
3761 0 : nexp(thisSpw)+=1;
3762 : // Initialize the antennaMap_ in svc
3763 0 : svc_p->clearMap();
3764 : // Get expected and data unflagged accumulation
3765 0 : svc_p->expectedUnflagged(sdbs);
3766 :
3767 :
3768 : #ifdef _OPENMP
3769 0 : Tgather+=(omp_get_wtime()-time0);
3770 0 : time0=omp_get_wtime();
3771 : #endif
3772 :
3773 :
3774 0 : if (sdbs.Ok()) {
3775 :
3776 : // Some unflagged data, so Attempting a solution
3777 0 : natt(thisSpw)+=1;
3778 :
3779 : // make phase- or amp-only, if necessary
3780 0 : sdbs.enforceAPonData(svc_p->apmode());
3781 :
3782 : // zero cross-hand weights, if necessary
3783 0 : sdbs.enforceSolveWeights(svc_p->phandonly());
3784 :
3785 : // Synchronize meta-info in SVC
3786 0 : svc_p->syncSolveMeta(sdbs);
3787 :
3788 : // Set or verify freqs in the caltable
3789 : //svc_p->setOrVerifyCTFrequencies(thisSpw);
3790 0 : svc_p->setCTFrequencies(thisSpw);
3791 :
3792 : // Size the solvePar arrays inside SVC
3793 : // (smart: if freqDepPar()=F, uses 1)
3794 : // returns the number of channel solutions to iterate over
3795 : //Int nChanSol=svc_p->sizeSolveParCurrSpw(sdbs.nChannels());
3796 0 : Int nChanSol=svc_p->sizeSolveParCurrSpw((svc_p->freqDepPar() ? sdbs.nChannels() : 1));
3797 :
3798 0 : if (svc_p->useGenericSolveOne()) {
3799 :
3800 : // We'll use the generic solver
3801 0 : VisCalSolver2 vcs(svc_p->solmode(),svc_p->rmsthresh());
3802 :
3803 : // Guess from the data
3804 0 : svc_p->guessPar(sdbs,corrDepFlags_);
3805 :
3806 0 : Bool totalGoodSol(False); // Will be set True if any channel is good
3807 : //for (Int ich=0;ich<nChanSol;++ich) {
3808 0 : for (Int ich=nChanSol-1;ich>-1;--ich) {
3809 0 : svc_p->markTimer();
3810 0 : svc_p->focusChan()=ich;
3811 :
3812 : // Execute the solve
3813 0 : Bool goodsol=vcs.solve(*ve_p,*svc_p,sdbs);
3814 :
3815 0 : if (goodsol) {
3816 0 : totalGoodSol=True;
3817 :
3818 0 : svc_p->formSolveSNR();
3819 0 : svc_p->applySNRThreshold();
3820 : }
3821 : else
3822 0 : svc_p->currMetaNote();
3823 :
3824 : // Record solution in its channel, good or bad
3825 0 : if (svc_p->freqDepPar())
3826 0 : svc_p->keep1(ich);
3827 :
3828 : } // ich
3829 :
3830 0 : if (totalGoodSol) {
3831 : // Keep this good solution, and count it
3832 0 : svc_p->keepNCT();
3833 0 : ++nGood;
3834 0 : nsuc(thisSpw)+=1;
3835 : }
3836 :
3837 0 : } // useGenericSolveOne
3838 : else {
3839 : // Use self-directed individual solve
3840 : // TBD: return T/F for goodness?
3841 : //cout << "Calling selfSolveOne()!!!!!!" << endl;
3842 0 : svc_p->selfSolveOne(sdbs);
3843 :
3844 : // Keep this good solution, and count it
3845 0 : svc_p->keepNCT();
3846 0 : ++nGood;
3847 0 : nsuc(thisSpw)+=1;
3848 : }
3849 :
3850 : } // sdbs.Ok()
3851 : else {
3852 : // Synchronize meta-info in SVC
3853 0 : svc_p->syncSolveMeta(sdbs);
3854 0 : cout << "Found no unflagged data at:";
3855 0 : svc_p->currMetaNote();
3856 : }
3857 : //cout << endl;
3858 : // Get all the antenna value counts
3859 0 : resultMap = svc_p->getAntennaMap();
3860 0 : calCounts->addAntennaCounts(thisSpw,msmc_p->nAnt(), svc_p->nPar(), resultMap);
3861 :
3862 : #ifdef _OPENMP
3863 0 : Tsolve+=(omp_get_wtime()-time0);
3864 : #endif
3865 :
3866 : // throw(AipsError("EARLY ESCAPE!!"));
3867 :
3868 0 : } // isol
3869 :
3870 : // Report nGood to logger
3871 : logSink() << " Found good "
3872 0 : << svc_p->typeName() << " solutions in "
3873 : << nGood << " solution intervals."
3874 0 : << LogIO::POST;
3875 :
3876 : #ifdef _OPENMP
3877 : #ifdef REPORT_CAL_TIMING
3878 : cout << "Calibrater::genericGatherAndSolve Timing: " << endl;
3879 : cout << " setup=" << Tsetup
3880 : << " gather=" << Tgather
3881 : << " (SDBadd=" << Tadd << ")"
3882 : << " solve=" << Tsolve
3883 : << " total=" << Tsetup+Tgather+Tsolve
3884 : // << " tick=" << omp_get_wtick()
3885 : << endl;
3886 : #endif
3887 : #endif
3888 :
3889 : // Report spws that were seen but not solved
3890 0 : Vector<Bool> unsolspw=(spwwts==0.0f);
3891 0 : summarize_uncalspws(unsolspw, "solv");
3892 :
3893 : // throw(AipsError("EARLY ESCAPE!!"));
3894 :
3895 0 : if (nGood>0) {
3896 0 : if (svc_p->typeName()!="BPOLY") { // needed?
3897 : // apply refant, etc.
3898 0 : svc_p->clearRefantMap();
3899 0 : svc_p->globalPostSolveTinker();
3900 :
3901 : // write to disk
3902 0 : svc_p->storeNCT();
3903 : }
3904 : }
3905 : else {
3906 : logSink() << "No output calibration table written."
3907 0 : << LogIO::POST;
3908 : }
3909 :
3910 : // Collect and update the refants
3911 0 : std::map<Int, std::map<Int, Int>> refantMap;
3912 0 : refantMap = svc_p->getRefantMap();
3913 0 : calCounts->updateRefants(msmc_p->nSpw(), msmc_p->nAnt(), refantMap);
3914 : // Compile all accumulated counts into a record
3915 0 : resRec_ = calCounts->makeRecord(msmc_p->nAnt(), vi.nPolarizationIds());
3916 :
3917 : // print a matrix to the logger
3918 0 : logSink() << "----- PER ANTENNA INFO -----" << LogIO::POST;
3919 : // print the antenna list
3920 0 : logSink() << " ";
3921 0 : for (Int ant=0; ant<msmc_p->nAnt(); ant++) {
3922 0 : logSink() << " ANT: " << ant << " ";
3923 : }
3924 0 : logSink() << LogIO::POST;
3925 :
3926 0 : for (Int spw=0; spw < msmc_p->nSpw(); spw++) {
3927 :
3928 0 : logSink() << LogIO::NORMAL << "SPW: " << spw;
3929 :
3930 0 : for (Int ant=0; ant < msmc_p->nAnt(); ant++) {
3931 0 : logSink() << " " << calCounts->antMapVal(spw, ant, "above_minsnr") << " ";
3932 : }
3933 0 : logSink() << LogIO::POST;
3934 : }
3935 :
3936 0 : logSink() << "----- PER SPW INFO -----" << LogIO::POST;
3937 : // print the result fields
3938 0 : logSink() << " ";
3939 0 : logSink() << " expected data_unflagged above_minblperant above_minsnr";
3940 0 : logSink() << LogIO::POST;
3941 0 : for (Int spw = 0; spw < msmc_p->nSpw(); spw++) {
3942 0 : logSink() << LogIO::NORMAL << "SPW: " << spw << " ";
3943 0 : logSink() << calCounts->spwMapVal(spw, "expected") << " ";
3944 0 : logSink() << calCounts->spwMapVal(spw, "data_unflagged") << " ";
3945 0 : logSink() << calCounts->spwMapVal(spw, "above_minblperant") << " ";
3946 0 : logSink() << calCounts->spwMapVal(spw, "above_minsnr");
3947 0 : logSink() << LogIO::POST;
3948 : }
3949 :
3950 0 : logSink() << "----- GLOBAL INFO -----" << LogIO::POST;
3951 0 : logSink() << "expected data_unflagged above_minblperant above_minsnr";
3952 0 : logSink() << LogIO::POST;
3953 0 : logSink() << calCounts->totalMapVal("expected") << " ";
3954 0 : logSink() << calCounts->totalMapVal("data_unflagged") << " ";
3955 0 : logSink() << calCounts->totalMapVal("above_minblperant") << " ";
3956 0 : logSink() << calCounts->totalMapVal("above_minsnr");
3957 0 : logSink() << LogIO::POST;
3958 :
3959 : // Fill activity record
3960 :
3961 0 : actRec_ = Record();
3962 0 : actRec_.define("origin","Calibrater::genericGatherAndSolve");
3963 0 : actRec_.define("nExpected",nexp);
3964 0 : actRec_.define("nAttempt",natt);
3965 0 : actRec_.define("nSucceed",nsuc);
3966 :
3967 : //cout << nexp << ", " << natt << ", " << nsuc << endl;
3968 :
3969 : {
3970 0 : Record solveRec=svc_p->solveActionRec();
3971 0 : if (solveRec.nfields()>0)
3972 0 : actRec_.merge(solveRec);
3973 0 : }
3974 :
3975 : // Reach here, all is good
3976 0 : return True;
3977 :
3978 0 : }
3979 :
3980 :
3981 0 : void Calibrater::writeHistory(LogIO& /*os*/, Bool /*cliCommand*/)
3982 : {
3983 : // Disabled 2016/04/19: avoid direct MS.HISTORY updates from
3984 : // below the python level, FOR NOW
3985 :
3986 0 : return;
3987 : /*
3988 : if (!historytab_p.isNull()) {
3989 : if (histLockCounter_p == 0) {
3990 : historytab_p.lock(false);
3991 : }
3992 : ++histLockCounter_p;
3993 :
3994 : os.postLocally();
3995 : if (cliCommand) {
3996 : hist_p->cliCommand(os);
3997 : } else {
3998 : hist_p->addMessage(os);
3999 : }
4000 :
4001 : if (histLockCounter_p == 1) {
4002 : historytab_p.unlock();
4003 : }
4004 : if (histLockCounter_p > 0) {
4005 : --histLockCounter_p;
4006 : }
4007 : } else {
4008 : os << LogIO::SEVERE << "calibrater is not yet initialized" << LogIO::POST;
4009 : }
4010 : */
4011 : }
4012 :
4013 0 : void Calibrater::setCalFilterConfiguration(String const &type,
4014 : Record const &config) {
4015 : // currently only SDDoubleCircleGainCal requires data filtering
4016 0 : if (type.startsWith("SDGAIN_OTFD")) {
4017 0 : calFilterConfig_p.define("mode", "SDGAIN_OTFD");
4018 0 : if (config.isDefined("smooth")) {
4019 0 : calFilterConfig_p.define("smooth", config.asBool("smooth"));
4020 : }
4021 0 : if (config.isDefined("radius")) {
4022 0 : calFilterConfig_p.define("radius", config.asString("radius"));
4023 : }
4024 : }
4025 0 : }
4026 :
4027 : // *********************************************
4028 : // OldCalibrater implementations that use vs_p
4029 :
4030 0 : OldCalibrater::OldCalibrater():
4031 : Calibrater(),
4032 0 : vs_p(0),
4033 0 : rawvs_p(0)
4034 : {
4035 : // cout << "This is the OLD VI2-aware Calibrater" << endl;
4036 0 : }
4037 :
4038 0 : OldCalibrater::OldCalibrater(String msname):
4039 : Calibrater(msname),
4040 0 : vs_p(0),
4041 0 : rawvs_p(0)
4042 : {
4043 0 : }
4044 :
4045 : /*
4046 : OldCalibrater::OldCalibrater(const OldCalibrater & other) :
4047 : Calibrater(other)
4048 : {
4049 : operator=(other);
4050 : }
4051 :
4052 : OldCalibrater& OldCalibrater::operator=(const OldCalibrater & other)
4053 : {
4054 : Calibrater::operator=(other); // copy parental units
4055 : vs_p=other.vs_p;
4056 : rawvs_p=other.rawvs_p;
4057 : return *this;
4058 : }
4059 : */
4060 :
4061 0 : OldCalibrater::~OldCalibrater()
4062 : {
4063 0 : OldCalibrater::cleanupVisSet();
4064 0 : }
4065 :
4066 :
4067 : // Select data (using MSSelection syntax)
4068 0 : void OldCalibrater::selectvis(const String& time,
4069 : const String& spw,
4070 : const String& scan,
4071 : const String& field,
4072 : const String& intent,
4073 : const String& obsIDs,
4074 : const String& baseline,
4075 : const String& uvrange,
4076 : const String& chanmode,
4077 : const Int& nchan,
4078 : const Int& start,
4079 : const Int& step,
4080 : const MRadialVelocity& mStart,
4081 : const MRadialVelocity& mStep,
4082 : const String& msSelect)
4083 : {
4084 : // Define primary measurement set selection criteria
4085 : // Inputs:
4086 : // time
4087 : // spw
4088 : // scan
4089 : // field
4090 : // intent
4091 : // obsIDs
4092 : // baseline
4093 : // uvrange
4094 : // chanmode const String& Frequency/velocity selection mode
4095 : // ("channel", "velocity" or
4096 : // "opticalvelocity")
4097 : // nchan const Int& No of channels to select
4098 : // start const Int& Start channel to select
4099 : // step const Int& Channel increment
4100 : // mStart const MRadialVelocity& Start radial vel. to select
4101 : // mStep const MRadialVelocity& Radial velocity increment
4102 : // msSelect const String& MS selection string (TAQL)
4103 : // Output to private data:
4104 : //
4105 0 : logSink() << LogOrigin("Calibrater","selectvis") << LogIO::NORMAL3;
4106 :
4107 : try {
4108 :
4109 : /*
4110 : cout << "time = " << time << " " << time.length() <<endl;
4111 : cout << "spw = " << spw << " " << spw.length() <<endl;
4112 : cout << "scan = " << scan << " " << scan.length() <<endl;
4113 : cout << "field = " << field << " " << field.length() <<endl;
4114 : cout << "baseline = " << baseline << " " << baseline.length() << endl;
4115 : cout << "uvrange = " << uvrange << " " << uvrange.length() << endl;
4116 : */
4117 :
4118 0 : logSink() << "Selecting data" << LogIO::POST;
4119 :
4120 : // Apply selection to the original MeasurementSet
4121 0 : logSink() << "Performing selection on MeasurementSet" << endl;
4122 :
4123 0 : if (mssel_p) {
4124 0 : delete mssel_p;
4125 0 : mssel_p=0;
4126 : };
4127 :
4128 : // Report non-trivial user selections
4129 0 : if (time!="")
4130 0 : logSink() << " Selecting on time: '" << time << "'" << endl;
4131 0 : if (spw!="")
4132 0 : logSink() << " Selecting on spw: '" << spw << "'" << endl;
4133 0 : if (scan!="")
4134 0 : logSink() << " Selecting on scan: '" << scan << "'" << endl;
4135 0 : if (field!="")
4136 0 : logSink() << " Selecting on field: '" << field << "'" << endl;
4137 0 : if (intent!="")
4138 0 : logSink() << " Selecting on intent: '" << intent << "'" << endl;
4139 0 : if(obsIDs != "")
4140 0 : logSink() << " Selecting by observation IDs: '" << obsIDs << "'" << endl;
4141 0 : if (baseline!="")
4142 0 : logSink() << " Selecting on antenna/baseline: '" << baseline << "'" << endl;
4143 0 : if (uvrange!="")
4144 0 : logSink() << " Selecting on uvrange: '" << uvrange << "'" << endl;
4145 0 : if (msSelect!="")
4146 0 : logSink() << " Selecting with TaQL: '" << msSelect << "'" << endl;
4147 0 : logSink() << LogIO::POST;
4148 :
4149 :
4150 : // Assume no selection, for starters
4151 : // gmoellen 2012/01/30 mssel_p = new MeasurementSet(sorted, ms_p);
4152 0 : mssel_p = new MeasurementSet(*ms_p);
4153 :
4154 : // Apply user-supplied selection
4155 0 : Bool nontrivsel=false;
4156 : // gmoellen 2012/01/30 nontrivsel= mssSetData(MeasurementSet(sorted, ms_p),
4157 :
4158 : // Ensure use of a fresh MSSelection object
4159 0 : if (mss_p) { delete mss_p; mss_p=NULL; }
4160 0 : mss_p=new MSSelection();
4161 0 : nontrivsel= mssSetData(*ms_p,
4162 0 : *mssel_p,"",
4163 : time,baseline,
4164 : field,spw,
4165 : uvrange,msSelect,
4166 : "",scan,"",intent, obsIDs,mss_p);
4167 :
4168 : // Keep any MR status for the MS
4169 0 : mssel_p->setMemoryResidentSubtables(ms_p->getMrsEligibility());
4170 :
4171 : // If non-trivial MSSelection invoked and nrow reduced:
4172 0 : if(nontrivsel && mssel_p->nrow()<ms_p->nrow()) {
4173 :
4174 : // Escape if no rows selected
4175 0 : if (mssel_p->nrow()==0)
4176 0 : throw(AipsError("Specified selection selects zero rows!"));
4177 :
4178 : // ...otherwise report how many rows are selected
4179 0 : logSink() << "By selection " << ms_p->nrow()
4180 0 : << " rows are reduced to " << mssel_p->nrow()
4181 0 : << LogIO::POST;
4182 : }
4183 : else {
4184 : // Selection did nothing:
4185 0 : logSink() << "Selection did not drop any rows" << LogIO::POST;
4186 : }
4187 :
4188 : // Now, re-create the associated VisSet
4189 0 : if(vs_p) delete vs_p; vs_p=0;
4190 0 : Block<int> sort(0);
4191 0 : Matrix<Int> noselection;
4192 : // gmoellen 2012/01/30 vs_p = new VisSet(*mssel_p,sort,noselection);
4193 0 : vs_p = new VisSet(*mssel_p,sort,noselection,false,0.0,false,false);
4194 0 : AlwaysAssert(vs_p, AipsError);
4195 :
4196 : // Attempt to use MSSelection for channel selection
4197 : // if user not using the old way
4198 0 : if (chanmode=="none") {
4199 0 : selectChannel(spw);
4200 : }
4201 : else {
4202 : // Reluctantly use the old-fashioned way
4203 : logSink() << LogIO::WARN
4204 : << "You have used the old-fashioned mode parameter" << endl
4205 : << "for channel selection. It still works, for now," << endl
4206 : << "but this will be eliminated in the near future." << endl
4207 : << "Please begin using the new channel selection" << endl
4208 0 : << "syntax in the spw parameter." << LogIO::POST;
4209 0 : selectChannel(chanmode,nchan,start,step,mStart,mStep);
4210 : }
4211 :
4212 0 : }
4213 0 : catch (MSSelectionError& x) {
4214 : // Re-initialize with the existing MS
4215 0 : logSink() << LogOrigin("Calibrater","selectvis",WHERE)
4216 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4217 0 : << LogIO::POST;
4218 : // jagonzal (CAS-4110): I guess it is not necessary to create these columns when the selection is empty
4219 0 : initialize(*ms_p,false,false,false);
4220 0 : throw(AipsError("Error in data selection specification: " + x.getMesg()));
4221 0 : }
4222 0 : catch (AipsError x) {
4223 : // Re-initialize with the existing MS
4224 0 : logSink() << LogOrigin("Calibrater","selectvis",WHERE)
4225 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4226 0 : << LogIO::POST;
4227 : // jagonzal (CAS-4110): I guess it is not necessary to create these columns when the selection is empty.
4228 0 : initialize(*ms_p,false,false,false);
4229 0 : throw(AipsError("Error in Calibrater::selectvis(): " + x.getMesg()));
4230 0 : }
4231 0 : };
4232 :
4233 :
4234 0 : Bool OldCalibrater::setapply (const String& type,
4235 : const Record& applypar)
4236 : {
4237 0 : logSink() << LogOrigin("Calibrater", "setapply(type, applypar)");
4238 :
4239 : // First try to create the requested VisCal object
4240 0 : VisCal *vc(NULL);
4241 :
4242 : try {
4243 :
4244 0 : if(!ok())
4245 0 : throw(AipsError("Calibrater not prepared for setapply."));
4246 :
4247 0 : String upType=type;
4248 0 : upType.upcase();
4249 :
4250 : logSink() << LogIO::NORMAL
4251 : << "Arranging to APPLY:"
4252 0 : << LogIO::POST;
4253 :
4254 : // Add a new VisCal to the apply list
4255 0 : vc = createVisCal(upType,*vs_p);
4256 :
4257 0 : vc->setApply(applypar);
4258 :
4259 : logSink() << LogIO::NORMAL << ". "
4260 0 : << vc->applyinfo()
4261 0 : << LogIO::POST;
4262 :
4263 0 : } catch (AipsError x) {
4264 : logSink() << LogIO::SEVERE << x.getMesg()
4265 : << " Check inputs and try again."
4266 0 : << LogIO::POST;
4267 0 : if (vc) delete vc;
4268 0 : throw(AipsError("Error in Calibrater::setapply."));
4269 : return false;
4270 0 : }
4271 :
4272 : // Creation apparently successful, so add to the apply list
4273 : // TBD: consolidate with above?
4274 : try {
4275 :
4276 0 : uInt napp=vc_p.nelements();
4277 0 : vc_p.resize(napp+1,false,true);
4278 0 : vc_p[napp] = vc;
4279 0 : vc=NULL;
4280 :
4281 : // Maintain sort of apply list
4282 0 : ve_p->setapply(vc_p);
4283 :
4284 0 : return true;
4285 :
4286 0 : } catch (AipsError x) {
4287 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4288 0 : << LogIO::POST;
4289 0 : if (vc) delete vc;
4290 0 : throw(AipsError("Error in Calibrater::setapply."));
4291 : return false;
4292 0 : }
4293 : return false;
4294 : }
4295 :
4296 : // Set up apply-able calibration via a Cal Library
4297 0 : Bool OldCalibrater::setcallib(Record callib) {
4298 :
4299 0 : logSink() << LogOrigin("Calibrater", "setcallib(callib)");
4300 :
4301 : // cout << "Calibrater::setcallib: callib.isFixed() = " << boolalpha << callib.isFixed() << endl;
4302 :
4303 0 : uInt ntab=callib.nfields();
4304 :
4305 : // cout << "callib.nfields() = " << ntab << endl;
4306 :
4307 : // Do some preliminary per-table verification
4308 0 : for (uInt itab=0;itab<ntab;++itab) {
4309 :
4310 0 : String tabname=callib.name(itab);
4311 :
4312 : // Insist that the table exists on disk
4313 0 : if (!Table::isReadable(tabname))
4314 0 : throw(AipsError("Caltable "+tabname+" does not exist."));
4315 :
4316 0 : }
4317 :
4318 : // Tables exist, so deploy them...
4319 :
4320 0 : for (uInt itab=0;itab<ntab;++itab) {
4321 :
4322 0 : String tabname=callib.name(itab);
4323 :
4324 : // Get the type from the table
4325 0 : String upType=calTableType(tabname);
4326 0 : upType.upcase();
4327 :
4328 : // Add table name to the record
4329 0 : Record thistabrec=callib.asrwRecord(itab);
4330 0 : thistabrec.define("tablename",tabname);
4331 :
4332 : // First try to create the requested VisCal object
4333 0 : VisCal *vc(NULL);
4334 :
4335 : try {
4336 :
4337 0 : if(!ok())
4338 0 : throw(AipsError("Calibrater not prepared for setapply."));
4339 :
4340 : logSink() << LogIO::NORMAL
4341 : << "Arranging to APPLY:"
4342 0 : << LogIO::POST;
4343 :
4344 : // Add a new VisCal to the apply list
4345 0 : vc = createVisCal(upType,*vs_p);
4346 :
4347 :
4348 : // ingest this table according to its callib
4349 0 : vc->setCallib(thistabrec,*mssel_p);
4350 :
4351 0 : } catch (AipsError x) {
4352 : logSink() << LogIO::SEVERE << x.getMesg()
4353 : << " Check inputs and try again."
4354 0 : << LogIO::POST;
4355 0 : if (vc) delete vc;
4356 0 : throw(AipsError("Error in Calibrater::setapply."));
4357 : return false;
4358 0 : }
4359 :
4360 : // Creation apparently successful, so add to the apply list
4361 : // TBD: consolidate with above?
4362 : try {
4363 :
4364 0 : uInt napp=vc_p.nelements();
4365 0 : vc_p.resize(napp+1,false,true);
4366 0 : vc_p[napp] = vc;
4367 0 : vc=NULL;
4368 :
4369 : // Maintain sort of apply list
4370 0 : ve_p->setapply(vc_p);
4371 :
4372 0 : } catch (AipsError x) {
4373 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4374 0 : << LogIO::POST;
4375 0 : if (vc) delete vc;
4376 0 : throw(AipsError("Error in Calibrater::setapply."));
4377 : return false;
4378 0 : }
4379 0 : }
4380 :
4381 : // All ok, if we get this far!
4382 0 : return true;
4383 :
4384 : }
4385 :
4386 :
4387 : // Set up apply-able calibration via a Cal Library
4388 0 : Bool OldCalibrater::setcallib2(Record callib, const casacore::MeasurementSet* ms) {
4389 :
4390 0 : logSink() << LogOrigin("OldCalibrater", "setcallib2(callib)");
4391 :
4392 : // cout << "Calibrater::setcallib2(callib) : " << boolalpha << callib << endl;
4393 :
4394 0 : uInt ntab=callib.nfields();
4395 :
4396 : // cout << "callib.nfields() = " << ntab << endl;
4397 :
4398 : // Do some preliminary per-table verification
4399 0 : for (uInt itab=0;itab<ntab;++itab) {
4400 :
4401 0 : String tabname=callib.name(itab);
4402 :
4403 : // Trap parang
4404 : // TBD...
4405 : // if (tabname=="<parang>")
4406 : // continue;
4407 :
4408 : // Insist that the table exists on disk
4409 0 : if (!Table::isReadable(tabname))
4410 0 : throw(AipsError("Caltable "+tabname+" does not exist."));
4411 :
4412 0 : }
4413 :
4414 : // Tables exist, so deploy them...
4415 :
4416 : // Local MS object for callib parsing (only)
4417 : // MeasurementSet lms(msname_p,Table::Update);
4418 : //cout << "OLD lms" << endl;
4419 :
4420 :
4421 : // Local const MS object for callib parsing (only)
4422 0 : const MeasurementSet *lmsp(0);
4423 0 : if (ms) {
4424 : // Use supplied MS (from outside), if specified...
4425 : // TBD: should we verify same base MS as ms_p/mssel_p?
4426 0 : lmsp=ms;
4427 : }
4428 : else {
4429 : // ...use internal one instead
4430 0 : lmsp=mssel_p;
4431 : }
4432 : // Reference for use below
4433 0 : const MeasurementSet &lms(*lmsp);
4434 :
4435 : // Get some global shape info:
4436 0 : Int MSnAnt = lms.antenna().nrow();
4437 0 : Int MSnSpw = lms.spectralWindow().nrow();
4438 :
4439 0 : for (uInt itab=0;itab<ntab;++itab) {
4440 :
4441 0 : String tabname=callib.name(itab);
4442 :
4443 : // Get the type from the table
4444 0 : String upType=calTableType(tabname);
4445 0 : upType.upcase();
4446 :
4447 : // Add table name to the record
4448 0 : Record thistabrec=callib.asrwRecord(itab);
4449 0 : thistabrec.define("tablename",tabname);
4450 :
4451 : // First try to create the requested VisCal object
4452 0 : VisCal *vc(NULL);
4453 :
4454 : try {
4455 :
4456 : // if(!ok())
4457 : // throw(AipsError("Calibrater not prepared for setapply."));
4458 :
4459 : logSink() << LogIO::NORMAL
4460 : << "Arranging to APPLY:"
4461 0 : << LogIO::POST;
4462 :
4463 : // Add a new VisCal to the apply list
4464 0 : vc = createVisCal(upType,msname_p,MSnAnt,MSnSpw);
4465 :
4466 : // ingest this table according to its callib
4467 0 : vc->setCallib(thistabrec,lms);
4468 :
4469 0 : } catch (AipsError x) {
4470 : logSink() << LogIO::SEVERE << x.getMesg()
4471 : << " Check inputs and try again."
4472 0 : << LogIO::POST;
4473 0 : if (vc) delete vc;
4474 0 : throw(AipsError("Error in Calibrater::callib2."));
4475 : return false;
4476 0 : }
4477 :
4478 : // Creation apparently successful, so add to the apply list
4479 : // TBD: consolidate with above?
4480 : try {
4481 :
4482 0 : uInt napp=vc_p.nelements();
4483 0 : vc_p.resize(napp+1,false,true);
4484 0 : vc_p[napp] = vc;
4485 0 : vc=NULL;
4486 :
4487 : // Maintain sort of apply list
4488 0 : ve_p->setapply(vc_p);
4489 :
4490 0 : } catch (AipsError x) {
4491 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4492 0 : << LogIO::POST;
4493 0 : if (vc) delete vc;
4494 0 : throw(AipsError("Error in Calibrater::setapply."));
4495 : return false;
4496 0 : }
4497 0 : }
4498 : // All ok, if we get this far!
4499 0 : return true;
4500 :
4501 : }
4502 :
4503 :
4504 0 : Bool OldCalibrater::setsolve (const String& type,
4505 : const Record& solvepar) {
4506 :
4507 : // Attempt to create the solvable object
4508 0 : SolvableVisCal *svc(NULL);
4509 : try {
4510 :
4511 0 : if(!ok())
4512 0 : throw(AipsError("Calibrater not prepared for setsolve."));
4513 :
4514 0 : String upType = type;
4515 0 : upType.upcase();
4516 :
4517 : // Clean out any old solve that was lying around
4518 0 : unsetsolve();
4519 :
4520 : logSink() << LogIO::NORMAL
4521 : << "Arranging to SOLVE:"
4522 0 : << LogIO::POST;
4523 :
4524 : // Create the new SolvableVisCal
4525 0 : svc = createSolvableVisCal(upType,*vs_p);
4526 0 : svc->setSolve(solvepar);
4527 :
4528 : logSink() << LogIO::NORMAL << ". "
4529 0 : << svc->solveinfo()
4530 0 : << LogIO::POST;
4531 :
4532 : // Creation apparently successful, keep it
4533 0 : svc_p=svc;
4534 0 : svc=NULL;
4535 :
4536 0 : return true;
4537 :
4538 0 : } catch (AipsError x) {
4539 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4540 0 : << LogIO::POST;
4541 0 : unsetsolve();
4542 0 : if (svc) delete svc;
4543 0 : throw(AipsError("Error in Calibrater::setsolve."));
4544 : return false;
4545 0 : }
4546 : return false;
4547 : }
4548 :
4549 0 : Bool OldCalibrater::correct(String mode)
4550 : {
4551 0 : logSink() << LogOrigin("Calibrater","correct") << LogIO::NORMAL;
4552 :
4553 0 : Bool retval = true;
4554 :
4555 : try {
4556 :
4557 : // make mode all-caps
4558 0 : String upmode=mode;
4559 0 : upmode.upcase();
4560 :
4561 : // If trialmode=T, only the flags will be set
4562 : // (and only written if not TRIAL)
4563 0 : Bool trialmode=(upmode.contains("TRIAL") ||
4564 0 : upmode.contains("FLAGONLY"));
4565 :
4566 : // Set up VisSet and its VisibilityIterator.
4567 :
4568 0 : VisibilityIterator::DataColumn whichOutCol = configureForCorrection ();
4569 :
4570 0 : VisIter& vi(vs_p->iter());
4571 0 : VisBufferAutoPtr vb (vi);
4572 0 : vi.origin();
4573 :
4574 : // Pass each timestamp (VisBuffer) to VisEquation for correction
4575 :
4576 0 : Vector<Bool> uncalspw(vi.numberSpw()); // Used to accumulate error messages
4577 0 : uncalspw.set(false); // instead of bombing the user
4578 0 : uncalspw.set(False); // instead of bombing the user
4579 : // in a loop.
4580 :
4581 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
4582 :
4583 0 : for (vi.origin(); vi.more(); vi++) {
4584 :
4585 0 : uInt spw = vb->spectralWindow();
4586 :
4587 : // Re-initialize weights from sigma column
4588 0 : vb->resetWeightMat();
4589 :
4590 : // If we can calibrate this vb, do it...
4591 0 : if (ve_p->spwOK(spw)){
4592 :
4593 : // throws exception if nothing to apply
4594 0 : ve_p->correct(*vb,trialmode);
4595 :
4596 : }
4597 : // ...else don't, prepare warning, and possibly set flags
4598 : else{
4599 :
4600 : // set uncalspw for warning message
4601 0 : uncalspw[spw] = true;
4602 : // set the flags, if we are being strict
4603 0 : if (upmode.contains("STRICT"))
4604 : // set the flags
4605 : // (don't touch the data/weights, which are initialized)
4606 0 : vb->flag().set(true);
4607 : }
4608 :
4609 : // Only if not a trial run, trigger write to disk
4610 0 : if (!upmode.contains("TRIAL")) {
4611 :
4612 0 : if (upmode.contains("CAL")) {
4613 0 : vi.setVis (vb->visCube(), whichOutCol);
4614 0 : vi.setWeightMat(vb->weightMat());
4615 : }
4616 :
4617 0 : if (upmode.contains("FLAG"))
4618 0 : vi.setFlag (vb->flag());
4619 :
4620 : }
4621 :
4622 : }
4623 : }
4624 :
4625 0 : vs_p->flush (); // Flush to disk
4626 :
4627 : // Now that we're out of the loop, summarize any errors.
4628 :
4629 0 : retval = summarize_uncalspws(uncalspw, "correct",
4630 0 : upmode.contains("STRICT"));
4631 :
4632 0 : actRec_=Record();
4633 0 : actRec_.define("origin","Calibrater::correct");
4634 0 : actRec_.defineRecord("VisEquation",ve_p->actionRec());
4635 :
4636 0 : }
4637 0 : catch (AipsError x) {
4638 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4639 0 : << LogIO::POST;
4640 :
4641 0 : logSink() << "Resetting all calibration application settings." << LogIO::POST;
4642 0 : unsetapply();
4643 :
4644 0 : throw(AipsError("Error in Calibrater::correct."));
4645 : retval = false; // Not that it ever gets here...
4646 0 : }
4647 0 : return retval;
4648 : }
4649 :
4650 0 : Bool OldCalibrater::corrupt() {
4651 :
4652 0 : logSink() << LogOrigin("Calibrater","corrupt") << LogIO::NORMAL;
4653 0 : Bool retval = true;
4654 :
4655 : try {
4656 :
4657 0 : if (!ok())
4658 0 : throw(AipsError("Calibrater not prepared for corrupt!"));
4659 :
4660 : // Nominally, we write out to the MODEL_DATA, unless absent
4661 0 : VisibilityIterator::DataColumn whichOutCol(VisibilityIterator::Model);
4662 :
4663 0 : if (!ms_p->tableDesc().isColumn("MODEL_DATA"))
4664 0 : throw(AipsError("MODEL_DATA column unexpectedly absent. Cannot corrupt."));
4665 :
4666 : // Ensure apply list non-zero and properly sorted
4667 0 : ve_p->setapply(vc_p);
4668 :
4669 : // Report the types that will be applied
4670 0 : applystate();
4671 :
4672 : // Arrange for iteration over data
4673 0 : Block<Int> columns;
4674 : // include scan iteration
4675 0 : columns.resize(5);
4676 0 : columns[0]=MS::ARRAY_ID;
4677 0 : columns[1]=MS::SCAN_NUMBER;
4678 0 : columns[2]=MS::FIELD_ID;
4679 0 : columns[3]=MS::DATA_DESC_ID;
4680 0 : columns[4]=MS::TIME;
4681 0 : vs_p->resetVisIter(columns,0.0);
4682 0 : VisIter& vi(vs_p->iter());
4683 0 : VisBuffer vb(vi);
4684 :
4685 : // Pass each timestamp (VisBuffer) to VisEquation for corruption.
4686 0 : Vector<Bool> uncalspw(vi.numberSpw()); // Used to accumulate error messages
4687 0 : uncalspw.set(false); // instead of bombing the user
4688 : // in a loop.
4689 0 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
4690 0 : Int spw = vi.spectralWindow();
4691 :
4692 : // Only proceed if spw can be calibrated
4693 0 : if (ve_p->spwOK(spw)) {
4694 :
4695 0 : for (vi.origin(); vi.more(); vi++) {
4696 :
4697 : // Corrupt the MODEL_DATA
4698 : // (note we are not treating weights and flags)
4699 0 : ve_p->corrupt(vb); // throws exception if nothing to apply
4700 0 : vi.setVis(vb.modelVisCube(),whichOutCol);
4701 :
4702 : }
4703 : }
4704 : else
4705 0 : uncalspw[spw] = true;
4706 : }
4707 : // Flush to disk
4708 0 : vs_p->flush();
4709 :
4710 : // Now that we're out of the loop, summarize any errors.
4711 0 : retval = summarize_uncalspws(uncalspw, "corrupt");
4712 0 : }
4713 0 : catch (AipsError x) {
4714 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
4715 0 : << LogIO::POST;
4716 :
4717 0 : logSink() << "Resetting all calibration application settings." << LogIO::POST;
4718 0 : unsetapply();
4719 :
4720 0 : throw(AipsError("Error in Calibrater::corrupt."));
4721 : retval = false; // Not that it ever gets here...
4722 0 : }
4723 0 : return retval;
4724 : }
4725 :
4726 0 : Bool OldCalibrater::initWeightsWithTsys(String wtmode, Bool dowtsp,
4727 : String tsystable, String gainfield, String interp, Vector<Int> spwmap) {
4728 :
4729 0 : logSink() << LogOrigin("Calibrater", "initWeightsWithTsys")
4730 0 : << LogIO::NORMAL;
4731 0 : Bool retval = true;
4732 :
4733 : try {
4734 :
4735 0 : if (!ok())
4736 0 : throw(AipsError("Calibrater not prepared for initWeights!"));
4737 :
4738 0 : String uptype = calTableType(tsystable);
4739 0 : if (!uptype.contains("TSYS")) {
4740 0 : throw(AipsError(
4741 0 : "Invalid calibration table type for Tsys weighting."));
4742 : }
4743 : // Set record format for calibration table application information
4744 0 : RecordDesc applyparDesc;
4745 0 : applyparDesc.addField("t", TpDouble);
4746 0 : applyparDesc.addField("table", TpString);
4747 0 : applyparDesc.addField("interp", TpString);
4748 0 : applyparDesc.addField("spw", TpArrayInt);
4749 0 : applyparDesc.addField("fieldstr", TpString);
4750 0 : applyparDesc.addField("calwt", TpBool);
4751 0 : applyparDesc.addField("spwmap", TpArrayInt);
4752 0 : applyparDesc.addField("opacity", TpArrayDouble);
4753 :
4754 : // Create record with the requisite field values
4755 0 : Record applypar(applyparDesc);
4756 0 : applypar.define("t", 0.0);
4757 0 : applypar.define("table", tsystable);
4758 0 : applypar.define("interp", interp);
4759 0 : applypar.define("spw", getSpwIdx(""));
4760 0 : applypar.define("fieldstr", gainfield);
4761 0 : applypar.define("calwt", true);
4762 0 : applypar.define("spwmap", spwmap);
4763 0 : applypar.define("opacity", Vector<Double>(1, 0.0));
4764 :
4765 0 : if (vc_p.nelements() > 0) {
4766 0 : logSink() << LogIO::WARN << "Resetting all calibration application settings." << LogIO::POST;
4767 0 : unsetapply();
4768 : }
4769 0 : logSink() << LogIO::NORMAL << "Weight initialization does not support selection. Resetting MS selection." << LogIO::POST;
4770 0 : selectvis();
4771 0 : StandardTsys vc = StandardTsys(*vs_p);
4772 0 : vc.setApply(applypar);
4773 :
4774 0 : logSink() << LogIO::NORMAL << ". " << vc.applyinfo() << LogIO::POST;
4775 0 : PtrBlock<VisCal*> vcb(1, &vc);
4776 : // Maintain sort of apply list
4777 0 : ve_p->setapply(vcb);
4778 :
4779 : // Detect WEIGHT_SPECTRUM and SIGMA_SPECTRUM
4780 0 : TableDesc mstd = ms_p->actualTableDesc();
4781 0 : String colWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
4782 0 : Bool wtspexists = mstd.isColumn(colWtSp);
4783 0 : String colSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
4784 0 : Bool sigspexists = mstd.isColumn(colSigSp);
4785 0 : Bool addsigsp = (dowtsp && !sigspexists);
4786 :
4787 : // Some log info
4788 0 : bool use_exposure = false;
4789 0 : if (wtmode == "tsys") {
4790 : logSink()
4791 : << "Initializing SIGMA and WEIGHT according to channel bandwidth and Tsys. NOTE this is an expert mode."
4792 0 : << LogIO::WARN << LogIO::POST;
4793 0 : } else if (wtmode == "tinttsys") {
4794 : logSink()
4795 : << "Initializing SIGMA and WEIGHT according to channel bandwidth, integration time, and Tsys. NOTE this is an expert mode."
4796 0 : << LogIO::WARN << LogIO::POST;
4797 0 : use_exposure = true;
4798 : } else {
4799 0 : throw(AipsError("Unrecognized wtmode specified: " + wtmode));
4800 : }
4801 :
4802 : // Force dowtsp if the column already exists
4803 0 : if (wtspexists && !dowtsp) {
4804 : logSink() << "Found WEIGHT_SPECTRUM; will force its initialization."
4805 0 : << LogIO::POST;
4806 0 : dowtsp = true;
4807 : }
4808 :
4809 : // Report that we are initializing the WEIGHT_SPECTRUM, and prepare to do so.
4810 0 : if (dowtsp) {
4811 :
4812 : // Ensure WEIGHT_SPECTRUM really exists at all
4813 : // (often it exists but is empty)
4814 0 : if (!wtspexists) {
4815 0 : logSink() << "Creating WEIGHT_SPECTRUM." << LogIO::POST;
4816 :
4817 : // Nominal defaulttileshape
4818 0 : IPosition dts(3, 4, 32, 1024);
4819 :
4820 : // Discern DATA's default tile shape and use it
4821 0 : const Record dminfo = ms_p->dataManagerInfo();
4822 0 : for (uInt i = 0; i < dminfo.nfields(); ++i) {
4823 0 : Record col = dminfo.asRecord(i);
4824 : //if (upcase(col.asString("NAME"))=="TILEDDATA") {
4825 0 : if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
4826 0 : dts = IPosition(
4827 0 : col.asRecord("SPEC").asArrayInt(
4828 0 : "DEFAULTTILESHAPE"));
4829 : //cout << "Found DATA's default tile: " << dts << endl;
4830 0 : break;
4831 : }
4832 0 : }
4833 :
4834 : // Add the column
4835 0 : String colWtSp = MS::columnName(MS::WEIGHT_SPECTRUM);
4836 0 : TableDesc tdWtSp;
4837 0 : tdWtSp.addColumn(
4838 0 : ArrayColumnDesc<Float>(colWtSp, "weight spectrum", 2));
4839 0 : TiledShapeStMan wtSpStMan("TiledWgtSpectrum", dts);
4840 0 : ms_p->addColumn(tdWtSp, wtSpStMan);
4841 0 : } else
4842 0 : logSink() << "Found WEIGHT_SPECTRUM." << LogIO::POST;
4843 : // Ensure WEIGHT_SPECTRUM really exists at all
4844 : // (often it exists but is empty)
4845 0 : if (!sigspexists) {
4846 0 : logSink() << "Creating SIGMA_SPECTRUM." << LogIO::POST;
4847 :
4848 : // Nominal defaulttileshape
4849 0 : IPosition dts(3, 4, 32, 1024);
4850 :
4851 : // Discern DATA's default tile shape and use it
4852 0 : const Record dminfo = ms_p->dataManagerInfo();
4853 0 : for (uInt i = 0; i < dminfo.nfields(); ++i) {
4854 0 : Record col = dminfo.asRecord(i);
4855 : //if (upcase(col.asString("NAME"))=="TILEDDATA") {
4856 0 : if (anyEQ(col.asArrayString("COLUMNS"), String("DATA"))) {
4857 0 : dts = IPosition(
4858 0 : col.asRecord("SPEC").asArrayInt(
4859 0 : "DEFAULTTILESHAPE"));
4860 : //cout << "Found DATA's default tile: " << dts << endl;
4861 0 : break;
4862 : }
4863 0 : }
4864 :
4865 : // Add the column
4866 0 : String colSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
4867 0 : TableDesc tdSigSp;
4868 0 : tdSigSp.addColumn(
4869 0 : ArrayColumnDesc<Float>(colSigSp, "sigma spectrum", 2));
4870 0 : TiledShapeStMan sigSpStMan("TiledSigtSpectrum", dts);
4871 0 : ms_p->addColumn(tdSigSp, sigSpStMan);
4872 : {
4873 0 : TableDesc loctd = ms_p->actualTableDesc();
4874 0 : String loccolSigSp = MS::columnName(MS::SIGMA_SPECTRUM);
4875 0 : AlwaysAssert(loctd.isColumn(loccolSigSp),AipsError);
4876 0 : }
4877 0 : }
4878 : }
4879 : else {
4880 0 : if (sigspexists) {
4881 0 : logSink() << "Removing SIGMA_SPECTRUM for non-channelized weight." << LogIO::POST;
4882 : if (true || ms_p->canRemoveColumn(colSigSp)) {
4883 0 : ms_p->removeColumn(colSigSp);
4884 : }
4885 : else
4886 : logSink() << LogIO::WARN << "Failed to remove SIGMA_SPECTRUM column. Values in SIGMA and SIGMA_SPECTRUM columns may be inconsistent after the operation." << LogIO::POST;
4887 : }
4888 : }
4889 :
4890 : // Arrange for iteration over data
4891 : // TBD: Be sure this sort is optimal for creating WS?
4892 0 : Block<Int> columns;
4893 : // include scan iteration
4894 0 : columns.resize(5);
4895 0 : columns[0] = MS::ARRAY_ID;
4896 0 : columns[1] = MS::SCAN_NUMBER;
4897 0 : columns[2] = MS::FIELD_ID;
4898 0 : columns[3] = MS::DATA_DESC_ID;
4899 0 : columns[4] = MS::TIME;
4900 :
4901 0 : vi::SortColumns sc(columns);
4902 0 : vi::VisibilityIterator2 vi2(*ms_p, sc, true);
4903 0 : vi::VisBuffer2 *vb = vi2.getVisBuffer();
4904 :
4905 0 : MSColumns mscol(*ms_p);
4906 0 : const MSSpWindowColumns& msspw(mscol.spectralWindow());
4907 0 : uInt nSpw = msspw.nrow();
4908 0 : Vector<Double> effChBw(nSpw, 0.0);
4909 0 : for (uInt ispw = 0; ispw < nSpw; ++ispw) {
4910 0 : effChBw[ispw] = mean(msspw.effectiveBW()(ispw));
4911 : }
4912 :
4913 0 : Int ivb(0);
4914 0 : for (vi2.originChunks(); vi2.moreChunks(); vi2.nextChunk()) {
4915 :
4916 0 : for (vi2.origin(); vi2.more(); vi2.next(), ++ivb) {
4917 :
4918 0 : Int spw = vb->spectralWindows()(0);
4919 :
4920 0 : auto nrow = vb->nRows();
4921 0 : Int nchan = vb->nChannels();
4922 0 : Int ncor = vb->nCorrelations();
4923 :
4924 : // Prepare for WEIGHT_SPECTRUM and SIGMA_SPECTRUM, if nec.
4925 0 : Cube<Float> newwtsp(0, 0, 0), newsigsp(0, 0, 0);
4926 0 : if (dowtsp) {
4927 0 : newwtsp.resize(ncor, nchan, nrow);
4928 0 : newwtsp.set(1.0);
4929 0 : newsigsp.resize(ncor, nchan, nrow);
4930 0 : newsigsp.set(1.0);
4931 : }
4932 :
4933 0 : if (ve_p->spwOK(spw)) {
4934 :
4935 : // Re-initialize weight info from sigma info
4936 : // This is smart wrt spectral weights, etc.
4937 : // (this makes W and WS, if present, "dirty" in the vb)
4938 : // TBD: only do this if !trial (else: avoid the I/O)
4939 : // vb->resetWeightsUsingSigma();
4940 : // Handle non-trivial modes
4941 : // Init WEIGHT, SIGMA from bandwidth & time
4942 0 : Matrix<Float> newwt(ncor, nrow), newsig(ncor, nrow);
4943 0 : newwt.set(1.0);
4944 0 : newsig.set(1.0);
4945 :
4946 : // Detect ACs
4947 0 : const Vector<Int> a1(vb->antenna1());
4948 0 : const Vector<Int> a2(vb->antenna2());
4949 0 : Vector<Bool> ac(a1 == a2);
4950 :
4951 : // XCs need an extra factor of 2
4952 0 : Vector<Float> xcfactor(nrow, 2.0);
4953 0 : xcfactor(ac) = 1.0; // (but not ACs)
4954 :
4955 : // The row-wise integration time
4956 0 : Vector<Float> expo(nrow);
4957 0 : convertArray(expo, vb->exposure());
4958 :
4959 : // Set weights to channel bandwidth first.
4960 0 : newwt.set(Float(effChBw(spw)));
4961 :
4962 : // For each correlation, apply exposure and xcfactor
4963 0 : for (Int icor = 0; icor < ncor; ++icor) {
4964 :
4965 0 : Vector<Float> wt(newwt.row(icor));
4966 0 : if (use_exposure) {
4967 0 : wt *= expo;
4968 : }
4969 0 : wt *= xcfactor;
4970 0 : if (dowtsp) {
4971 0 : for (Int ich = 0; ich < nchan; ++ich) {
4972 : Vector<Float> wtspi(
4973 0 : newwtsp(Slice(icor, 1, 1),
4974 0 : Slice(ich, 1, 1), Slice()).nonDegenerate(
4975 0 : IPosition(1, 2)));
4976 0 : wtspi = wt;
4977 0 : }
4978 : }
4979 0 : }
4980 : // Handle SIGMA_SPECTRUM
4981 0 : if (dowtsp) {
4982 0 : newsigsp = 1.0f / sqrt(newwtsp);
4983 : }
4984 : // sig from wt is inverse sqrt
4985 0 : newsig = 1.0f / sqrt(newwt);
4986 :
4987 : // Arrange write-back of both SIGMA and WEIGHT
4988 0 : vb->setSigma(newsig);
4989 0 : vb->setWeight(newwt);
4990 0 : if (dowtsp) {
4991 0 : vb->initWeightSpectrum(newwtsp);
4992 0 : vb->initSigmaSpectrum(newsigsp);
4993 : }
4994 : // Force writeback to disk (need to initialize weight/sigma before applying cal table)
4995 0 : vb->writeChangesBack();
4996 :
4997 : // Arrange for _in-place_ apply on CORRECTED_DATA (init from DATA)
4998 : // (this makes CD "dirty" in the vb)
4999 : // TBD: only do this if !trial (else: avoid the I/O)
5000 0 : vb->setVisCubeCorrected(vb->visCube());
5001 :
5002 : // Make flagcube dirty in the vb
5003 : // NB: we must _always_ do this I/O (even trial mode)
5004 0 : vb->setFlagCube(vb->flagCube());
5005 :
5006 : // Make all vb "not dirty"; we'll carefully arrange the writeback below
5007 0 : vb->dirtyComponentsClear();
5008 :
5009 : // throws exception if nothing to apply
5010 0 : ve_p->correct2(*vb, false, dowtsp);
5011 :
5012 0 : if (dowtsp) {
5013 0 : vb->setWeightSpectrum(vb->weightSpectrum());
5014 : // If WS was calibrated, set W to its channel-axis median
5015 0 : vb->setWeight( partialMedians(vb->weightSpectrum(), IPosition(1, 1)) );
5016 0 : newsigsp = 1.0f / sqrt(vb->weightSpectrum());
5017 0 : vb->initSigmaSpectrum(newsigsp);
5018 0 : vb->setSigma( partialMedians(newsigsp, IPosition(1, 1)) );
5019 : } else {
5020 0 : vb->setWeight(vb->weight());
5021 0 : newsig = 1.0f / sqrt(vb->weight());
5022 0 : vb->setSigma(newsig);
5023 : }
5024 : // Force writeback to disk
5025 0 : vb->writeChangesBack();
5026 :
5027 0 : } else {//Not calibrating the spw
5028 0 : if (dowtsp && !wtspexists) {
5029 : // newly created WS Need to initialize
5030 0 : vb->initWeightSpectrum(newwtsp);
5031 : }
5032 0 : if (addsigsp) {
5033 : // newly created SS Need to initialize
5034 0 : vb->initSigmaSpectrum(newsigsp);
5035 0 : vb->writeChangesBack();
5036 : }
5037 : }
5038 0 : }
5039 : }
5040 : // clear-up Tsys caltable from list of apply
5041 0 : unsetapply();
5042 :
5043 0 : } catch (AipsError x) {
5044 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg()
5045 0 : << LogIO::POST;
5046 :
5047 : logSink() << "Resetting all calibration application settings."
5048 0 : << LogIO::POST;
5049 0 : unsetapply();
5050 :
5051 0 : throw(AipsError("Error in Calibrater::initWeights."));
5052 : retval = false; // Not that it ever gets here...
5053 0 : }
5054 0 : return retval;
5055 : }
5056 :
5057 0 : Bool OldCalibrater::solve() {
5058 :
5059 0 : logSink() << LogOrigin("Calibrater","solve") << LogIO::NORMAL3;
5060 :
5061 : try {
5062 :
5063 0 : if (!ok())
5064 0 : throw(AipsError("Calibrater not prepared for solve."));
5065 :
5066 : // Handle nothing-to-solve-for case
5067 0 : if (!svc_p)
5068 0 : throw(AipsError("Please run setsolve before attempting to solve."));
5069 :
5070 : // Handle specified caltable
5071 0 : if (true && svc_p) {
5072 :
5073 : /*
5074 : cout << "name: " << svc_p->calTableName() << endl;
5075 : cout << boolalpha;
5076 : cout << "append? " << svc_p->append() << endl;
5077 : cout << "opened? " << Table::isOpened(svc_p->calTableName()) << endl;
5078 : cout << "readable? " << Table::isReadable(svc_p->calTableName()) << endl;
5079 : // cout << "writable? " << Table::isWritable(svc_p->calTableName()) << endl;
5080 : cout << "canDelete? " << Table::canDeleteTable(svc_p->calTableName(),true) << endl;
5081 : */
5082 :
5083 : // If table exists (readable) and not deletable
5084 : // we have to abort (append=T requires deletable)
5085 0 : if ( Table::isReadable(svc_p->calTableName()) &&
5086 0 : !TableUtil::canDeleteTable(svc_p->calTableName()) ) {
5087 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?)."));
5088 : }
5089 : }
5090 :
5091 : // Arrange VisEquation for solve
5092 0 : ve_p->setsolve(*svc_p);
5093 :
5094 : // Ensure apply list properly sorted w.r.t. solvable term
5095 0 : ve_p->setapply(vc_p);
5096 :
5097 : // Report what is being applied and solved-for
5098 0 : applystate();
5099 0 : solvestate();
5100 :
5101 :
5102 : // Report correct/corrupt apply order
5103 : // ve_p->state();
5104 :
5105 : // Set the channel mask
5106 0 : svc_p->setChanMask(chanmask_);
5107 :
5108 : // Generally use standard solver
5109 0 : if (svc_p->useGenericGatherForSolve())
5110 0 : genericGatherAndSolve(); // using VisBuffGroupAcc
5111 : else {
5112 : //cout << "Fully self-directed data gather and solve" << endl;
5113 : // Fully self-directed data gather and solve
5114 0 : svc_p->selfGatherAndSolve(*vs_p,*ve_p);
5115 : }
5116 :
5117 0 : svc_p->clearChanMask();
5118 :
5119 0 : } catch (AipsError x) {
5120 0 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
5121 :
5122 0 : logSink() << "Reseting entire solve/apply state." << LogIO::POST;
5123 0 : reset();
5124 :
5125 0 : throw(AipsError("Error in Calibrater::solve."));
5126 : return false;
5127 0 : }
5128 :
5129 0 : return true;
5130 :
5131 : }
5132 :
5133 0 : Vector<Double> OldCalibrater::modelfit(const Int& niter,
5134 : const String& stype,
5135 : const Vector<Double>& par,
5136 : const Vector<Bool>& vary,
5137 : const String& file) {
5138 :
5139 : /*
5140 : cout << "Calibrater::modelfit" << endl;
5141 : cout << " niter = " << niter << endl;
5142 : cout << " stype = " << stype << endl;
5143 : cout << " par = " << par << endl;
5144 : cout << " vary = " << vary << endl;
5145 : cout << " file = " << file << endl;
5146 : */
5147 : // logSink() << LogOrigin("Calibrater","modelfit") << LogIO::NORMAL;
5148 :
5149 : try {
5150 0 : if(!ok()) throw(AipsError("Calibrater not ok()"));
5151 :
5152 : // Construct UVMod with the VisSet
5153 0 : UVMod uvmod(*vs_p);
5154 :
5155 0 : if (stype=="P")
5156 0 : uvmod.setModel(ComponentType::POINT, par, vary);
5157 0 : else if (stype=="G")
5158 0 : uvmod.setModel(ComponentType::GAUSSIAN, par, vary);
5159 0 : else if (stype=="D")
5160 0 : uvmod.setModel(ComponentType::DISK, par, vary);
5161 : else
5162 0 : throw(AipsError("Unrecognized component type in Calibrater::modelfit."));
5163 :
5164 : // Run the fit
5165 0 : uvmod.modelfit(niter,file);
5166 :
5167 : // Return the parameter vector
5168 0 : return uvmod.par();
5169 :
5170 0 : } catch (AipsError x) {
5171 0 : logSink() << LogIO::SEVERE << "Caught exception: " << x.getMesg() << LogIO::POST;
5172 0 : throw(AipsError("Error in Calibrater::modelfit."));
5173 :
5174 : return Vector<Double>();
5175 0 : }
5176 :
5177 : }
5178 :
5179 0 : void OldCalibrater::fluxscale(const String& infile,
5180 : const String& outfile,
5181 : const Vector<Int>& refField,
5182 : const Vector<Int>& refSpwMap,
5183 : const Vector<Int>& tranField,
5184 : const Bool& append,
5185 : const Float& inGainThres,
5186 : const String& antSel,
5187 : const String& timerangeSel,
5188 : const String& scanSel,
5189 : SolvableVisCal::fluxScaleStruct& oFluxScaleFactor,
5190 : const String& oListFile,
5191 : const Bool& incremental,
5192 : const Int& fitorder,
5193 : const Bool& display) {
5194 :
5195 : // throw(AipsError("Method 'fluxscale' is temporarily disabled."));
5196 :
5197 : // TBD: write inputs to MSHistory
5198 0 : logSink() << LogOrigin("Calibrater","fluxscale") << LogIO::NORMAL3;
5199 :
5200 0 : SolvableVisCal *fsvj_(NULL);
5201 : try {
5202 : // If infile is Calibration table
5203 0 : if (Table::isReadable(infile) &&
5204 0 : TableUtil::tableInfo(infile).type()=="Calibration") {
5205 :
5206 : // get calibration type
5207 0 : String caltype;
5208 0 : caltype = TableUtil::tableInfo(infile).subType();
5209 : logSink() << "Table " << infile
5210 : << " is of type: "<< caltype
5211 0 : << LogIO::POST;
5212 0 : String message="Table "+infile+" is of type: "+caltype;
5213 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
5214 :
5215 : // form selection
5216 0 : String select="";
5217 : // Selection is empty for case of no tran specification
5218 0 : if (tranField.nelements()>0) {
5219 :
5220 : // All selected fields
5221 0 : Vector<Int> allflds = concatenateArray(refField,tranField);
5222 :
5223 : // Assemble TaQL
5224 0 : ostringstream selectstr;
5225 0 : selectstr << "FIELD_ID IN [";
5226 0 : for (Int iFld=0; iFld<allflds.shape(); iFld++) {
5227 0 : if (iFld>0) selectstr << ", ";
5228 0 : selectstr << allflds(iFld);
5229 : }
5230 0 : selectstr << "]";
5231 0 : select=selectstr.str();
5232 0 : }
5233 :
5234 : // Construct proper SVC object
5235 0 : if (caltype == "G Jones") {
5236 0 : fsvj_ = createSolvableVisCal("G",*vs_p);
5237 0 : } else if (caltype == "T Jones") {
5238 0 : fsvj_ = createSolvableVisCal("T",*vs_p);
5239 : } else {
5240 : // Can't process other than G and T (add B?)
5241 0 : ostringstream typeErr;
5242 0 : typeErr << "Type " << caltype
5243 0 : << " not supported in fluxscale.";
5244 :
5245 0 : throw(AipsError(typeErr.str()));
5246 0 : }
5247 :
5248 : // fill table with selection
5249 0 : RecordDesc applyparDesc;
5250 0 : applyparDesc.addField ("table", TpString);
5251 0 : applyparDesc.addField ("select", TpString);
5252 0 : Record applypar(applyparDesc);
5253 0 : applypar.define ("table", infile);
5254 0 : applypar.define ("select", select);
5255 0 : fsvj_->setApply(applypar);
5256 :
5257 : //Bool incremental=false;
5258 : // Make fluxscale calculation
5259 0 : Vector<String> fldnames(MSFieldColumns(ms_p->field()).name().getColumn());
5260 : //fsvj_->fluxscale(refField,tranField,refSpwMap,fldnames,oFluxScaleFactor,
5261 0 : fsvj_->fluxscale(outfile,refField,tranField,refSpwMap,fldnames,inGainThres,antSel,
5262 : timerangeSel,scanSel,oFluxScaleFactor, oListFile,incremental,fitorder,display);
5263 : // oListFile);
5264 :
5265 : // If no outfile specified, use infile (overwrite!)
5266 0 : String out(outfile);
5267 0 : if (out.length()==0)
5268 0 : out = infile;
5269 :
5270 : // Store result
5271 0 : if (append) {
5272 0 : logSink() << "Appending result to " << out << LogIO::POST;
5273 0 : String message="Appending result to "+out;
5274 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
5275 0 : } else {
5276 0 : logSink() << "Storing result in " << out << LogIO::POST;
5277 0 : String message="Storing result in "+out;
5278 0 : MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
5279 0 : }
5280 0 : fsvj_->storeNCT(out,append);
5281 :
5282 : // Clean up
5283 0 : delete fsvj_;
5284 :
5285 0 : } else {
5286 : // Table not found/unreadable, or not Calibration
5287 0 : ostringstream tabErr;
5288 0 : tabErr << "File " << infile
5289 0 : << " does not exist or is not a Calibration Table.";
5290 :
5291 0 : throw(AipsError(tabErr.str()));
5292 :
5293 0 : }
5294 0 : } catch (AipsError x) {
5295 :
5296 : logSink() << LogIO::SEVERE
5297 : << "Caught Exception: "
5298 : << x.getMesg()
5299 0 : << LogIO::POST;
5300 :
5301 : // Clean up
5302 0 : if (fsvj_) delete fsvj_;
5303 :
5304 : // Write to MS History table
5305 : // String message="Caught Exception: "+x.getMesg();
5306 : // MSHistoryHandler::addMessage(*ms_p, message, "calibrater", "", "calibrater::fluxscale()");
5307 :
5308 0 : throw(AipsError("Error in Calibrater::fluxscale."));
5309 :
5310 : return;
5311 :
5312 0 : }
5313 0 : return;
5314 :
5315 :
5316 : }
5317 :
5318 0 : void OldCalibrater::accumulate(const String& intab,
5319 : const String& incrtab,
5320 : const String& outtab,
5321 : const String& fields,
5322 : const String& calFields,
5323 : const String& interp,
5324 : const Double& t,
5325 : const Vector<Int>& spwmap) {
5326 :
5327 : // logSink() << LogOrigin("Calibrater","accumulate") << LogIO::NORMAL;
5328 :
5329 0 : logSink() << "Beginning accumulate." << LogIO::POST;
5330 :
5331 : // SVJ objects:
5332 0 : SolvableVisCal *incal_(NULL), *incrcal_(NULL);
5333 :
5334 : try {
5335 :
5336 : /*
5337 : cout << "intab = " << intab << endl;
5338 : cout << "incrtab = " << incrtab << endl;
5339 : cout << "outtab = " << outtab << endl;
5340 : cout << "fields = " << fields << endl;
5341 : cout << "calFields = " << calFields << endl;
5342 : cout << "interp = " << interp << endl;
5343 : cout << "t = " << t << endl;
5344 : */
5345 :
5346 : // Incremental table's type sets the type we are dealing with
5347 0 : String caltype=calTableType(incrtab);
5348 :
5349 : // If no input cumulative timescale specified, then
5350 : // a valid input cumulative table must be specified
5351 0 : if (t < 0.0) {
5352 :
5353 0 : String intype=calTableType(intab);
5354 :
5355 0 : if (intype!=caltype) {
5356 :
5357 0 : ostringstream typeErr;
5358 0 : typeErr << "Table " << intab
5359 0 : << " is not the same type as "
5360 0 : << incrtab << " (" << caltype << ")";
5361 0 : throw(AipsError(typeErr.str()));
5362 0 : }
5363 0 : }
5364 :
5365 : // At this point all is ok; we will:
5366 : // o fill from intab and accumulate to it (t<0), OR
5367 : // o create a new cumulative table from scratch (t>0)
5368 :
5369 : // If creating a new cumulative table, it must span the whole dataset,
5370 : // so reset data selection to whole MS, and setup iterator
5371 0 : if (t>0.0) {
5372 0 : selectvis();
5373 0 : Block<Int> columns;
5374 0 : columns.resize(4);
5375 0 : columns[0]=MS::ARRAY_ID;
5376 0 : columns[1]=MS::TIME;
5377 0 : columns[2]=MS::FIELD_ID;
5378 0 : columns[3]=MS::DATA_DESC_ID;
5379 0 : vs_p->resetVisIter(columns,t);
5380 0 : }
5381 :
5382 : // logSink() << "Table " << infile
5383 : // << " is of type: "<< caltype
5384 : // << LogIO::POST;
5385 :
5386 0 : incal_ = createSolvableVisCal(caltype,*vs_p);
5387 0 : incrcal_ = createSolvableVisCal(caltype,*vs_p);
5388 :
5389 : // TBD: move to svj.setAccumulate?
5390 0 : if ( !(incal_->accumulatable()) ) {
5391 0 : ostringstream typeErr;
5392 0 : typeErr << "Type " << caltype
5393 0 : << " not yet supported in accumulate.";
5394 0 : throw(AipsError(typeErr.str()));
5395 0 : }
5396 :
5397 : // At this point, accumulation apparently supported,
5398 : // so continue...
5399 :
5400 : // initialize the cumulative solutions
5401 0 : incal_->setAccumulate(*vs_p,intab,"",t,-1);
5402 :
5403 :
5404 : /*
5405 : // form selection on incr table
5406 : String incrSel="";
5407 : if (calFields.shape()>0) {
5408 :
5409 : // Assemble TaQL
5410 : ostringstream selectstr;
5411 : selectstr << "FIELD_ID IN [";
5412 : for (Int iFld=0; iFld<calFields.shape(); iFld++) {
5413 : if (iFld>0) selectstr << ", ";
5414 : selectstr << calFields(iFld);
5415 : }
5416 : selectstr << "]";
5417 : incrSel=selectstr.str();
5418 : }
5419 : */
5420 :
5421 : // fill incr table with selection
5422 : logSink() << "Preparing to accumulate calibration from table: "
5423 : << incrtab
5424 0 : << LogIO::POST;
5425 :
5426 : // Set record format for calibration table application information
5427 0 : RecordDesc applyparDesc;
5428 0 : applyparDesc.addField ("t", TpDouble);
5429 0 : applyparDesc.addField ("table", TpString);
5430 : // applyparDesc.addField ("select", TpString);
5431 0 : applyparDesc.addField ("fieldstr", TpString);
5432 0 : applyparDesc.addField ("interp", TpString);
5433 0 : applyparDesc.addField ("spwmap",TpArrayInt);
5434 :
5435 : // Create record with the requisite field values
5436 0 : Record applypar(applyparDesc);
5437 0 : applypar.define ("t", t);
5438 0 : applypar.define ("table", incrtab);
5439 : // applypar.define ("select", incrSel);
5440 0 : applypar.define ("fieldstr", calFields);
5441 0 : applypar.define ("interp", interp);
5442 0 : applypar.define ("spwmap",spwmap);
5443 :
5444 0 : incrcal_->setApply(applypar);
5445 :
5446 0 : Vector<Int> fldidx(0);
5447 0 : if (fields.length()>0)
5448 0 : fldidx=getFieldIdx(fields);
5449 :
5450 : // All ready, now do the accumulation
5451 0 : incal_->accumulate(incrcal_,fldidx);
5452 :
5453 : // ...and store the result
5454 : logSink() << "Storing accumulated calibration in table: "
5455 : << outtab
5456 0 : << LogIO::POST;
5457 :
5458 0 : if (outtab != "")
5459 0 : incal_->calTableName()=outtab;
5460 :
5461 0 : incal_->storeNCT();
5462 :
5463 0 : delete incal_;
5464 0 : delete incrcal_;
5465 :
5466 : logSink() << "Finished accumulation."
5467 0 : << LogIO::POST;
5468 :
5469 0 : } catch (AipsError x) {
5470 : logSink() << LogIO::SEVERE
5471 : << "Caught Exception: "
5472 : << x.getMesg()
5473 0 : << LogIO::POST;
5474 :
5475 0 : if (incal_) delete incal_;
5476 0 : if (incrcal_) delete incrcal_;
5477 :
5478 0 : throw(AipsError("Error in Calibrater::accumulate."));
5479 : return;
5480 0 : }
5481 0 : return;
5482 :
5483 : }
5484 :
5485 0 : void OldCalibrater::specifycal(const String& type,
5486 : const String& caltable,
5487 : const String& time,
5488 : const String& spw,
5489 : const String& antenna,
5490 : const String& pol,
5491 : const Vector<Double>& parameter,
5492 : const String& infile,
5493 : const Bool& uniform) {
5494 :
5495 0 : logSink() << LogOrigin("Calibrater","specifycal") << LogIO::NORMAL;
5496 :
5497 : // SVJ objects:
5498 0 : SolvableVisCal *cal_(NULL);
5499 :
5500 : try {
5501 :
5502 : // Set record format for calibration table application information
5503 0 : RecordDesc specifyDesc;
5504 0 : specifyDesc.addField ("caltable", TpString);
5505 0 : specifyDesc.addField ("time", TpString);
5506 0 : specifyDesc.addField ("spw", TpArrayInt);
5507 0 : specifyDesc.addField ("antenna", TpArrayInt);
5508 0 : specifyDesc.addField ("pol", TpString);
5509 0 : specifyDesc.addField ("parameter", TpArrayDouble);
5510 0 : specifyDesc.addField ("caltype",TpString);
5511 0 : specifyDesc.addField ("infile",TpString);
5512 0 : specifyDesc.addField ("uniform",TpBool);
5513 :
5514 : // Create record with the requisite field values
5515 0 : Record specify(specifyDesc);
5516 0 : specify.define ("caltable", caltable);
5517 0 : specify.define ("time", time);
5518 0 : if (spw=="*")
5519 0 : specify.define ("spw",Vector<Int>(1,-1));
5520 : else
5521 0 : specify.define ("spw",getSpwIdx(spw));
5522 0 : if (antenna=="*")
5523 0 : specify.define ("antenna",Vector<Int>(1,-1) );
5524 : else
5525 0 : specify.define ("antenna",getAntIdx(antenna));
5526 0 : specify.define ("pol",pol);
5527 0 : specify.define ("parameter",parameter);
5528 0 : specify.define ("caltype",type);
5529 0 : specify.define ("infile",infile);
5530 0 : specify.define ("uniform",uniform);
5531 :
5532 : // Now do it
5533 0 : String utype=upcase(type);
5534 0 : if (utype=="G" || utype.contains("AMP") || utype.contains("PH"))
5535 0 : cal_ = createSolvableVisCal("G",*vs_p);
5536 0 : else if (utype=='K' || utype.contains("SBD") || utype.contains("DELAY"))
5537 0 : cal_ = createSolvableVisCal("K",*vs_p);
5538 0 : else if (utype.contains("MBD"))
5539 0 : cal_ = createSolvableVisCal("K",*vs_p); // as of 5.3, KMBD is just K
5540 0 : else if (utype.contains("ANTPOS"))
5541 0 : cal_ = createSolvableVisCal("KANTPOS",*vs_p);
5542 0 : else if (utype.contains("TSYS"))
5543 0 : cal_ = createSolvableVisCal("TSYS",*vs_p);
5544 0 : else if (utype.contains("EVLAGAIN") ||
5545 0 : utype.contains("SWP") ||
5546 0 : utype.contains("RQ"))
5547 0 : cal_ = createSolvableVisCal("EVLASWP",*vs_p);
5548 0 : else if (utype.contains("OPAC"))
5549 0 : cal_ = createSolvableVisCal("TOPAC",*vs_p);
5550 0 : else if (utype.contains("GC") || utype.contains("EFF"))
5551 0 : cal_ = createSolvableVisCal("GAINCURVE",*vs_p);
5552 0 : else if (utype.contains("TEC"))
5553 0 : cal_ = createSolvableVisCal("TEC",*vs_p);
5554 0 : else if (utype.contains("SDSKY_PS"))
5555 0 : cal_ = createSolvableVisCal("SDSKY_PS",*vs_p);
5556 0 : else if (utype.contains("SDSKY_RASTER"))
5557 0 : cal_ = createSolvableVisCal("SDSKY_RASTER",*vs_p);
5558 0 : else if (utype.contains("SDSKY_OTF"))
5559 0 : cal_ = createSolvableVisCal("SDSKY_OTF",*vs_p);
5560 : else
5561 0 : throw(AipsError("Unrecognized caltype."));
5562 :
5563 : // set up for specification (set up the CalSet)
5564 0 : cal_->setSpecify(specify);
5565 :
5566 : // fill with specified values
5567 0 : cal_->specify(specify);
5568 :
5569 : // Store result
5570 0 : cal_->storeNCT();
5571 :
5572 0 : delete cal_;
5573 :
5574 0 : } catch (AipsError x) {
5575 : logSink() << LogIO::SEVERE
5576 : << "Caught Exception: "
5577 : << x.getMesg()
5578 0 : << LogIO::POST;
5579 :
5580 0 : if (cal_) delete cal_;
5581 :
5582 0 : throw(AipsError("Error in Calibrater::specifycal."));
5583 : return;
5584 0 : }
5585 0 : return;
5586 :
5587 : }
5588 :
5589 0 : Bool OldCalibrater::smooth(const String& infile,
5590 : String& outfile, // const Bool& append,
5591 : const String& smoothtype,
5592 : const Double& smoothtime,
5593 : const String& fields)
5594 : {
5595 :
5596 : // TBD: support append?
5597 : // TBD: spw selection?
5598 :
5599 0 : logSink() << LogOrigin("Calibrater","smooth") << LogIO::NORMAL;
5600 :
5601 0 : logSink() << "Beginning smoothing/interpolating method." << LogIO::POST;
5602 :
5603 :
5604 : // A pointer to an SVC
5605 0 : SolvableVisCal *svc(NULL);
5606 :
5607 : try {
5608 :
5609 : // Handle no in file
5610 0 : if (infile=="")
5611 0 : throw(AipsError("Please specify an input calibration table."));
5612 :
5613 : // Handle bad smoothtype
5614 0 : if (smoothtype!="mean" && smoothtype!="median")
5615 0 : throw(AipsError("Unrecognized smooth type!"));
5616 :
5617 : // Handle bad smoothtime
5618 0 : if (smoothtime<=0)
5619 0 : throw(AipsError("Please specify a strictly positive smoothtime."));
5620 :
5621 : // Handle no outfile
5622 0 : if (outfile=="") {
5623 0 : outfile=infile;
5624 : logSink() << "Will overwrite input file with smoothing result."
5625 0 : << LogIO::POST;
5626 : }
5627 :
5628 :
5629 0 : svc = createSolvableVisCal(calTableType(infile),*vs_p);
5630 :
5631 0 : if (svc->smoothable()) {
5632 :
5633 : // Fill calibration table using setApply
5634 0 : RecordDesc applyparDesc;
5635 0 : applyparDesc.addField ("table", TpString);
5636 0 : Record applypar(applyparDesc);
5637 0 : applypar.define ("table", infile);
5638 0 : svc->setApply(applypar);
5639 :
5640 : // Convert refFields/transFields to index lists
5641 0 : Vector<Int> fldidx(0);
5642 0 : if (fields.length()>0)
5643 0 : fldidx=getFieldIdx(fields);
5644 :
5645 : // Delegate to SVC
5646 0 : svc->smooth(fldidx,smoothtype,smoothtime);
5647 :
5648 : // Store the result on disk
5649 : // if (append) logSink() << "Appending result to " << outfile << LogIO::POST;
5650 : //else
5651 0 : logSink() << "Storing result in " << outfile << LogIO::POST;
5652 :
5653 :
5654 0 : if (outfile != "")
5655 0 : svc->calTableName()=outfile;
5656 0 : svc->storeNCT();
5657 :
5658 : // Clean up
5659 0 : if (svc) delete svc; svc=NULL;
5660 :
5661 : // Apparently, it worked
5662 0 : return true;
5663 :
5664 0 : }
5665 : else
5666 0 : throw(AipsError("This type ("+svc->typeName()+") does not support smoothing."));
5667 :
5668 0 : } catch (AipsError x) {
5669 :
5670 : logSink() << LogIO::SEVERE
5671 : << "Caught Exception: "
5672 : << x.getMesg()
5673 0 : << LogIO::POST;
5674 : // Clean up
5675 0 : if (svc) delete svc; svc=NULL;
5676 :
5677 0 : throw(AipsError("Error in Calibrater::smooth."));
5678 :
5679 : return false;
5680 0 : }
5681 : return false;
5682 : }
5683 :
5684 : // List a calibration table
5685 0 : Bool OldCalibrater::listCal(const String& infile,
5686 : const String& field,
5687 : const String& antenna,
5688 : const String& spw,
5689 : const String& listfile,
5690 : const Int& pagerows) {
5691 :
5692 0 : SolvableVisCal *svc(NULL);
5693 0 : logSink() << LogOrigin("Calibrater","listCal");
5694 :
5695 : try {
5696 :
5697 : // Trap (currently) unsupported types
5698 0 : if (upcase(calTableType(infile))=="GSPLINE" ||
5699 0 : upcase(calTableType(infile))=="BPOLY")
5700 0 : throw(AipsError("GSPLINE and BPOLY tables cannot currently be listed."));
5701 :
5702 : // Get user's selected fields, ants
5703 0 : Vector<Int> ufldids=getFieldIdx(field);
5704 0 : Vector<Int> uantids=getAntIdx(antenna);
5705 :
5706 0 : String newSpw = spw;
5707 0 : Bool defaultSelect = false;
5708 0 : if (spw.empty()) { // list all channels (default)
5709 0 : defaultSelect = true;
5710 0 : newSpw = "*";
5711 : logSink() << LogIO::NORMAL1 << "Spws selected: ALL" << endl
5712 0 : << "Channels selected: ALL" << LogIO::POST;
5713 : }
5714 : // Get user's selected spw and channels
5715 0 : Vector<Int> uspwids=getSpwIdx(newSpw);
5716 0 : Matrix<Int> uchanids=getChanIdx(newSpw);
5717 0 : if (!defaultSelect) {
5718 : logSink() << LogIO::NORMAL1 << "Spw and Channel selection matrix: "
5719 : << endl << "Each rows shows: [ Spw , Start Chan , Stop Chan , Chan Step ]"
5720 0 : << endl << uchanids << LogIO::POST;
5721 : }
5722 : logSink() << LogIO::DEBUG2
5723 : << "uspwids = " << uspwids << endl
5724 0 : << "uchanids = " << uchanids << LogIO::POST;
5725 :
5726 : // By default, do first spw, first chan
5727 0 : if (uspwids.nelements()==0) {
5728 0 : uchanids.resize(1,4);
5729 0 : uchanids=0;
5730 : }
5731 :
5732 : // Set record format for calibration table application information
5733 0 : RecordDesc applyparDesc;
5734 0 : applyparDesc.addField ("table", TpString);
5735 :
5736 : // Create record with the requisite field values
5737 0 : Record applypar(applyparDesc);
5738 0 : applypar.define ("table", infile);
5739 :
5740 : // Generate the VisCal to be listed
5741 0 : svc = createSolvableVisCal(calTableType(infile),*vs_p);
5742 0 : svc->setApply(applypar);
5743 :
5744 : // list it
5745 0 : svc->listCal(ufldids,uantids,uchanids, //uchanids(0,0),uchanids(0,1),
5746 : listfile,pagerows);
5747 :
5748 0 : if (svc) delete svc; svc=NULL;
5749 :
5750 0 : return true;
5751 :
5752 0 : } catch (AipsError x) {
5753 :
5754 : logSink() << LogIO::SEVERE
5755 : << "Caught Exception: "
5756 : << x.getMesg()
5757 0 : << LogIO::POST;
5758 : // Clean up
5759 0 : if (svc) delete svc; svc=NULL;
5760 :
5761 0 : throw(AipsError("Error in Calibrater::listCal."));
5762 :
5763 : return false;
5764 0 : }
5765 : return false;
5766 :
5767 : }
5768 :
5769 0 : Bool OldCalibrater::initialize(MeasurementSet& inputMS,
5770 : Bool compress,
5771 : Bool addScratch, Bool addModel) {
5772 :
5773 0 : logSink() << LogOrigin("Calibrater","") << LogIO::NORMAL3;
5774 :
5775 : try {
5776 0 : timer_p.mark();
5777 :
5778 : // Set pointer ms_p from input MeasurementSet
5779 0 : if (ms_p) {
5780 0 : *ms_p=inputMS;
5781 : } else {
5782 0 : ms_p = new MeasurementSet(inputMS);
5783 0 : AlwaysAssert(ms_p,AipsError);
5784 : };
5785 :
5786 : // Disabled 2016/04/19 (gmoellen): avoid direct MS.HISTORY
5787 : // updates from below the python level, FOR NOW
5788 :
5789 : /*
5790 :
5791 : // Setup to write LogIO to HISTORY Table in MS
5792 : if(!(Table::isReadable(ms_p->historyTableName()))){
5793 : // create a new HISTORY table if its not there
5794 : TableRecord &kws = ms_p->rwKeywordSet();
5795 : SetupNewTable historySetup(ms_p->historyTableName(),
5796 : MSHistory::requiredTableDesc(),Table::New);
5797 : kws.defineTable(MS::keywordName(MS::HISTORY), Table(historySetup));
5798 : MSHistoryHandler::addMessage(*ms_p, "HISTORY Table added by Calibrater",
5799 : "Calibrater","","Calibrater::initialize()");
5800 : }
5801 : historytab_p=Table(ms_p->historyTableName(),
5802 : TableLock(TableLock::UserNoReadLocking), Table::Update);
5803 : // jagonzal (CAS-4110): When the selectvis method throws an exception the initialize method
5804 : // is called again to leave the calibrater in a proper state, and since there was a previous
5805 : // initialization the history handler was already created, and has to be destroyed before
5806 : // creating a new one to avoid leaveing the HISTORY table opened.
5807 : if (hist_p) delete hist_p;
5808 : hist_p= new MSHistoryHandler(*ms_p, "calibrater");
5809 :
5810 : // (2016/04/19) */
5811 :
5812 : // Remember the ms's name
5813 0 : msname_p=ms_p->tableName();
5814 :
5815 :
5816 : // Add/init scr cols, if requested (init is hard-wired)
5817 0 : if (addScratch || addModel) {
5818 0 : Bool alsoinit=true;
5819 0 : VisSetUtil::addScrCols(*ms_p,addModel,addScratch,alsoinit,compress);
5820 : }
5821 :
5822 : // Set the selected MeasurementSet to be the same initially
5823 : // as the input MeasurementSet
5824 0 : if (mssel_p) delete mssel_p;
5825 0 : mssel_p=new MeasurementSet(*ms_p);
5826 :
5827 : logSink() << LogIO::NORMAL
5828 : << "Initializing nominal selection to the whole MS."
5829 0 : << LogIO::POST;
5830 :
5831 :
5832 : // Create a VisSet with no selection
5833 : // (gmoellen 2012/02/06: this merely makes a VisIter now)
5834 0 : if (vs_p) {
5835 0 : delete vs_p;
5836 0 : vs_p=0;
5837 : };
5838 0 : Block<Int> nosort(0);
5839 0 : Matrix<Int> noselection;
5840 0 : Double timeInterval=0;
5841 : // gmoellen 2012/02/06 vs_p=new VisSet(*ms_p,nosort,noselection,addScratch,timeInterval,compress, addModel);
5842 0 : vs_p=new VisSet(*ms_p,nosort,noselection,false,timeInterval,false,false);
5843 :
5844 : // Size-up the chanmask PB
5845 0 : initChanMask();
5846 :
5847 : // Create the associated VisEquation
5848 : // TBD: move to ctor and make it non-pointer
5849 0 : if (ve_p) {
5850 0 : delete ve_p;
5851 0 : ve_p=0;
5852 : };
5853 0 : ve_p=new VisEquation();
5854 :
5855 : // Reset the apply/solve VisCals
5856 0 : reset(true,true);
5857 :
5858 0 : return true;
5859 :
5860 0 : } catch (AipsError x) {
5861 0 : logSink() << LogOrigin("Calibrater","initialize",WHERE)
5862 : << LogIO::SEVERE << "Caught exception: " << x.getMesg()
5863 0 : << LogIO::POST;
5864 0 : cleanup();
5865 0 : cleanupVisSet();
5866 0 : if (ms_p) delete ms_p; ms_p=NULL;
5867 0 : if (hist_p) delete hist_p; hist_p=NULL;
5868 :
5869 0 : throw(AipsError("Error in Calibrater::initialize()"));
5870 : return false;
5871 0 : }
5872 : return false;
5873 : }
5874 :
5875 0 : Bool OldCalibrater::initCalSet(const Int& calSet)
5876 : {
5877 :
5878 : // logSink() << LogOrigin("Calibrater","initCalSet") << LogIO::NORMAL3;
5879 :
5880 0 : if (vs_p) {
5881 :
5882 0 : Block<Int> columns;
5883 : // include scan iteration, for more optimal iteration
5884 0 : columns.resize(5);
5885 0 : columns[0]=MS::ARRAY_ID;
5886 0 : columns[1]=MS::SCAN_NUMBER;
5887 0 : columns[2]=MS::FIELD_ID;
5888 0 : columns[3]=MS::DATA_DESC_ID;
5889 0 : columns[4]=MS::TIME;
5890 0 : vs_p->resetVisIter(columns,0.0);
5891 :
5892 0 : vs_p->initCalSet(calSet);
5893 0 : return true;
5894 0 : }
5895 : else {
5896 0 : throw(AipsError("Calibrater cannot initCalSet"));
5897 : return false;
5898 : }
5899 : }
5900 :
5901 0 : Bool OldCalibrater::cleanupVisSet() {
5902 :
5903 : // logSink() << LogOrigin("OldCalibrater","cleanupVisSet") << LogIO::NORMAL;
5904 :
5905 0 : if(vs_p) delete vs_p; vs_p=0;
5906 :
5907 : // Delete chanmask
5908 0 : initChanMask();
5909 :
5910 0 : return true;
5911 :
5912 : }
5913 :
5914 0 : VisibilityIterator::DataColumn OldCalibrater::configureForCorrection ()
5915 : {
5916 0 : if (!ok())
5917 0 : throw(AipsError("Calibrater not prepared for correct!"));
5918 :
5919 : // Nominally, we write out to the CORRECTED_DATA, unless absent
5920 0 : VisibilityIterator::DataColumn whichOutCol(VisibilityIterator::Corrected);
5921 :
5922 0 : if (!ms_p->tableDesc().isColumn("CORRECTED_DATA"))
5923 0 : throw(AipsError("CORRECTED_DATA column unexpectedly absent. Cannot correct."));
5924 :
5925 : // Ensure apply list non-zero and properly sorted
5926 0 : ve_p->setapply(vc_p);
5927 :
5928 : // Report the types that will be applied
5929 0 : applystate();
5930 :
5931 : // Arrange for iteration over data
5932 0 : Block<Int> columns;
5933 : // include scan iteration
5934 0 : columns.resize(5);
5935 0 : columns[0]=MS::ARRAY_ID;
5936 0 : columns[1]=MS::SCAN_NUMBER;
5937 0 : columns[2]=MS::FIELD_ID;
5938 0 : columns[3]=MS::DATA_DESC_ID;
5939 0 : columns[4]=MS::TIME;
5940 :
5941 : // Reset the VisibilityIterator in the VisSet.
5942 0 : vs_p->resetVisIter (columns, 0.0);
5943 :
5944 0 : return whichOutCol;
5945 0 : }
5946 :
5947 0 : void OldCalibrater::selectChannel(const String& spw) {
5948 :
5949 : // Initialize the chanmask_
5950 0 : initChanMask();
5951 :
5952 0 : if (mss_p && mssel_p) {
5953 :
5954 : // Refresh the frequencySelections object to feed to VI2, if relevant
5955 0 : frequencySelections_p.reset(new vi::FrequencySelections());
5956 :
5957 0 : vi::FrequencySelectionUsingChannels usingChannels;
5958 0 : usingChannels.add(*mss_p,mssel_p);
5959 0 : frequencySelections_p->add(usingChannels);
5960 :
5961 : // cout << usingChannels.toString() << endl;
5962 : // cout << "FS.size() = " << frequencySelections_p->size() << endl;
5963 :
5964 0 : }
5965 :
5966 0 : Matrix<Int> chansel = getChanIdx(spw);
5967 0 : uInt nselspw=chansel.nrow();
5968 :
5969 0 : if (nselspw==0)
5970 : logSink() << "Frequency selection: Selecting all channels in all spws."
5971 0 : << LogIO::POST;
5972 : else {
5973 :
5974 0 : logSink() << "Frequency selection: " << LogIO::POST;
5975 :
5976 : // Trap non-unit step (for now)
5977 0 : if (ntrue(chansel.column(3)==1)!=nselspw) {
5978 : logSink() << LogIO::WARN
5979 : << "Calibration does not support non-unit channel stepping; "
5980 : << "using step=1."
5981 0 : << LogIO::POST;
5982 0 : chansel.column(3)=1;
5983 : }
5984 :
5985 0 : Int nspw=vs_p->numberSpw();
5986 0 : Vector<Int> nChan0;
5987 0 : nChan0 = vs_p->numberChan();
5988 :
5989 0 : Vector<Int> uspw(chansel.column(0));
5990 0 : Vector<Int> ustart(chansel.column(1));
5991 0 : Vector<Int> uend(chansel.column(2));
5992 :
5993 0 : Vector<Int> start(nspw,INT_MAX);
5994 0 : Vector<Int> end(nspw,-INT_MAX);
5995 0 : logSink() << LogIO::NORMAL;
5996 0 : for (uInt i=0;i<nselspw;++i) {
5997 :
5998 0 : Int& spw=uspw(i);
5999 :
6000 : // Initialize this spw mask, if necessary (def = masked)
6001 0 : if (!chanmask_[spw])
6002 0 : chanmask_[spw]=new Vector<Bool>(nChan0(spw),true);
6003 :
6004 : // revise net start/end/nchan
6005 0 : start(spw)=min(start(spw),ustart(i));
6006 0 : end(spw)=max(end(spw),uend(i));
6007 0 : Int nchan=end(spw)-start(spw)+1; // net inclusive nchan
6008 :
6009 : // User's
6010 0 : Int step=chansel(i,3);
6011 0 : Int unchan=uend(i)-ustart(i)+1;
6012 :
6013 : // Update the mask (false = valid)
6014 0 : (*chanmask_[spw])(Slice(ustart(i),unchan))=false;
6015 :
6016 :
6017 : logSink() << ". Spw " << spw << ":"
6018 0 : << ustart(i) << "~" << uend(i)
6019 0 : << " (" << uend(i)-ustart(i)+1 << " channels,"
6020 : << " step by " << step << ")"
6021 0 : << endl;
6022 :
6023 : /*
6024 : cout << i << " " << spw << " {"
6025 : << start(spw) << " [" << ustart(i) << " "
6026 : << uend(i) << "] " << end(spw) << "}" << endl;
6027 : cout << "chanmask = ";
6028 : for (Int j=0;j<nChan0(spw);++j) cout << (*chanmask_[spw])(j);
6029 : cout << endl << endl;
6030 : */
6031 :
6032 : // Call via VisSet (avoid call to VisIter::origin)
6033 0 : vs_p->selectChannel(1,start(spw),nchan,step,spw,false);
6034 :
6035 : } // i
6036 0 : logSink() << LogIO::POST;
6037 :
6038 0 : } // non-triv spw selection
6039 :
6040 : // For testing:
6041 : if (false) {
6042 :
6043 : VisIter& vi(vs_p->iter());
6044 : VisBuffer vb(vi);
6045 :
6046 : // Pass each timestamp (VisBuffer) to VisEquation for correction
6047 : for (vi.originChunks(); vi.moreChunks(); vi.nextChunk()) {
6048 : vi.origin();
6049 : // for (vi.origin(); vi.more(); vi++)
6050 : cout << vb.spectralWindow() << " "
6051 : << vb.nChannel() << " "
6052 : << vb.channel() << " "
6053 : << vb.visCube().shape()
6054 : << endl;
6055 : }
6056 : }
6057 :
6058 0 : }
6059 :
6060 0 : void OldCalibrater::initChanMask() {
6061 :
6062 0 : for (uInt i=0;i<chanmask_.nelements();++i)
6063 0 : if (chanmask_[i])
6064 0 : delete chanmask_[i];
6065 0 : if (vs_p) {
6066 0 : chanmask_.resize(vs_p->numberSpw(),true);
6067 0 : chanmask_=NULL;
6068 : }
6069 : else {
6070 : // throw(AipsError("Trouble sizing chanmask!"));
6071 : // just don't support channel masking:
6072 0 : chanmask_.resize(0,true);
6073 : }
6074 :
6075 0 : }
6076 :
6077 : // Select on channel in the VisSet
6078 0 : void OldCalibrater::selectChannel(const String& mode,
6079 : const Int& nchan,
6080 : const Int& start, const Int& step,
6081 : const MRadialVelocity& mStart,
6082 : const MRadialVelocity& mStep) {
6083 :
6084 : // Set data selection variables
6085 0 : dataMode_p=mode;
6086 0 : dataNchan_p=nchan;
6087 0 : if (dataNchan_p<0) dataNchan_p=0;
6088 0 : dataStart_p=start;
6089 0 : if (dataStart_p<0) dataNchan_p=0;
6090 0 : dataStep_p=step;
6091 0 : if (dataStep_p<1) dataNchan_p=1;
6092 :
6093 0 : mDataStart_p=mStart;
6094 0 : mDataStep_p=mStep;
6095 :
6096 : // Select on frequency channel
6097 0 : if(dataMode_p=="channel") {
6098 : // *** this bit here is temporary till we unifomize data selection
6099 : //Getting the selected SPWs
6100 0 : MSMainColumns msc(*mssel_p);
6101 0 : Vector<Int> dataDescID = msc.dataDescId().getColumn();
6102 : Bool dum;
6103 0 : Sort sort( dataDescID.getStorage(dum),sizeof(Int) );
6104 0 : sort.sortKey((uInt)0,TpInt);
6105 0 : Vector<uInt> index,uniq;
6106 0 : sort.sort(index,dataDescID.nelements());
6107 0 : uInt nSpw = sort.unique(uniq,index);
6108 :
6109 0 : Vector<Int> selectedSpw(nSpw);
6110 0 : Vector<Int> nChan(nSpw);
6111 0 : for (uInt k=0; k < nSpw; ++k) {
6112 0 : selectedSpw[k]=dataDescID[index[uniq[k]]];
6113 0 : nChan[k]=vs_p->numberChan()(selectedSpw[k]);
6114 :
6115 : }
6116 0 : if(dataNchan_p==0) dataNchan_p=vs_p->numberChan()(selectedSpw[0]);
6117 0 : if(dataStart_p<0) {
6118 : logSink() << LogIO::SEVERE << "Illegal start pixel = "
6119 0 : << dataStart_p << LogIO::POST;
6120 : }
6121 0 : Int end = Int(dataStart_p) + Int(dataNchan_p) * Int(dataStep_p);
6122 0 : for (uInt k=0; k < selectedSpw.nelements() ; ++k){
6123 0 : if(end < 1 || end > nChan[k]) {
6124 : logSink() << LogIO::SEVERE << "Illegal step pixel = " << dataStep_p
6125 0 : << " in Spw " << selectedSpw[k]
6126 0 : << LogIO::POST;
6127 : }
6128 : logSink() << "Selecting "<< dataNchan_p
6129 : << " channels, starting at visibility channel "
6130 : << dataStart_p << " stepped by "
6131 0 : << dataStep_p << " in Spw " << selectedSpw[k] << LogIO::POST;
6132 :
6133 : // Set frequency channel selection for all spectral window id's
6134 : Int nch;
6135 : //Vector<Int> nChan=vs_p->numberChan();
6136 : //Int nSpw=vs_p->numberSpw();
6137 0 : if (dataNchan_p==0) {
6138 0 : nch=nChan(k);
6139 : }else {
6140 0 : nch=dataNchan_p;
6141 : };
6142 0 : vs_p->selectChannel(1,dataStart_p,nch,dataStep_p,selectedSpw[k]);
6143 :
6144 : }
6145 0 : }
6146 : // Select on velocity
6147 0 : else if (dataMode_p=="velocity") {
6148 0 : MVRadialVelocity mvStart(mDataStart_p.get("m/s"));
6149 0 : MVRadialVelocity mvStep(mDataStep_p.get("m/s"));
6150 : MRadialVelocity::Types
6151 0 : vType((MRadialVelocity::Types)mDataStart_p.getRefPtr()->getType());
6152 : logSink() << "Selecting "<< dataNchan_p
6153 : << " channels, starting at radio velocity " << mvStart
6154 : << " stepped by " << mvStep << ", reference frame is "
6155 0 : << MRadialVelocity::showType(vType) << LogIO::POST;
6156 0 : vs_p->iter().selectVelocity(Int(dataNchan_p), mvStart, mvStep,
6157 : vType, MDoppler::RADIO);
6158 0 : }
6159 :
6160 : // Select on optical velocity
6161 0 : else if (dataMode_p=="opticalvelocity") {
6162 0 : MVRadialVelocity mvStart(mDataStart_p.get("m/s"));
6163 0 : MVRadialVelocity mvStep(mDataStep_p.get("m/s"));
6164 : MRadialVelocity::Types
6165 0 : vType((MRadialVelocity::Types)mDataStart_p.getRefPtr()->getType());
6166 : logSink() << "Selecting "<< dataNchan_p
6167 : << " channels, starting at optical velocity " << mvStart
6168 : << " stepped by " << mvStep << ", reference frame is "
6169 0 : << MRadialVelocity::showType(vType) << LogIO::POST;
6170 0 : vs_p->iter().selectVelocity(Int(dataNchan_p), mvStart, mvStep,
6171 : vType, MDoppler::OPTICAL);
6172 0 : }
6173 :
6174 :
6175 0 : }
6176 :
6177 0 : Bool OldCalibrater::ok() {
6178 :
6179 0 : if(vs_p && ms_p && mssel_p && ve_p) {
6180 0 : return true;
6181 : }
6182 : else {
6183 0 : logSink() << "Calibrater is not yet initialized" << LogIO::POST;
6184 0 : return false;
6185 : }
6186 : }
6187 :
6188 0 : Bool OldCalibrater::genericGatherAndSolve() {
6189 :
6190 : #ifdef _OPENMP
6191 0 : Double Tsetup(0.0),Tgather(0.0),Tsolve(0.0);
6192 0 : Double time0=omp_get_wtime();
6193 : #endif
6194 :
6195 : //cout << "Generic gather and solve." << endl;
6196 :
6197 : // Create the solver
6198 0 : VisCalSolver vcs;
6199 :
6200 : // Inform logger/history
6201 0 : logSink() << "Solving for " << svc_p->typeName()
6202 0 : << LogIO::POST;
6203 :
6204 : // Initialize the svc according to current VisSet
6205 : // (this counts intervals, sizes CalSet)
6206 0 : Vector<Int> nChunkPerSol;
6207 0 : Int nSol = svc_p->sizeUpSolve(*vs_p,nChunkPerSol);
6208 :
6209 : // Create the in-memory (New)CalTable
6210 0 : svc_p->createMemCalTable();
6211 :
6212 : // The iterator, VisBuffer
6213 0 : VisIter& vi(vs_p->iter());
6214 0 : VisBuffer vb(vi);
6215 :
6216 0 : Vector<Int> slotidx(vs_p->numberSpw(),-1);
6217 :
6218 : // We will remember which spws couldn't be processed
6219 0 : Vector<Bool> unsolspw(vi.numberSpw());
6220 0 : unsolspw.set(false);
6221 :
6222 : // Manage verbosity of partial channel averaging
6223 0 : Vector<Bool> verb(vi.numberSpw(),true);
6224 :
6225 0 : Vector<Int64> nexp(vi.numberSpw(),0), natt(vi.numberSpw(),0),nsuc(vi.numberSpw(),0);
6226 :
6227 : #ifdef _OPENMP
6228 0 : Tsetup+=(omp_get_wtime()-time0);
6229 : #endif
6230 :
6231 0 : Int nGood(0);
6232 0 : vi.originChunks();
6233 0 : for (Int isol=0;isol<nSol && vi.moreChunks();++isol) {
6234 :
6235 : #ifdef _OPENMP
6236 0 : time0=omp_get_wtime();
6237 : #endif
6238 :
6239 0 : nexp(vi.spectralWindow())+=1;
6240 :
6241 : // capture obs, scan info so we can set it later
6242 : // (and not rely on what the VB averaging code can't properly do)
6243 0 : Vector<Int> scv,obsv;
6244 0 : Int solscan=vi.scan(scv)(0),solobs=vi.observationId(obsv)(0);
6245 :
6246 : // Arrange to accumulate
6247 : // VisBuffAccumulator vba(vs_p->numberAnt(),svc_p->preavg(),false);
6248 0 : VisBuffGroupAcc vbga(vs_p->numberAnt(),vs_p->numberSpw(),vs_p->numberFld(),svc_p->preavg());
6249 :
6250 0 : for (Int ichunk=0;ichunk<nChunkPerSol(isol);++ichunk) {
6251 :
6252 : // Current _chunk_'s spw
6253 0 : Int spw(vi.spectralWindow());
6254 :
6255 : // Only accumulate for solve if we can pre-calibrate
6256 0 : if (ve_p->spwOK(spw)) {
6257 :
6258 : // Collapse each timestamp in this chunk according to VisEq
6259 : // with calibration and averaging
6260 0 : for (vi.origin(); vi.more(); vi++) {
6261 :
6262 : // Force read of the field Id
6263 0 : vb.fieldId();
6264 :
6265 : // Apply the channel mask (~no-op, if unnecessary)
6266 0 : svc_p->applyChanMask(vb);
6267 :
6268 : // This forces the data/model/wt I/O, and applies
6269 : // any prior calibrations
6270 0 : ve_p->collapse(vb);
6271 :
6272 : // If permitted/required by solvable component, normalize
6273 0 : if (svc_p->normalizable())
6274 0 : vb.normalize();
6275 :
6276 : // If this solve not freqdep, and channels not averaged yet, do so
6277 0 : if (!svc_p->freqDepMat() && vb.nChannel()>1)
6278 0 : vb.freqAveCubes();
6279 :
6280 0 : if (svc_p->freqDepPar() &&
6281 0 : svc_p->fsolint()!="none" &&
6282 0 : svc_p->fintervalCh()>0.0) {
6283 : // cout << "svc_p->currSpw() = " << svc_p->currSpw() << endl;
6284 0 : if (verb(spw))
6285 : logSink() << " Reducing nchan in spw "
6286 : << spw
6287 0 : << " from " << vb.nChannel();
6288 0 : vb.channelAve(svc_p->chanAveBounds(spw));
6289 :
6290 : // Kludge for 3.4 to reset corr-indep flag to correct channel axis shape
6291 : // (because we use vb.flag() below, rather than vb.flagCube())
6292 0 : vb.flag().assign(operator>(partialNTrue(vb.flagCube(),IPosition(1,0)),0UL));
6293 :
6294 0 : if (verb(spw)) {
6295 : logSink() << " to "
6296 0 : << vb.nChannel() << LogIO::POST;
6297 0 : verb(spw)=false; // suppress future verbosity in this spw
6298 : }
6299 : }
6300 :
6301 : // Accumulate collapsed vb in a time average
6302 : // (only if the vb contains any unflagged data)
6303 0 : if (nfalse(vb.flag())>0)
6304 0 : vbga.accumulate(vb);
6305 :
6306 : }
6307 : }
6308 : else
6309 : // This spw not accumulated for solve
6310 0 : unsolspw(spw)=true;
6311 :
6312 : // Advance the VisIter, if possible
6313 0 : if (vi.moreChunks()) vi.nextChunk();
6314 :
6315 : }
6316 :
6317 : // Finalize the averged VisBuffer
6318 0 : vbga.finalizeAverage();
6319 :
6320 : // Establish meta-data for this interval
6321 : // (some of this may be used _during_ solve)
6322 : // (this sets currSpw() in the SVC)
6323 0 : Bool vbOk=(vbga.nBuf()>0 && svc_p->syncSolveMeta(vbga));
6324 :
6325 0 : svc_p->overrideObsScan(solobs,solscan);
6326 :
6327 : #ifdef _OPENMP
6328 0 : Tgather+=(omp_get_wtime()-time0);
6329 0 : time0=omp_get_wtime();
6330 : #endif
6331 :
6332 0 : if (vbOk) {
6333 :
6334 : // Use spw of first VB in vbga
6335 : // TBD: (currSpw==thisSpw) here?? (I.e., use svc_p->currSpw()? currSpw is prot!)
6336 0 : Int thisSpw=svc_p->spwMap()(vbga(0).spectralWindow());
6337 :
6338 0 : natt(thisSpw)+=1;
6339 :
6340 0 : slotidx(thisSpw)++;
6341 :
6342 : // Make data amp- or phase-only, if needed
6343 0 : vbga.enforceAPonData(svc_p->apmode());
6344 :
6345 : // Select on correlation via weights, according to the svc
6346 0 : vbga.enforceSolveCorrWeights(svc_p->phandonly());
6347 :
6348 0 : if (svc_p->useGenericSolveOne()) {
6349 : // generic individual solve
6350 :
6351 : //cout << "Generic individual solve: isol=" << isol << endl;
6352 :
6353 : // First guess
6354 0 : svc_p->guessPar(vbga(0));
6355 :
6356 : // Solve for each parameter channel (in curr Spw)
6357 :
6358 : // (NB: force const version of nChanPar() [why?])
6359 : // for (Int ich=0;ich<((const SolvableVisCal*)svc_p)->nChanPar();++ich) {
6360 0 : Bool totalGoodSol(false);
6361 0 : for (Int ich=((const SolvableVisCal*)svc_p)->nChanPar()-1;ich>-1;--ich) {
6362 : // for (Int ich=0;ich<((const SolvableVisCal*)svc_p)->nChanPar();++ich) {
6363 :
6364 : // If pars chan-dep, SVC mechanisms for only one channel at a time
6365 0 : svc_p->markTimer();
6366 0 : svc_p->focusChan()=ich;
6367 :
6368 : // svc_p->state();
6369 :
6370 : // cout << "Starting solution..." << endl;
6371 :
6372 : // Pass VE, SVC, VB to solver
6373 0 : Bool goodSoln=vcs.solve(*ve_p,*svc_p,vbga);
6374 :
6375 : // cout << "goodSoln= " << boolalpha << goodSoln << endl;
6376 :
6377 : // If good...
6378 0 : if (goodSoln) {
6379 0 : totalGoodSol=true;
6380 :
6381 0 : svc_p->formSolveSNR();
6382 0 : svc_p->applySNRThreshold();
6383 :
6384 : // ..and file this solution in the correct slot
6385 0 : if (svc_p->freqDepPar())
6386 0 : svc_p->keep1(ich);
6387 : // svc_p->keep(slotidx(thisSpw));
6388 : // Int n=svc_p->nSlots(thisSpw);
6389 : // svc_p->printActivity(n,slotidx(thisSpw),vi.fieldId(),thisSpw,nGood);
6390 :
6391 : }
6392 : else {
6393 : // report where this failure occured
6394 0 : svc_p->currMetaNote();
6395 0 : if (svc_p->freqDepPar())
6396 : // We must record a flagged solution for this channel
6397 0 : svc_p->keep1(ich);
6398 : }
6399 :
6400 : } // parameter channels
6401 :
6402 0 : if (totalGoodSol) {
6403 0 : svc_p->keepNCT();
6404 0 : nsuc(thisSpw)+=1;
6405 : }
6406 :
6407 : // Count good solutions.
6408 0 : if (totalGoodSol) nGood++;
6409 :
6410 : }
6411 : else {
6412 : //cout << "Self-directed individual solve: isol=" << isol << endl;
6413 : // self-directed individual solve
6414 : // TBD: selfSolveOne should return T/F for "good"
6415 0 : svc_p->selfSolveOne(vbga);
6416 :
6417 : // File this solution in the correct slot of the CalSet
6418 0 : svc_p->keepNCT();
6419 :
6420 0 : nGood++;
6421 0 : nsuc(thisSpw)+=1;
6422 : }
6423 :
6424 : } // vbOK
6425 :
6426 : #ifdef _OPENMP
6427 0 : Tsolve+=(omp_get_wtime()-time0);
6428 : #endif
6429 :
6430 0 : } // isol
6431 :
6432 : logSink() << " Found good "
6433 0 : << svc_p->typeName() << " solutions in "
6434 : << nGood << " slots."
6435 0 : << LogIO::POST;
6436 : #ifdef _OPENMP
6437 : #ifdef REPORT_CAL_TIMING
6438 : cout << "OldCalibrater::genericGatherAndSolve Timing: " << endl;
6439 : cout << " setup=" << Tsetup
6440 : << " gather=" << Tgather
6441 : << " solve=" << Tsolve
6442 : << " total=" << Tsetup+Tgather+Tsolve
6443 : << endl;
6444 : #endif
6445 : #endif
6446 :
6447 0 : summarize_uncalspws(unsolspw, "solv");
6448 :
6449 : // Store whole of result in a caltable
6450 0 : if (nGood==0) {
6451 : logSink() << "No output calibration table written."
6452 0 : << LogIO::POST;
6453 : }
6454 : else {
6455 :
6456 : // TBD: Remove BPOLY specificity here
6457 0 : if (svc_p->typeName()!="BPOLY") {
6458 : // Do global post-solve tinkering (e.g., phase-only, normalization, etc.)
6459 0 : svc_p->globalPostSolveTinker();
6460 :
6461 : // write the table
6462 0 : svc_p->storeNCT();
6463 :
6464 : }
6465 : }
6466 :
6467 : // Fill activity record
6468 0 : actRec_=Record();
6469 0 : actRec_.define("origin","Calibrater::genericGatherAndSolve");
6470 0 : actRec_.define("nExpected",nexp);
6471 0 : actRec_.define("nAttempt",natt);
6472 0 : actRec_.define("nSucceed",nsuc);
6473 :
6474 : {
6475 0 : Record solveRec=svc_p->solveActionRec();
6476 0 : if (solveRec.nfields()>0)
6477 0 : actRec_.merge(solveRec);
6478 0 : }
6479 :
6480 :
6481 :
6482 0 : return true;
6483 :
6484 0 : }
6485 :
6486 :
6487 : } //# NAMESPACE CASA - END
|
