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