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