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