Line data Source code
1 : //# FringeJones.cc: Implementation of FringeJones
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2011
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This library is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU Library General Public License as published by
7 : //# the Free Software Foundation; either version 2 of the License, or (at your
8 : //# option) any later version.
9 : //#
10 : //# This library 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 Library General Public
13 : //# License for more details.
14 : //#
15 : //# You should have received a copy of the GNU Library General Public License
16 : //# along with this library; if not, write to the Free Software Foundation,
17 : //# Inc., 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 :
27 : #include <synthesis/MeasurementComponents/FringeJones.h>
28 : #include <msvis/MSVis/VisBuffer.h>
29 : #include <msvis/MSVis/VisBuffAccumulator.h>
30 : #include <casacore/ms/MeasurementSets/MSColumns.h>
31 : #include <synthesis/CalTables/CTIter.h>
32 : #include <synthesis/MeasurementEquations/VisEquation.h> // *
33 : #include <synthesis/MeasurementComponents/SolveDataBuffer.h>
34 : #include <synthesis/MeasurementComponents/MSMetaInfoForCal.h>
35 : #include <casacore/lattices/Lattices/ArrayLattice.h>
36 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
37 : #include <casacore/scimath/Mathematics/FFTServer.h>
38 :
39 : #include <casacore/casa/Arrays/ArrayMath.h>
40 : #include <casacore/casa/Arrays/MatrixMath.h>
41 : #include <casacore/casa/Arrays/ArrayLogical.h>
42 : #include <casacore/casa/BasicSL/String.h>
43 : #include <casacore/casa/Utilities/Assert.h>
44 : #include <casacore/casa/Exceptions/Error.h>
45 : #include <casacore/casa/System/Aipsrc.h>
46 :
47 : #include <sstream>
48 :
49 : #include <casacore/measures/Measures/MCBaseline.h>
50 : #include <casacore/measures/Measures/MDirection.h>
51 : #include <casacore/measures/Measures/MEpoch.h>
52 : #include <casacore/measures/Measures/MeasTable.h>
53 :
54 : #include <casacore/casa/Logging/LogMessage.h>
55 : #include <casacore/casa/Logging/LogSink.h>
56 :
57 : #include <casacore/casa/Arrays/MaskedArray.h>
58 : #include <casacore/casa/Arrays/MaskArrMath.h>
59 :
60 : #include <gsl/gsl_rng.h>
61 : #include <gsl/gsl_randist.h>
62 : #include <gsl/gsl_vector.h>
63 : #include <gsl/gsl_blas.h>
64 : #include <gsl/gsl_spblas.h>
65 : #include <gsl/gsl_multilarge_nlinear.h>
66 : #include <gsl/gsl_multimin.h>
67 : #include <gsl/gsl_linalg.h>
68 : #include <iomanip> // needed for setprecision
69 :
70 : // DEVDEBUG gates the development debugging information to standard
71 : // error; it should be set to 0 for production.
72 : #define DEVDEBUG false
73 : #define KDISPSCALE 1e6
74 :
75 : using namespace casa::vi;
76 : using namespace casacore;
77 :
78 : namespace casa { //# NAMESPACE CASA - BEGIN
79 :
80 : // Start of GSL compliant solver
81 :
82 : class AuxParamBundle {
83 : public:
84 : SDBList &sdbs;
85 : size_t nCalls;
86 : private:
87 : // We make sure there are no copy or default constructors to
88 : // preserve the integrity of our reference member.
89 : AuxParamBundle();
90 : AuxParamBundle(AuxParamBundle const&);
91 : AuxParamBundle const& operator=(AuxParamBundle const&);
92 :
93 : size_t refant;
94 : size_t nCorrelations;
95 : size_t corrStep;
96 : Double t0;
97 : Double reftime;
98 : std::map< Int, std::set< Int > > activeAntennas;
99 : std::map< Int, Int > antennaIndexMap;
100 : // Can't I just have a vector, which maps indices to values anyway?
101 : std::vector<bool> parameterFlags;
102 : Int nParams;
103 : std::map< Int, Int > parameterMap;
104 : Int activeCorr;
105 : public:
106 219 : AuxParamBundle(SDBList& sdbs_, size_t refant, const std::map< Int, std::set<Int> >& activeAntennas_, Vector<Bool> paramActive) :
107 219 : sdbs(sdbs_),
108 219 : nCalls(0),
109 219 : refant(refant),
110 219 : nCorrelations(sdbs.nCorrelations() > 1 ? 2 : 1),
111 219 : corrStep(sdbs.nCorrelations() > 2 ? 3 : 1),
112 219 : activeAntennas(activeAntennas_),
113 219 : parameterFlags(),
114 219 : parameterMap(),
115 219 : activeCorr(-1)
116 : // corrStep(3)
117 : {
118 219 : Int last_index = sdbs.nSDB() - 1 ;
119 219 : t0 = sdbs(0).time()(0);
120 219 : Double tlast = sdbs(last_index).time()(0);
121 219 : reftime = 0.5*(t0 + tlast);
122 :
123 219 : parameterFlags = paramActive.tovector();
124 219 : Int j = 0; // The CASA parameter index (0=peculiar phase, 1=delay, 2=rate, 3=dispersive)
125 219 : Int i = 0; // the Least Squares parameter vector index, depending on what's being solved for
126 1095 : for (auto p=parameterFlags.begin(); p!=parameterFlags.end(); p++) {
127 876 : if (*p) {
128 653 : parameterMap.insert(std::pair<Int, Int>(j, i));
129 653 : i++;
130 : }
131 876 : j++;
132 : }
133 219 : if (i==0) {
134 0 : throw(AipsError("No parameters specified!"));
135 : }
136 219 : nParams = i; // There's always at least one parameter!
137 : // cerr << "AuxParamBundle reftime " << reftime << " t0 " << t0 <<" dt " << tlast - t0 << endl;
138 219 : }
139 :
140 3078672694 : Int nParameters() {
141 3078672694 : return nParams;
142 : }
143 6104 : Double get_t0() {
144 6104 : return t0;
145 : }
146 : Double
147 : get_ref_time() {
148 : return reftime;
149 : }
150 : size_t
151 575 : get_num_corrs() {
152 : //return sdbs.nCorrelations() > 1 ? 2 : 1;
153 575 : return nCorrelations;
154 : }
155 : size_t
156 1068 : get_num_antennas() {
157 1068 : if (activeCorr < 0) {
158 0 : throw(AipsError("Correlation out of range."));
159 : }
160 1068 : std::set< Int > ants = activeAntennas.find(activeCorr)->second;
161 2136 : return (size_t) ants.size();
162 1068 : }
163 : size_t
164 5512 : get_max_antenna_index() {
165 5512 : if (activeCorr < 0) {
166 0 : throw(AipsError("Correlation out of range."));
167 : }
168 5512 : return *(activeAntennas.find(activeCorr)->second.rbegin());
169 : }
170 : // Sometimes there is Int, sometimes size_t; the following ones are casacore::Int.
171 : Int
172 : get_num_data_points() {
173 : Int nTotalRows = 0;
174 : for (Int i = 0; i != sdbs.nSDB(); i++) {
175 : nTotalRows += sdbs(i).nRows();
176 : }
177 : return nTotalRows * sdbs.nChannels();
178 : }
179 : Int
180 356 : get_actual_num_data_points() {
181 356 : Int nTotalRows = 0;
182 42572 : for (Int i = 0; i != sdbs.nSDB(); i++) {
183 42216 : SolveDataBuffer& s (sdbs(i));
184 772166 : for (Int irow=0; irow!=s.nRows(); irow++) {
185 729950 : if (s.flagRow()(irow)) continue;
186 727956 : nTotalRows++;
187 : }
188 : }
189 356 : return nTotalRows * sdbs.nChannels();
190 : }
191 : size_t
192 14264156 : get_data_corr_index(size_t icorr) {
193 14264156 : if (icorr > nCorrelations) {
194 0 : throw(AipsError("Correlation out of range."));
195 : }
196 14264156 : size_t dcorr = icorr * corrStep;
197 14264156 : return dcorr;
198 : }
199 : bool
200 29709442 : isActive(size_t iant) {
201 29709442 : std::set<Int> ants = activeAntennas.find(activeCorr)->second;
202 29709442 : if (iant == refant) return true;
203 25866054 : else return (ants.find(iant) != ants.end());
204 29709442 : }
205 : Int
206 4870651592 : get_param_corr_param_index(size_t iant0, size_t ipar) {
207 4870651592 : if (iant0 == refant) return -1;
208 4109409584 : int iant1 = antennaIndexMap[iant0];
209 4109409584 : if (iant1 > antennaIndexMap[refant]) {
210 4109409584 : iant1 -= 1;
211 : }
212 : int ipar1;
213 4109409584 : auto p = parameterMap.find(ipar);
214 4109409584 : if (p==parameterMap.end()) {
215 1030739046 : ipar1 = -1;
216 : }
217 : else {
218 3078670538 : ipar1 = (iant1 * nParameters()) + p->second;
219 : }
220 4109409584 : return ipar1;
221 : }
222 : size_t
223 14264156 : get_active_corr() {
224 14264156 : return activeCorr;
225 : }
226 : void
227 356 : set_active_corr(size_t icorr) {
228 356 : activeCorr = icorr;
229 356 : antennaIndexMap.clear();
230 356 : Int i = 0;
231 356 : std::set<Int>::iterator it;
232 356 : std::set<Int> ants = activeAntennas.find(activeCorr)->second;
233 2512 : for (it = ants.begin(); it != ants.end(); it++) {
234 2156 : antennaIndexMap[*it] = i++;
235 : }
236 356 : }
237 : };
238 :
239 :
240 :
241 : void
242 0 : print_baselines(std::set<std::pair< Int, Int > > baselines) {
243 0 : cerr << "Baselines encountered ";
244 0 : std::set<std::pair< Int, Int > >::iterator it;
245 0 : for (it=baselines.begin(); it != baselines.end(); ++it) {
246 0 : cerr << "(" << it->first << ", " << it->second << ") ";
247 : }
248 0 : cerr << endl;
249 0 : }
250 :
251 :
252 : int
253 2874 : expb_f(const gsl_vector *param, void *d, gsl_vector *f)
254 : {
255 2874 : AuxParamBundle *bundle = (AuxParamBundle *)d;
256 2874 : SDBList& sdbs = bundle->sdbs;
257 2874 : Double refTime = bundle->get_t0();
258 :
259 2874 : gsl_vector_set_zero(f);
260 : // Vector<Double> freqs = sdbs.freqs();
261 :
262 2874 : const Double reffreq0=sdbs(0).freqs()(0); // First freq in first SDB
263 :
264 2874 : size_t count = 0; // This is the master index.
265 :
266 2874 : Double sumwt = 0.0;
267 2874 : Double xi_squared = 0.0;
268 :
269 319815 : for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++) {
270 316941 : SolveDataBuffer& s (sdbs(ibuf));
271 316941 : if (!s.Ok()) continue;
272 :
273 316637 : const Vector<Double> freqs(s.freqs()); // This ibuf's freqs
274 316637 : Float fmin_ = min(freqs);
275 316637 : Float fmax = max(freqs);
276 :
277 316637 : const Cube<Complex>& v(s.visCubeCorrected());
278 316637 : const Cube<Bool>& fl(s.flagCube());
279 316637 : const Cube<Float>& weights(s.weightSpectrum());
280 :
281 7640770 : for (Int irow=0; irow!=s.nRows(); irow++) {
282 7324133 : if (s.flagRow()(irow)) continue;
283 :
284 7316830 : Int ant1(s.antenna1()(irow));
285 7316830 : Int ant2(s.antenna2()(irow));
286 7316830 : if (!bundle->isActive(ant1) || !bundle->isActive(ant2))
287 2864 : continue;
288 7313966 : if (ant1==ant2) continue;
289 :
290 : // VisBuffer.h seems to suggest that a vb.visCube may have shape
291 : // (nCorr(), nChannel(), nRow())
292 6809389 : size_t icorr0 = bundle->get_active_corr();
293 6809389 : size_t dcorr = bundle->get_data_corr_index(icorr0);
294 : // We also need to get the right parameters for this,
295 : // polarization (icorr is an encoding of the
296 : // polarization of the correlation products).
297 :
298 : Double phi0, tau, r, disp;
299 : {
300 : Int i;
301 : Double phi0_1, tau1, r1, disp1;
302 : Double phi0_2, tau2, r2, disp2;
303 :
304 6809389 : phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
305 6809389 : tau1 = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
306 6809389 : r1 = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
307 6809389 : disp1 = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
308 :
309 6809389 : phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
310 6809389 : tau2 = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
311 6809389 : r2 = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
312 6809389 : disp2 = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
313 :
314 6809389 : phi0 = phi0_2 - phi0_1;
315 6809389 : tau = tau2 - tau1;
316 6809389 : r = r2-r1;
317 6809389 : disp = disp2-disp1;
318 : }
319 :
320 75203829 : for (size_t ichan = 0; ichan != v.ncolumn(); ichan++) {
321 68394440 : if (fl(dcorr, ichan, irow)) continue;
322 :
323 67827326 : Float freq = freqs(ichan);
324 67827326 : Float k_disp = KDISPSCALE*C::_2pi*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
325 :
326 67827326 : Complex vis = v(dcorr, ichan, irow);
327 67827326 : Double w0 = weights(dcorr, ichan, irow);
328 : // FIXME: what should we use to scale the weights?
329 : // Double weightScale = norm(vis);
330 : // Double weightScale = abs(vis);
331 : // Double weightScale = 1;
332 : // Double weightScale = 1/sqrt(w0); // Actually AIPS 0, not AIPS 1!
333 : // Double weightScale = pow(w0, -0.75); // AIPS 2
334 : // Double weightScale = 1/w0; // AIPS 3
335 : // Double weightScale = norm(vis); // Casa 1, I guess
336 67827326 : Double w = sqrt(w0);
337 67827326 : sumwt += w*w;
338 67827326 : if (fabs(w) < FLT_EPSILON) continue;
339 : // We have to turn the delay back into seconds from nanoseconds.
340 : // Freq difference is in Hz, which comes out typically as 1e6 bands
341 : //Double wDf = C::_2pi*(freqs(ichan) - freqs(0))*1e-9;
342 67827326 : Double wDf = C::_2pi*(freqs(ichan) - reffreq0)*1e-9;
343 : //
344 67827326 : Double t1 = s.time()(0);
345 : // FIXME: Remind me why we *do* scale wDf with 1e-9
346 : // but do *not* do that with ref_freq?
347 : // I have a theory which is mine:
348 : // this is because tau is in nanoseconds.
349 : //Double ref_freq = freqs(0);
350 : //Double wDt = C::_2pi*(t1 - refTime) * ref_freq;
351 67827326 : Double wDt = C::_2pi*(t1 - refTime) * reffreq0;
352 :
353 67827326 : Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp);
354 67827326 : Double vtheta = arg(vis);
355 :
356 67827326 : Double c_r = w*(cos(mtheta) - cos(vtheta));
357 67827326 : Double c_i = w*(sin(mtheta) - sin(vtheta));
358 67827326 : gsl_vector_set(f, count, c_r);
359 67827326 : gsl_vector_set(f, count+1, c_i);
360 :
361 67827326 : count += 2;
362 67827326 : xi_squared += c_r*c_r + c_i*c_i;
363 : }
364 : }
365 316637 : }
366 : // cerr << "Residual xi-squared = " << xi_squared << endl;
367 2874 : return GSL_SUCCESS;
368 : }
369 :
370 :
371 : int
372 2874 : expb_df(CBLAS_TRANSPOSE_t TransJ, const gsl_vector* param, const gsl_vector *u, void *bundle_, gsl_vector *v, gsl_matrix *JTJ)
373 : {
374 :
375 : // param is the current vector for which we're finding the jacobian.
376 : // if TransJ is true, evaluate J^T u and store in v.
377 : // Also store J^T . J in lower half of JTJ.
378 2874 : std::set <std::pair < Int, Int> > baselines;
379 2874 : AuxParamBundle *bundle = (AuxParamBundle *)bundle_;
380 :
381 2874 : SDBList& sdbs = bundle->sdbs;
382 2874 : const Double reffreq0=sdbs(0).freqs()(0); // First freq in first SDB
383 :
384 2874 : size_t count = 0; // This is the master index.
385 :
386 2874 : gsl_vector_set_zero(v);
387 2874 : gsl_matrix_set_zero(JTJ);
388 :
389 2874 : Double refTime = bundle->get_t0();
390 :
391 319815 : for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++) {
392 316941 : SolveDataBuffer& s (sdbs(ibuf));
393 316941 : if (!s.Ok()) continue;
394 :
395 316637 : const Vector<Double>& freqs(s.freqs()); // This ibuf's freqs
396 316637 : Float fmin_ = min(freqs);
397 316637 : Float fmax = max(freqs);
398 316637 : const Cube<Complex>& vis(s.visCubeCorrected());
399 316637 : const Cube<Bool>& fl(s.flagCube());
400 316637 : const Cube<Float>& weights(s.weightSpectrum());
401 :
402 316637 : Double t1 = s.time()(0);
403 7640770 : for (Int irow=0; irow!=s.nRows(); irow++) {
404 7831574 : if (s.flagRow()(irow)) continue;
405 :
406 7316830 : Int ant1(s.antenna1()(irow));
407 7316830 : Int ant2(s.antenna2()(irow));
408 :
409 7316830 : if (ant1==ant2) continue;
410 6809389 : if (!bundle->isActive(ant1) || !bundle->isActive(ant2)) {
411 0 : continue;
412 : }
413 :
414 : // VisBuffer.h seems to suggest that a vb.visCube may have shape
415 : // (nCorr(), nChannel(), nRow())
416 :
417 6809389 : size_t icorr0 = bundle->get_active_corr();
418 6809389 : size_t dcorr = bundle->get_data_corr_index(icorr0);
419 : // We also need to get the right parameters for this
420 : // polarization (icorr is an encoding of the
421 : // polarization of the correlation products).
422 :
423 : Double phi0, tau, r, disp;
424 : {
425 : Int i;
426 : Double phi0_1, tau1, r1, disp1;
427 : Double phi0_2, tau2, r2, disp2;
428 :
429 6809389 : phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
430 6809389 : tau1 = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
431 6809389 : r1 = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
432 6809389 : disp1 = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
433 :
434 6809389 : phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
435 6809389 : tau2 = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
436 6809389 : r2 = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
437 6809389 : disp2 = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
438 :
439 6809389 : phi0 = phi0_2 - phi0_1;
440 6809389 : tau = tau2 - tau1;
441 6809389 : r = r2-r1;
442 6809389 : disp = disp2-disp1;
443 : }
444 :
445 :
446 : //Double ref_freq = freqs(0);
447 : //Double wDt = C::_2pi*(t1 - refTime) * ref_freq;
448 6809389 : Double wDt = C::_2pi*(t1 - refTime) * reffreq0;
449 6809389 : bool found_data = false;
450 :
451 75203829 : for (size_t ichan = 0; ichan != vis.ncolumn(); ichan++) {
452 68394440 : if (fl(dcorr, ichan, irow)) continue;
453 67827326 : Double w0 = weights(dcorr, ichan, irow);
454 67827326 : Double w = sqrt(w0);
455 67827326 : if (fabs(w) < FLT_EPSILON) continue;
456 67827326 : found_data = true;
457 :
458 67827326 : Float freq = freqs(ichan);
459 67827326 : Float k_disp = KDISPSCALE*C::_2pi*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
460 :
461 : // Add a 1e-9 factor because tau parameter is in nanoseconds.
462 : //Double wDf = C::_2pi*(freqs(ichan) - freqs(0))*1e-9;
463 67827326 : Double wDf = C::_2pi*(freqs(ichan) - reffreq0)*1e-9;
464 : //
465 67827326 : Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp);
466 67827326 : Double ws = sin(mtheta);
467 67827326 : Double wc = cos(mtheta);
468 :
469 67827326 : Double p0 = 1.0;
470 67827326 : Double p1 = wDf;
471 67827326 : Double p2 = wDt;
472 67827326 : Double p3 = k_disp;
473 :
474 67827326 : Vector<Double> dterm2(4);
475 67827326 : dterm2(0) = -p0;
476 67827326 : dterm2(1) = -p1;
477 67827326 : dterm2(2) = -p2;
478 67827326 : dterm2(3) = -p3;
479 :
480 67827326 : Vector<Double> dterm1(4);
481 67827326 : dterm1(0) = p0;
482 67827326 : dterm1(1) = p1;
483 67827326 : dterm1(2) = p2;
484 67827326 : dterm1(3) = p3;
485 :
486 : /*
487 : What we want to express is just:
488 : J[count + 0, iparam2 + 0] = w*-ws*-1.0;
489 : J[count + 1, iparam2 + 0] = w*+wc*-1.0;
490 : J[count + 0, iparam2 + 1] = w*-ws*-wDf;
491 : J[count + 1, iparam2 + 1] = w*+wc*-wDf;
492 : J[count + 0, iparam2 + 2] = w*-ws*-wDt;
493 : J[count + 1, iparam2 + 2] = w*+wc*-wDt;
494 :
495 : But in the GSL multilarge framework we have to
496 : be ready to calculate either J*u for a given u
497 : or J^T*u, depending on the flag TransJ, and we also have to fill in the
498 :
499 : v[iparam + ...] = J[count + ..., iparam + ...] * u[iparam + ...]
500 :
501 : or
502 :
503 : v[iparam + ...] = J^T[iparam + ..., count + ...] * u[count + ...]
504 :
505 : <https://www.gnu.org/software/gsl/doc/html/nls.html#c.gsl_multifit_nlinear_default_parameters>
506 :
507 : "Additionally, the normal equations matrix J^T J should be stored in the lower half of JTJ."
508 :
509 : So we should also use
510 : JTJ[iparam + ..., iparam + ...] += J^T[iparam + ..., count + ...] J[count + ..., iparam + ...]
511 :
512 : */
513 67827326 : if (TransJ==CblasNoTrans) {
514 0 : for (Int di=0; di<4; di++) {
515 : Int i;
516 0 : if ((i = bundle->get_param_corr_param_index(ant2, di))>=0) {
517 0 : (*gsl_vector_ptr(v, count + 0)) += (w*-ws*dterm2(di)) * gsl_vector_get(u, i);
518 0 : (*gsl_vector_ptr(v, count + 1)) += (w*+wc*dterm2(di)) * gsl_vector_get(u, i);
519 : }
520 0 : if ((i = bundle->get_param_corr_param_index(ant1, di))>=0) {
521 0 : (*gsl_vector_ptr(v, count + 0)) += gsl_vector_get(u, i) * (w*-ws*dterm1(di));
522 0 : (*gsl_vector_ptr(v, count + 1)) += gsl_vector_get(u, i) * (w*+wc*dterm1(di));
523 : }
524 : }
525 : } else {
526 339136630 : for (Int di=0; di<4; di++) {
527 : Int i;
528 271309304 : if ((i = bundle->get_param_corr_param_index(ant2, di))>=0) {
529 202674764 : (*gsl_vector_ptr(v, i)) += (w*-ws*dterm2(di)) * gsl_vector_get(u, count + 0);
530 202674764 : (*gsl_vector_ptr(v, i)) += (w*+wc*dterm2(di)) * gsl_vector_get(u, count + 1);
531 : }
532 271309304 : if ((i = bundle->get_param_corr_param_index(ant1, di))>=0) {
533 137427776 : (*gsl_vector_ptr(v, i)) += gsl_vector_get(u, count + 0) * (w*-ws*dterm1(di));
534 137427776 : (*gsl_vector_ptr(v, i)) += gsl_vector_get(u, count + 1) * (w*+wc*dterm1(di));
535 : }
536 : }
537 : }
538 67827326 : if (JTJ) {
539 : Int i, j;
540 67827326 : Double wterm = (-ws) * (-ws) + (+wc) * (+wc);
541 67827326 : if (fabs(1-wterm) > 1e-15)
542 0 : throw AipsError("Insufficiently at one");
543 339136630 : for (Int di=0; di<4; di++) {
544 949582564 : for (Int dj=0; dj<=di; dj++) {
545 1086665224 : if (((i = bundle->get_param_corr_param_index(ant2, di))>=0) &&
546 408391964 : ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
547 405865916 : (*gsl_matrix_ptr(JTJ, i, j)) += w0*dterm2(di)*dterm2(dj);
548 : }
549 954169748 : if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
550 275896488 : ((j = bundle->get_param_corr_param_index(ant1, dj))>=0)) {
551 275054472 : (*gsl_matrix_ptr(JTJ, i, j)) += w0*dterm1(di)*dterm1(dj);
552 : }
553 : }
554 : }
555 : // iant1 != iant2, so we don't have to worry about collisions
556 339136630 : for (Int di=0; di<4; di++) {
557 1356546520 : for (Int dj=0; dj<4; dj++) {
558 : Int i0, j0;
559 1634948320 : if (((i0 = bundle->get_param_corr_param_index(ant1, di))>=0) &&
560 549711104 : ((j0 = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
561 412681168 : Int i1 = max(i0, j0);
562 412681168 : Int j1 = min(i0, j0);
563 412681168 : (*gsl_matrix_ptr(JTJ, i1, j1)) += w0*dterm2(di)*dterm1(dj);
564 : }
565 : }
566 : }
567 67827326 : count += 2;
568 : } // if JTJ
569 67827326 : } // loop over channels
570 6809389 : if (found_data) {
571 6809389 : std::pair<Int, Int> antpair = std::make_pair(ant1, ant2);
572 6809389 : bool newBaseline = (baselines.find(antpair) == baselines.end());
573 6809389 : if (newBaseline) {
574 76022 : baselines.insert(antpair);
575 : // cerr << "paramFlagging for antenna "<< ant1 << ": ";
576 : // for (size_t di=0; di<4; di++) {
577 : // cerr << (bundle->get_param_corr_param_index(ant1, di)>=0) << " ";
578 : // }
579 : // cerr << endl;
580 : // cerr << "indices for antenna "<< ant1 << ": ";
581 : // if (bundle->get_param_corr_param_index(ant1, 0) >= 0) {
582 : // for (size_t di=0; di<4; di++) {
583 : // cerr << bundle->get_param_corr_param_index(ant1, di) << " ";
584 : // }
585 : // cerr << endl;
586 : // }
587 : // cerr << "phi0 " << phi0 << " tau " << tau << " r " << r << endl;
588 : }
589 : }
590 : }
591 : }
592 : if (DEVDEBUG && 0) {
593 : print_baselines(baselines);
594 : cerr << "count " << count << endl;
595 : cerr << "v = ";
596 : for (size_t i=0; i!=v->size; i++) {
597 : cerr << gsl_vector_get(v, i) << " ";
598 : }
599 : cerr << endl;
600 : // if (JTJ) {
601 : // cerr <<"JTJ " << std::scientific << endl;
602 : // for (size_t i=0; i!=JTJ->size1; i++) {
603 : // for (size_t j=0; j!=JTJ->size2; j++) {
604 : // cerr << gsl_matrix_get(JTJ, i, j) << " ";
605 : // }
606 : // cerr << endl;
607 : // }
608 : // cerr << endl;
609 : // }
610 : }
611 2874 : return GSL_SUCCESS;
612 2874 : }
613 :
614 :
615 :
616 : int
617 356 : expb_hess(gsl_vector *param, AuxParamBundle *bundle, gsl_matrix *hess, Double xi_squared, gsl_vector *snr_vector, LogIO& logSink)
618 : {
619 : // We calculate the diagonal for the hessian as used by AIPS for
620 : // the signal to noise. The AIPS formulation, by Fred Schwab, is a
621 : // hand-rolled routine that solves a different problem to ours: by
622 : // using a triangular matrix for the Jacobian (requiring antenna i<
623 : // antenna j) the J^T * J term vanishes throughout and the Hessian
624 : // of the Xi^2 functional *only* includes the second-order
625 : // derivative terms, which are usually neglected.
626 : //
627 : // This is very clever, but it also means different covariance and
628 : // information matrices, and therefore a different SNR. Here we
629 : // use a generic least squares solver but we cheat slightly and use
630 : // the AIPS form for the Hessian and SNR calculations.
631 : //
632 : // FIXME: Is there any compelling reason to use gsl_vectors for
633 : // this, given that we're not really hooked in to the gsl
634 : // least squares framework by the time we do this?
635 :
636 356 : SDBList& sdbs = bundle->sdbs;
637 356 : Double refTime = bundle->get_t0();
638 :
639 : // Dimensions of (num_antennas); is the same dimension as
640 : // param vector here.
641 356 : gsl_matrix_set_zero(hess);
642 :
643 356 : const Double reffreq0=sdbs(0).freqs()(0); // First freq in first SDB
644 :
645 356 : size_t nobs = 0;
646 356 : Double sumwt = 0;
647 356 : size_t numpar = param->size;
648 :
649 42572 : for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++)
650 : {
651 42216 : SolveDataBuffer& s (sdbs(ibuf));
652 42216 : if (!s.Ok()) continue;
653 :
654 42152 : const Vector<Double> freqs(s.freqs()); // This ibuf's freqs
655 42152 : Float fmin_ = min(freqs);
656 42152 : Float fmax = max(freqs);
657 :
658 42152 : const Cube<Complex>& v(s.visCubeCorrected());
659 42152 : const Cube<Bool>& fl(s.flagCube());
660 42152 : const Cube<Float>& weights(s.weightSpectrum());
661 :
662 :
663 771718 : for (Int irow=0; irow!=s.nRows(); irow++) {
664 729566 : if (s.flagRow()(irow)) continue;
665 :
666 727892 : Int ant1(s.antenna1()(irow));
667 727892 : Int ant2(s.antenna2()(irow));
668 727892 : if (!bundle->isActive(ant1) || !bundle->isActive(ant2))
669 716 : continue;
670 727176 : if (ant1==ant2) continue;
671 :
672 : // VisBuffer.h seems to suggest that a vb.visCube may have shape
673 : // (nCorr(), nChannel(), nRow())
674 645378 : size_t icorr0 = bundle->get_active_corr();
675 645378 : size_t dcorr = bundle->get_data_corr_index(icorr0);
676 : // We also need to get the right parameters for this,
677 : // polarization (icorr is an encoding of the
678 : // polarization of the correlation products).
679 :
680 : Double phi0, tau, r, disp;
681 : {
682 : Int i;
683 : Double phi0_1, tau1, r1, disp1;
684 : Double phi0_2, tau2, r2, disp2;
685 :
686 645378 : phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
687 645378 : tau1 = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
688 645378 : r1 = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
689 645378 : disp1 = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
690 :
691 :
692 645378 : phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
693 645378 : tau2 = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
694 645378 : r2 = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
695 645378 : disp2 = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
696 :
697 645378 : phi0 = phi0_2 - phi0_1;
698 645378 : tau = tau2 - tau1;
699 645378 : r = r2-r1;
700 645378 : disp = disp2-disp1;
701 : }
702 :
703 7913458 : for (size_t ichan = 0; ichan != v.ncolumn(); ichan++) {
704 7268080 : if (fl(dcorr, ichan, irow)) continue;
705 7217052 : Complex vis = v(dcorr, ichan, irow);
706 : // Fixme: this isn't a square root.
707 7217052 : Double w0 = weights(dcorr, ichan, irow);
708 7217052 : sumwt += w0;
709 7217052 : Double w = w0;
710 7217052 : nobs++;
711 7217052 : if (fabs(w) < FLT_EPSILON) continue;
712 :
713 : // We have to turn the delay back into seconds from nanoseconds.
714 : // Freq difference is in Hz, which comes out typically as 1e6 bands
715 : //Double wDf = C::_2pi*(freqs(ichan) - freqs(0))*1e-9;
716 7217052 : Double wDf = C::_2pi*(freqs(ichan) - reffreq0)*1e-9;
717 : //
718 7217052 : Double t1 = s.time()(0);
719 :
720 : //Double ref_freq = freqs(0);
721 : //Double wDt = C::_2pi*(t1 - refTime) * ref_freq;
722 7217052 : Double wDt = C::_2pi*(t1 - refTime) * reffreq0;
723 :
724 7217052 : Float freq = freqs(ichan);
725 7217052 : Float k_disp = KDISPSCALE*C::_2pi*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
726 :
727 7217052 : Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp);
728 7217052 : Double vtheta = arg(vis);
729 :
730 : // Hold on a minute though!
731 7217052 : Double cx = w*cos(vtheta - mtheta);
732 :
733 7217052 : Matrix<Double> dterm(4,4);
734 7217052 : dterm(0, 0) = cx;
735 7217052 : dterm(0, 1) = wDf*cx;
736 7217052 : dterm(0, 2) = wDt*cx;
737 7217052 : dterm(0, 3) = k_disp*dterm(0, 1);
738 7217052 : dterm(1, 1) = wDf*dterm(0, 1);
739 7217052 : dterm(1, 2) = wDt*dterm(0, 1);
740 7217052 : dterm(1, 3) = wDf*dterm(0, 3);
741 7217052 : dterm(2, 2) = wDt*dterm(1, 2);
742 7217052 : dterm(2, 3) = wDt*dterm(1, 3);
743 7217052 : dterm(3, 3) = k_disp*dterm(2, 3);
744 :
745 : // Symmetry terms:
746 7217052 : dterm(1, 0) = dterm(0, 1);
747 7217052 : dterm(2, 0) = dterm(0, 2);
748 7217052 : dterm(3, 0) = dterm(0, 3);
749 7217052 : dterm(2, 1) = dterm(1, 2);
750 7217052 : dterm(3, 1) = dterm(1, 3);
751 7217052 : dterm(3, 2) = dterm(2, 3);
752 :
753 :
754 :
755 36085260 : for (Int di=0; di<4; di++) {
756 144341040 : for (Int dj=0; dj<4; dj++) {
757 : Int i, j;
758 168331904 : if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
759 52859072 : ((j = bundle->get_param_corr_param_index(ant1, dj))>=0)) {
760 39719696 : *gsl_matrix_ptr(hess, i, j) += dterm(di, dj);
761 : }
762 : // Exactly the same logic, but with antenna2
763 201447456 : if (((i = bundle->get_param_corr_param_index(ant2, di))>=0) &&
764 85974624 : ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
765 64675800 : *gsl_matrix_ptr(hess, i, j) += dterm(di, dj);
766 : }
767 : }
768 : }
769 : // FIXME: Not just diagonal terms any more!
770 : // Note that some of these are not in the lower
771 : // triangular part, even though they are copied
772 : // faithfully from AIPS which thinks it is filling
773 : // a triangular matrix and handles symmetry
774 : // later. Unless I've missed something (again).
775 36085260 : for (Int di=0; di<4; di++) {
776 144341040 : for (Int dj=0; dj<4; dj++) {
777 : Int i, j;
778 168331904 : if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
779 52859072 : ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
780 39719696 : *gsl_matrix_ptr(hess, i, j) -= dterm(di, dj);
781 39719696 : *gsl_matrix_ptr(hess, j, i) -= dterm(dj, di);
782 : } // if
783 : } // dj
784 : } // di
785 7217052 : } // ichan
786 : } // irow
787 42152 : } // ibuff
788 : // s is more tricky: it is the xi^2 term from exp_f
789 :
790 356 : xi_squared = max(xi_squared, 1e-25);
791 :
792 : if (DEVDEBUG && 0) {
793 : cerr << "The matrix is" << endl;
794 : cerr << setprecision(3) << scientific;
795 : for (size_t i = 0; i != hess->size1; i++) {
796 : for (size_t j = 0; j < hess->size2; j++) {
797 : cerr << gsl_matrix_get(hess,i,j) << " ";
798 : }
799 : cerr << endl;
800 : }
801 : cerr << endl;
802 : // str.unsetf(cerr:floatfield);
803 : cerr << std::defaultfloat;
804 : }
805 : //
806 356 : size_t hsize = hess->size1;
807 : int signum;
808 356 : gsl_permutation *perm = gsl_permutation_alloc (hsize);
809 356 : gsl_matrix *lu = gsl_matrix_alloc(hsize, hsize);
810 356 : gsl_matrix *inv_hess = gsl_matrix_alloc(hsize, hsize);
811 :
812 356 : gsl_linalg_LU_decomp(hess, perm, &signum);
813 356 : Double det = gsl_linalg_LU_det(hess, signum);
814 356 : if ((fabs(det) < GSL_DBL_EPSILON) || std::isnan(det)) {
815 0 : logSink << "Hessian matrix singular (determinant=" << det << "); setting signal-to-noise ratio to zero." << LogIO::POST;
816 : // Singular matrix; fill snrs with zero.
817 0 : for (size_t i=0; i < hess->size1; i+=bundle->nParameters()) {
818 0 : Double snr = 0;
819 0 : gsl_vector_set(snr_vector, i, snr);
820 : }
821 : }
822 : else {
823 356 : gsl_linalg_LU_invert(hess, perm, inv_hess);
824 :
825 356 : Double sigma2 = xi_squared / (nobs - numpar) * nobs / sumwt;
826 : // cerr << "xi_squared " << xi_squared << " Nobs " << nobs << " sumwt " << sumwt << " sigma2 " << sigma2 << endl;
827 2156 : for (size_t i=0; i < hess->size1; i+=bundle->nParameters()) {
828 1800 : Double h = gsl_matrix_get(inv_hess, i, i);
829 1800 : Double snr0 = sqrt(sigma2*h*0.5);
830 1800 : snr0 = min(snr0, 9999.999);
831 1800 : Double snr = (snr0 > 1e-20) ? 1.0/snr0 : snr0;
832 1800 : gsl_vector_set(snr_vector, i, snr);
833 : }
834 : }
835 356 : gsl_matrix_free(lu);
836 356 : gsl_matrix_free(inv_hess);
837 356 : gsl_permutation_free(perm);
838 : // SNR[i], according to aips, is 1/sqrt(sigma2*hess(i1,i1)*0.5);
839 : // Note that in AIPS i1 ranges over 1..NANT
840 : // We use 1 as a success code.
841 356 : return 1;
842 : }
843 :
844 :
845 : Int
846 219 : findRefAntWithData(SDBList& sdbs, Vector<Int>& refAntList, Int prtlev) {
847 219 : std::set<Int> activeAntennas;
848 26754 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
849 26535 : SolveDataBuffer& s(sdbs(ibuf));
850 26535 : if (!s.Ok())
851 32 : continue;
852 26503 : Cube<Bool> fl = s.flagCube();
853 495101 : for (Int irow=0; irow!=s.nRows(); irow++) {
854 468598 : if (s.flagRow()(irow))
855 1653 : continue;
856 466945 : Int a1(s.antenna1()(irow));
857 466945 : Int a2(s.antenna2()(irow));
858 : // Not using irow
859 466945 : Matrix<Bool> flr = fl.xyPlane(irow);
860 466945 : if (!allTrue(flr)) {
861 436448 : activeAntennas.insert(a1);
862 436448 : activeAntennas.insert(a2);
863 : }
864 466945 : }
865 26503 : }
866 219 : if (prtlev > 2) {
867 0 : cout << "[FringeJones.cc::findRefAntWithData] refantlist " << refAntList << endl;
868 0 : cout << "[FringeJones.cc::findRefAntWithData] activeAntennas: ";
869 0 : std::copy(
870 : activeAntennas.begin(),
871 : activeAntennas.end(),
872 0 : std::ostream_iterator<Int>(std::cout, " ")
873 : );
874 0 : cout << endl;
875 : }
876 219 : Int refAnt = -1;
877 438 : for (Vector<Int>::ConstIteratorSTL a = refAntList.begin(); a != refAntList.end(); a++) {
878 219 : if (activeAntennas.find(*a) != activeAntennas.end()) {
879 219 : if (prtlev > 2)
880 0 : cout << "[FringeJones.cc::findRefAntWithData] We are choosing refant " << *a << endl;
881 219 : refAnt = *a;
882 219 : break;
883 : } else {
884 0 : if (prtlev > 2)
885 0 : cout << "[FringeJones.cc::findRefAntWithData] No data for refant " << *a << endl;
886 : }
887 219 : }
888 219 : return refAnt;
889 219 : }
890 :
891 :
892 : // Stolen from SolveDataBuffer
893 : void
894 219 : aggregateTimeCentroid(SDBList& sdbs, Int refAnt, std::map<Int, Double>& aggregateTime) {
895 : // Weighted average of SDBs' timeCentroids
896 219 : std::map<Int, Double> aggregateWeight;
897 26754 : for (Int i=0; i < sdbs.nSDB(); ++i) {
898 26535 : SolveDataBuffer& s = sdbs(i);
899 26535 : Vector<Double> wtvD;
900 : // Sum over correlations and channels to get a vector of weights for each row
901 53070 : Vector<Float> wtv(partialSums(s.weightSpectrum(), IPosition(2, 0, 1)));
902 26535 : wtvD.resize(wtv.nelements());
903 26535 : convertArray(wtvD, wtv);
904 495325 : for (Int j=0; j<s.nRows(); j++) {
905 468790 : Int a1 = s.antenna1()(j);
906 468790 : Int a2 = s.antenna2()(j);
907 : Int ant;
908 468790 : if (a1 == refAnt) { ant = a2; }
909 338749 : else if (a2 == refAnt) { ant = a1; }
910 338029 : else continue;
911 130761 : Double w = wtv(j);
912 130761 : aggregateTime[ant] += w*s.timeCentroid()(j);
913 130761 : aggregateWeight[ant] += w;
914 : }
915 26535 : }
916 1464 : for (auto it=aggregateTime.begin(); it!=aggregateTime.end(); ++it) {
917 1245 : Int a = it->first;
918 1245 : it->second /= aggregateWeight[a];
919 : }
920 :
921 219 : }
922 :
923 :
924 : void
925 0 : print_gsl_vector(gsl_vector *v)
926 : {
927 0 : const size_t n = v->size;
928 0 : for (size_t i=0; i!=n; i++) {
929 0 : cerr << gsl_vector_get(v, i) << " ";
930 0 : if (i>0 && (i % 4)==3) cerr << endl;
931 : }
932 0 : cerr << endl;
933 0 : }
934 :
935 : void
936 0 : print_max_gsl3(gsl_vector *v)
937 : {
938 0 : double phi_max = 0.0;
939 0 : double del_max = 0.0;
940 0 : double rat_max = 0.0;
941 :
942 0 : const size_t n = v->size;
943 0 : for (size_t i=0; i!=n/3; i++) {
944 0 : if (fabs(gsl_vector_get(v, 3*i+0)) > fabs(phi_max)) phi_max = gsl_vector_get(v, 3*i+0);
945 0 : if (fabs(gsl_vector_get(v, 3*i+1)) > fabs(del_max)) del_max = gsl_vector_get(v, 3*i+1);
946 0 : if (fabs(gsl_vector_get(v, 3*i+2)) > fabs(rat_max)) rat_max = gsl_vector_get(v, 3*i+2);
947 : }
948 0 : cerr << "phi_max " << phi_max << " del_max " << del_max << " rat_max " << rat_max << endl;
949 0 : }
950 :
951 :
952 :
953 : /*
954 : gsl-2.4/multilarge_nlinear/fdf.c defines gsl_multilarge_nlinear_driver,
955 : which I have butchered for my purposes here into
956 : not_gsl_multilarge_nlinear_driver(). We still iterate the nonlinear
957 : least squares solver until completion, but we adopt a convergence
958 : criterion copied from AIPS.
959 :
960 : Inputs: maxiter - maximum iterations to allow
961 : w - workspace
962 :
963 : Additionally I've removed the info parameter, and I may yet regret it.
964 :
965 : Originally:
966 : info - (output) info flag on why iteration terminated
967 : 1 = stopped due to small step size ||dx|
968 : 2 = stopped due to small gradient
969 : 3 = stopped due to small change in f
970 : GSL_ETOLX = ||dx|| has converged to within machine
971 : precision (and xtol is too small)
972 : GSL_ETOLG = ||g||_inf is smaller than machine
973 : precision (gtol is too small)
974 : GSL_ETOLF = change in ||f|| is smaller than machine
975 : precision (ftol is too small)
976 :
977 : Return:
978 : GSL_SUCCESS if converged
979 : GSL_MAXITER if maxiter exceeded without converging
980 : GSL_ENOPROG if no accepted step found on first iteration
981 : */
982 :
983 : int
984 356 : least_squares_inner_driver (const size_t maxiter,
985 : gsl_multilarge_nlinear_workspace * w,
986 : AuxParamBundle *bundle)
987 : {
988 : int status;
989 356 : size_t iter = 0;
990 : /* call user callback function prior to any iterations
991 : * with initial system state */
992 : Double s;
993 356 : Double last_s = 1.0e30;
994 356 : Bool converged = false;
995 : do {
996 2518 : status = gsl_multilarge_nlinear_iterate (w);
997 : /*
998 : * If the solver reports no progress on the first iteration,
999 : * then it didn't find a single step to reduce the
1000 : * cost function and more iterations won't help so return.
1001 : *
1002 : * If we get a no progress flag on subsequent iterations,
1003 : * it means we did find a good step in a previous iteration,
1004 : * so continue iterating since the solver has now reset
1005 : * mu to its initial value.
1006 : */
1007 2518 : if (status == GSL_ENOPROG && iter == 0) {
1008 0 : return GSL_EMAXITER;
1009 : }
1010 :
1011 2518 : Double fnorm = gsl_blas_dnrm2(w->f);
1012 2518 : s = 0.5 * fnorm * fnorm;
1013 : if (DEVDEBUG) {
1014 : cerr << "Iter: " << iter << " ";
1015 : cerr << "Parameters: " << endl;
1016 : print_gsl_vector(w->x);
1017 : print_max_gsl3(w->dx);
1018 : cerr << "1 - s/last_s=" << 1 - s/last_s << endl;
1019 : }
1020 2518 : ++iter;
1021 2518 : if ((1 - s/last_s < 5e-6) && (iter > 1)) converged = true;
1022 2518 : last_s = s;
1023 : /* old test for convergence:
1024 : status = not_gsl_multilarge_nlinear_test(xtol, gtol, ftol, info, w); */
1025 2518 : } while (!converged && iter < maxiter);
1026 : /*
1027 : * the following error codes mean that the solution has converged
1028 : * to within machine precision, so record the error code in info
1029 : * and return success
1030 : */
1031 356 : if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG)
1032 : {
1033 0 : status = GSL_SUCCESS;
1034 : }
1035 : /* check if max iterations reached */
1036 356 : if (iter >= maxiter && status != GSL_SUCCESS)
1037 356 : status = GSL_EMAXITER; return status;
1038 : } /* gsl_multilarge_nlinear_driver() */
1039 :
1040 :
1041 :
1042 :
1043 : void
1044 219 : least_squares_driver(SDBList& sdbs, Matrix<Float>& casa_param, Matrix<Bool>& casa_flags, Matrix<Float>& casa_snr, Int refant,
1045 : const std::map< Int, std::set<Int> >& activeAntennas, Int maxits, Vector<Bool> paramActive, LogIO& logSink) {
1046 : // The variable casa_param is the Casa calibration framework's RParam matrix; we transcribe our results into it only at the end.
1047 : // n below is number of variables,
1048 : // p is number of parameters
1049 :
1050 : // We could pass in an AuxParamBundle instead I guess?
1051 219 : AuxParamBundle bundle(sdbs, refant, activeAntennas, paramActive);
1052 575 : for (size_t icor=0; icor != bundle.get_num_corrs(); icor++) {
1053 356 : bundle.set_active_corr(icor);
1054 356 : if (bundle.get_num_antennas() == 0) {
1055 0 : logSink << "No antennas for correlation " << icor << "; not running least-squares solver." << LogIO::POST;
1056 0 : continue;
1057 : }
1058 356 : if (bundle.get_num_antennas() == 1) {
1059 0 : logSink << "No baselines for correlation " << icor << "; not running least-squares solver." << LogIO::POST;
1060 0 : continue;
1061 : }
1062 : // Four parameters for every antenna, with dispersion
1063 356 : size_t p = bundle.nParameters() * (bundle.get_num_antennas() - 1);
1064 : // We need to store complex visibilities in a real matrix so we
1065 : // just store real and imaginary components separately.
1066 356 : size_t n = 2 * bundle.get_actual_num_data_points();
1067 :
1068 : if (DEVDEBUG) {
1069 : cerr << "bundle.nParameters() " << bundle.nParameters()
1070 : << " bundle.get_num_antennas() " <<bundle.get_num_antennas() << endl;
1071 : cerr << "p " << p << " n " << n << endl;
1072 : }
1073 : // Parameters for the least-squares solver.
1074 : // param_tol sets roughly the number of decimal places accuracy you want in the answer;
1075 : // I feel that 3 is probably plenty for fringe fitting.
1076 : // param_tol is not used
1077 : //const double param_tol = 1.0e-3;
1078 : // gtol is not used
1079 : // const double gtol = pow(GSL_DBL_EPSILON, 1.0/3.0);
1080 : // ftol is not used
1081 : // const double ftol = 1.0e-20;
1082 356 : const size_t max_iter = maxits ;
1083 :
1084 356 : const gsl_multilarge_nlinear_type *T = gsl_multilarge_nlinear_trust;
1085 356 : gsl_multilarge_nlinear_parameters params = gsl_multilarge_nlinear_default_parameters();
1086 : // the Moré scaling is the best equipped to handle very different scales;
1087 : //it should be the best choice to accommodate dispersive terms of O(f)
1088 356 : params.scale = gsl_multilarge_nlinear_scale_more;
1089 356 : params.trs = gsl_multilarge_nlinear_trs_lm;
1090 356 : params.solver = gsl_multilarge_nlinear_solver_cholesky;
1091 356 : gsl_multilarge_nlinear_workspace *w = gsl_multilarge_nlinear_alloc(T, ¶ms, n, p);
1092 : gsl_multilarge_nlinear_fdf f;
1093 :
1094 356 : f.f = &expb_f;
1095 : /* Can't set to NULL for finite-difference Jacobian in multilarge case. */
1096 356 : f.df = &expb_df;
1097 356 : f.n = n; /* number of data points */
1098 356 : f.p = p; /* number of parameters */
1099 356 : f.params = &bundle;
1100 :
1101 :
1102 : // Our original param is a matrix of (3*nCorr, nElem).
1103 : // We have to transcribe it to a vector.
1104 :
1105 356 : gsl_vector *gp = gsl_vector_alloc(p);
1106 356 : gsl_vector_set_zero(gp);
1107 :
1108 : // We transcribe Casa parameters into gsl vector format, as required by the solver.
1109 2756 : for (size_t iant=0; iant != bundle.get_max_antenna_index()+1; iant++) {
1110 2400 : if (!bundle.isActive(iant)) {
1111 244 : continue;
1112 : }
1113 10780 : for (int di=0; di<4; di++) {
1114 8624 : Int param_ind = bundle.get_param_corr_param_index(iant, di);
1115 8624 : if (param_ind < 0) continue;
1116 5416 : gsl_vector_set(gp, param_ind, casa_param(4*icor + di, iant));
1117 : }
1118 : }
1119 356 : gsl_vector *gp_orig = gsl_vector_alloc(p);
1120 : // Keep a copy of original parameters
1121 356 : gsl_vector_memcpy (gp_orig, gp);
1122 : // initialise workspace
1123 356 : gsl_multilarge_nlinear_init(gp, &f, w);
1124 :
1125 : // compute initial residual norm */
1126 356 : gsl_vector *res_f = gsl_multilarge_nlinear_residual(w);
1127 :
1128 : int info;
1129 356 : int status = least_squares_inner_driver(max_iter, w, &bundle);
1130 356 : double chi1 = gsl_blas_dnrm2(res_f);
1131 :
1132 356 : gsl_vector_sub(gp_orig, w->x);
1133 356 : gsl_vector *res = gsl_multilarge_nlinear_position(w);
1134 :
1135 : // We transcribe values back from gsl_vector to the param matrix
1136 :
1137 356 : gsl_matrix *hess = gsl_matrix_alloc(p,p);
1138 356 : gsl_vector *snr_vector = gsl_vector_alloc(p);
1139 356 : gsl_matrix_set_zero(hess);
1140 356 : gsl_vector_set_zero(snr_vector);
1141 356 : expb_hess(gp, &bundle, hess, chi1*chi1, snr_vector, logSink);
1142 :
1143 : // Transcribe parameters back into CASA arrays
1144 2756 : for (size_t iant=0; iant != bundle.get_max_antenna_index()+1; iant++) {
1145 2400 : if (!bundle.isActive(iant)) continue;
1146 2156 : Int iparam = bundle.get_param_corr_param_index(iant, 0);
1147 2156 : if (iparam<0) {
1148 356 : continue;
1149 : }
1150 : if (DEVDEBUG) {
1151 : logSink << "iparam " << iparam << endl;
1152 : logSink << "Old values for ant " << iant << " correlation " << icor
1153 : << " delay " << casa_param(4*icor + 1, iant) << " ns "
1154 : << " rate " << casa_param(4*icor + 2, iant)
1155 : << " angle " << casa_param(4*icor + 0, iant)
1156 : << endl;
1157 : logSink << "New values for ant " << iant << " correlation " << icor << ":";
1158 : int i;
1159 : if ((i=bundle.get_param_corr_param_index(iant, 0))>=0) {
1160 : logSink << " Angle " << gsl_vector_get(res, i);
1161 : }
1162 : if ((i=bundle.get_param_corr_param_index(iant, 1))>=0) {
1163 : logSink << " delay " << gsl_vector_get(res, i) << " ns ";
1164 : }
1165 : if ((i=bundle.get_param_corr_param_index(iant, 2))>=0) {
1166 : logSink << " rate " << gsl_vector_get(res, i);
1167 : }
1168 : logSink << "." << LogIO::POST;
1169 : }
1170 1800 : if (status==GSL_SUCCESS || status==GSL_EMAXITER) {
1171 : // Current policy is to assume that exceeding max
1172 : // number of iterations is not a deal-breaker, leave it
1173 : // to SNR calculation to decide if the results are
1174 : // useful.
1175 9000 : for (size_t di=0; di<4; di++) {
1176 7200 : int i=bundle.get_param_corr_param_index(iant, di);
1177 7200 : int i0 = bundle.get_param_corr_param_index(iant, 0);
1178 7200 : if (i>=0) {
1179 5416 : casa_param(4*icor + di, iant) = gsl_vector_get(res, i);
1180 5416 : casa_snr(4*icor + di, iant) = gsl_vector_get(snr_vector, i0);
1181 : } else {
1182 1784 : casa_param(4*icor + di, iant) = 0.0;
1183 1784 : casa_snr(4*icor + di, iant) = 0.0;
1184 : }
1185 : }
1186 1800 : } else { // gsl solver failed; flag data
1187 0 : logSink << "Least-squares solver failed to converge; flagging" << endl;
1188 0 : for (size_t di=0; di<4; di++) {
1189 0 : casa_flags(4*icor + di, iant) = false;
1190 : }
1191 : }
1192 : }
1193 :
1194 : logSink << "Least squares complete for correlation " << icor
1195 356 : << " after " << gsl_multilarge_nlinear_niter(w) << " iterations." << LogIO::POST;
1196 :
1197 : // << "reason for stopping: " << ((info == 1) ? "small step size" : "small gradient") << endl
1198 : // << "initial |f(x)| = " << chi0 << endl
1199 : // << "final |f(x)| = " << chi1 << endl
1200 : // << "final step taken = " << diffsize
1201 :
1202 : if (DEVDEBUG) {
1203 : cerr << "casa_param " << casa_param << endl;
1204 :
1205 : switch (info) {
1206 : case 1:
1207 : logSink << "Small step size." << endl;
1208 : break;
1209 : case 2:
1210 : logSink << "Flatness." << endl;
1211 : }
1212 : logSink << LogIO::POST;
1213 : }
1214 356 : gsl_vector_free(gp);
1215 356 : gsl_vector_free(gp_orig);
1216 356 : gsl_vector_free(snr_vector);
1217 356 : gsl_matrix_free(hess);
1218 356 : gsl_multilarge_nlinear_free(w);
1219 : }
1220 219 : }
1221 :
1222 :
1223 :
1224 :
1225 : // **********************************************************
1226 : // CTRateAwareTimeInterp1 Implementations
1227 : //
1228 :
1229 320 : CTRateAwareTimeInterp1::CTRateAwareTimeInterp1(NewCalTable& ct,
1230 : const casacore::String& timetype,
1231 : casacore::Array<Float>& result,
1232 320 : casacore::Array<Bool>& rflag) :
1233 320 : CTTimeInterp1(ct,timetype,result,rflag)
1234 320 : {}
1235 :
1236 : // Destructor (nothing to do locally)
1237 640 : CTRateAwareTimeInterp1::~CTRateAwareTimeInterp1() {}
1238 :
1239 48600 : Bool CTRateAwareTimeInterp1::interpolate(Double newtime) {
1240 : // Call generic first
1241 48600 : if (CTTimeInterp1::interpolate(newtime)) {
1242 : // Only if generic yields new calibration
1243 : // NB: lastWasExact_=exact in generic
1244 20029 : applyPhaseRate(timeType().contains("nearest") || lastWasExact_);
1245 20029 : return true;
1246 : }
1247 : else
1248 : // No change
1249 28571 : return false;
1250 :
1251 : }
1252 :
1253 : // Do the phase rate math
1254 20029 : void CTRateAwareTimeInterp1::applyPhaseRate(Bool single)
1255 : {
1256 :
1257 20029 : Int ispw=mcols_p->spwId()(0); // should only be one (sliced ct_)!
1258 20029 : MSSpectralWindow msSpw(ct_.spectralWindow());
1259 20029 : MSSpWindowColumns msCol(msSpw);
1260 : //Vector<Double> refFreqs;
1261 : //msCol.refFrequency().getColumn(refFreqs,True);
1262 :
1263 20029 : Vector<Double> freqs;
1264 20029 : msCol.chanFreq().get(ispw,freqs,True); // should only be 1
1265 20029 : Double centroidFreq=freqs(0);
1266 :
1267 : // cout << "time = " << (currTime_ - timeRef_) << endl;
1268 :
1269 20029 : if (single) {
1270 1083 : for (Int ipol=0;ipol<2;ipol++) {
1271 722 : Double dtime=(currTime_-timeRef_)-timelist_(currIdx_);
1272 722 : Double phase=result_(IPosition(2,ipol*4,0));
1273 722 : Double rate=result_(IPosition(2,ipol*4+2,0));
1274 722 : phase+=2.0*C::pi*rate*centroidFreq*dtime;
1275 722 : result_(IPosition(2,ipol*4,0))=phase;
1276 : }
1277 : } else {
1278 19668 : Vector<uInt> rows(2); indgen(rows); rows+=uInt(currIdx_);
1279 39336 : Cube<Float> r(mcols_p->fparamArray("",rows));
1280 :
1281 19668 : Vector<Double> dtime(2);
1282 19668 : dtime(0)=(currTime_-timeRef_)-timelist_(currIdx_);
1283 19668 : dtime(1)=(currTime_-timeRef_)-timelist_(currIdx_+1);
1284 19668 : Double wt=dtime(1) / (dtime(1)-dtime(0));
1285 :
1286 :
1287 59004 : for (Int ipol=0;ipol<2;ipol++) {
1288 39336 : Vector<Double> phase(2), rate(2);
1289 39336 : phase(0)=r.xyPlane(0)(IPosition(2,ipol*4,0));
1290 39336 : phase(1)=r.xyPlane(1)(IPosition(2,ipol*4,0));
1291 39336 : rate(0)=r.xyPlane(0)(IPosition(2,ipol*4+2,0));
1292 39336 : rate(1)=r.xyPlane(1)(IPosition(2,ipol*4+2,0));
1293 :
1294 39336 : phase(0)+=2.0*C::pi*rate(0)*centroidFreq*dtime(0);
1295 39336 : phase(1)+=2.0*C::pi*rate(1)*centroidFreq*dtime(1);
1296 :
1297 39336 : Vector<Complex> ph(2);
1298 39336 : ph(0)=Complex(cos(phase(0)),sin(phase(0)));
1299 39336 : ph(1)=Complex(cos(phase(1)),sin(phase(1)));
1300 39336 : ph(0)=Float(wt)*ph(0) + Float(1.0-wt)*ph(1);
1301 39336 : result_(IPosition(2,ipol*4,0))=arg(ph(0));
1302 39336 : }
1303 19668 : }
1304 20029 : }
1305 :
1306 :
1307 :
1308 :
1309 : // **********************************************************
1310 : // FringeJones Implementations
1311 : //
1312 0 : FringeJones::FringeJones(VisSet& vs) :
1313 : VisCal(vs), // virtual base
1314 : VisMueller(vs), // virtual base
1315 0 : GJones(vs) // immediate parent
1316 : {
1317 0 : if (prtlev()>2) cout << "FringeJones::FringeJones(vs)" << endl;
1318 0 : }
1319 :
1320 0 : FringeJones::FringeJones(String msname,Int MSnAnt,Int MSnSpw) :
1321 : VisCal(msname,MSnAnt,MSnSpw), // virtual base
1322 : VisMueller(msname,MSnAnt,MSnSpw), // virtual base
1323 0 : GJones(msname,MSnAnt,MSnSpw) // immediate parent
1324 : {
1325 0 : if (prtlev()>2) cout << "FringeJones::FringeJones(msname,MSnAnt,MSnSpw)" << endl;
1326 0 : }
1327 :
1328 34 : FringeJones::FringeJones(const MSMetaInfoForCal& msmc) :
1329 : VisCal(msmc), // virtual base
1330 : VisMueller(msmc), // virtual base
1331 34 : GJones(msmc) // immediate parent
1332 : {
1333 34 : if (prtlev()>2) cout << "FringeJones::FringeJones(msmc)" << endl;
1334 34 : }
1335 :
1336 0 : FringeJones::FringeJones(Int nAnt) :
1337 : VisCal(nAnt),
1338 : VisMueller(nAnt),
1339 0 : GJones(nAnt)
1340 : {
1341 0 : if (prtlev()>2) cout << "FringeJones::FringeJones(nAnt)" << endl;
1342 0 : }
1343 :
1344 68 : FringeJones::~FringeJones() {
1345 34 : if (prtlev()>2) cout << "FringeJones::~FringeJones()" << endl;
1346 68 : }
1347 :
1348 6 : void FringeJones::setApply(const Record& apply) {
1349 : // Call parent to do conventional things
1350 6 : GJones::setApply(apply);
1351 :
1352 6 : if (calWt())
1353 5 : logSink() << " (" << this->typeName() << ": Enforcing calWt()=false for phase/delay-like terms)" << LogIO::POST;
1354 :
1355 : // Enforce calWt() = false for delays
1356 6 : calWt()=false;
1357 :
1358 : /*
1359 : // Extract per-spw ref Freq for phase(delay) calculation
1360 : // from the CalTable
1361 : // TBD: revise as per refFreq decisions
1362 : MSSpectralWindow msSpw(ct_->spectralWindow());
1363 : MSSpWindowColumns msCol(msSpw);
1364 : msCol.refFrequency().getColumn(KrefFreqs_,true);
1365 : KrefFreqs_/=1.0e9; // in GHz
1366 :
1367 : /// Re-assign KrefFreq_ according spwmap (if any)
1368 : if (spwMap().nelements()>0) {
1369 : Vector<Double> tmpfreqs;
1370 : tmpfreqs.assign(KrefFreqs_);
1371 : for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
1372 : if (spwMap()(ispw)>-1)
1373 : KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
1374 : }
1375 : */
1376 :
1377 : // Use the "physical" (centroid) frequency, per spw
1378 6 : MSSpectralWindow msSpw(ct_->spectralWindow());
1379 6 : MSSpWindowColumns msCol(msSpw);
1380 6 : Vector<Double> tmpfreqs;
1381 6 : Vector<Double> chanfreq;
1382 6 : tmpfreqs.resize(msSpw.nrow());
1383 30 : for (rownr_t ispw=0;ispw<msSpw.nrow();++ispw) {
1384 24 : if (ispw < msSpw.nrow()) {
1385 24 : msCol.chanFreq().get(ispw,chanfreq,true); // reshape, if nec.
1386 24 : Int nch=chanfreq.nelements();
1387 24 : tmpfreqs(ispw)=chanfreq(nch/2);
1388 : }
1389 : }
1390 :
1391 6 : KrefFreqs_.resize(nSpw()); KrefFreqs_.set(0.0);
1392 30 : for (rownr_t ispw=0;ispw<(rownr_t)nSpw();++ispw) {
1393 24 : if (ispw < tmpfreqs.nelements())
1394 24 : KrefFreqs_(ispw)=tmpfreqs(ispw);
1395 : }
1396 :
1397 : /// Re-assign KrefFreq_ according spwmap (if any)
1398 6 : if (spwMap().nelements()>0) {
1399 12 : for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
1400 6 : if (spwMap()(ispw)>-1 && ispw < tmpfreqs.nelements())
1401 0 : KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
1402 : }
1403 :
1404 6 : KrefFreqs_/=1.0e9; // in GHz
1405 6 : }
1406 :
1407 0 : void FringeJones::setApply() {
1408 : // This was omitted in copying KJones. It shouldn't have been.
1409 :
1410 : // Call parent to do conventional things
1411 0 : GJones::setApply();
1412 :
1413 : // Enforce calWt() = false for delays
1414 0 : calWt()=false;
1415 :
1416 : // Set the ref freqs to something usable
1417 0 : KrefFreqs_.resize(nSpw());
1418 0 : KrefFreqs_.set(5.0);
1419 :
1420 0 : }
1421 :
1422 :
1423 2 : void FringeJones::setCallib(const Record& callib,
1424 : const MeasurementSet& selms) {
1425 :
1426 : // Call parent to do conventional things
1427 2 : SolvableVisCal::setCallib(callib,selms);
1428 :
1429 : /*
1430 : if (calWt())
1431 : logSink() << " (" << this->typeName() << ": Enforcing calWt()=false for phase/delay-like terms)" << LogIO::POST;
1432 : */
1433 : // Enforce calWt() = false for delays
1434 2 : calWt()=false;
1435 :
1436 : /*
1437 : // Extract per-spw ref Freq for phase(delay) calculation
1438 : // from the CalTable
1439 : KrefFreqs_.assign(cpp_->refFreqIn());
1440 : KrefFreqs_/=1.0e9; // in GHz
1441 :
1442 : // Re-assign KrefFreq_ according spwmap (if any)
1443 : if (spwMap().nelements()>0) {
1444 : Vector<Double> tmpfreqs;
1445 : tmpfreqs.assign(KrefFreqs_);
1446 : for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
1447 : if (spwMap()(ispw)>-1)
1448 : KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
1449 : }
1450 : */
1451 :
1452 : // Use the "physical" (centroid) frequency, per spw
1453 2 : KrefFreqs_.resize(nSpw());
1454 10 : for (Int ispw=0;ispw<nSpw();++ispw) {
1455 8 : const Vector<Double>& f(cpp_->freqIn(ispw));
1456 8 : Int nf=f.nelements();
1457 8 : KrefFreqs_[ispw]=f[nf/2]; // center (usually this will be same as [0])
1458 : }
1459 2 : KrefFreqs_/=1.0e9; // In GHz
1460 :
1461 : // Re-assign KrefFreq_ according spwmap (if any)
1462 2 : if (spwMap().nelements()>0) {
1463 2 : Vector<Double> tmpfreqs;
1464 2 : tmpfreqs.assign(KrefFreqs_);
1465 4 : for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
1466 2 : if (spwMap()(ispw)>-1)
1467 0 : KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
1468 2 : }
1469 :
1470 2 : }
1471 :
1472 26 : void FringeJones::setSolve(const Record& solve) {
1473 :
1474 : // Call parent to do conventional things
1475 26 : GJones::setSolve(solve);
1476 26 : preavg() = -DBL_MAX;
1477 :
1478 : // refant isn't properly set until selfSolveOne. We set it to a
1479 : // known value here so that it can be checked in debugging code.
1480 26 : refant() = -1;
1481 26 : if (prtlev() > 2) {
1482 0 : cout << "Before GJones::setSolve" << endl
1483 0 : << "FringeJones::setSolve()" <<endl
1484 0 : << "FringeJones::refant() = "<< refant() <<endl
1485 0 : << "FringeJones::refantlist() = "<< refantlist() <<endl;
1486 : }
1487 26 : if (solve.isDefined("zerorates")) {
1488 26 : zeroRates() = solve.asBool("zerorates");
1489 : }
1490 26 : if (solve.isDefined("globalsolve")) {
1491 26 : globalSolve() = solve.asBool("globalsolve");
1492 : }
1493 26 : if (solve.isDefined("delaywindow")) {
1494 26 : Array<Double> dw = solve.asArrayDouble("delaywindow");
1495 26 : delayWindow() = dw;
1496 26 : } else {
1497 0 : cerr << "No delay window!" << endl;
1498 : }
1499 26 : if (solve.isDefined("ratewindow")) {
1500 26 : rateWindow() = solve.asArrayDouble("ratewindow");
1501 : } else {
1502 0 : cerr << "No rate window!" << endl;
1503 : }
1504 26 : if (solve.isDefined("niter")) {
1505 26 : maxits() = solve.asInt("niter");
1506 : }
1507 26 : if (solve.isDefined("paramactive")) {
1508 26 : paramActive() = solve.asArrayBool("paramactive");
1509 : }
1510 26 : if (solve.isDefined("concatspws")) {
1511 26 : concatSPWs() = solve.asBool("concatspws");
1512 : }
1513 26 : if (solve.isDefined("corrcomb")) {
1514 : //cerr << "FringeJones::setsolve() Corrcomb is set! To:"
1515 : // << solve.asString("corrcomb")
1516 : // << endl;
1517 26 : corrcomb() = solve.asString("corrcomb");
1518 : }
1519 26 : }
1520 :
1521 : // Note: this was previously omitted
1522 0 : void FringeJones::specify(const Record& specify) {
1523 :
1524 0 : return SolvableVisCal::specify(specify);
1525 :
1526 : }
1527 :
1528 4680 : void FringeJones::calcAllJones() {
1529 :
1530 4680 : if (prtlev()>6) cout << " FringeJones::calcAllJones()" << endl;
1531 :
1532 : // Should handle OK flags in this method, and only
1533 : // do Jones calc if OK
1534 :
1535 4680 : Vector<Complex> oneJones;
1536 4680 : Vector<Bool> oneJOK;
1537 4680 : Vector<Float> onePar;
1538 4680 : Vector<Bool> onePOK;
1539 :
1540 4680 : ArrayIterator<Complex> Jiter(currJElem(),1);
1541 4680 : ArrayIterator<Bool> JOKiter(currJElemOK(),1);
1542 4680 : ArrayIterator<Float> Piter(currRPar(),1);
1543 4680 : ArrayIterator<Bool> POKiter(currParOK(),1);
1544 :
1545 : Double phase;
1546 :
1547 53280 : for (Int iant=0; iant<nAnt(); iant++) {
1548 48600 : onePar.reference(Piter.array());
1549 48600 : onePOK.reference(POKiter.array());
1550 :
1551 48600 : Float fmin_ = min(currFreq());
1552 48600 : Float fmax = max(currFreq());
1553 1638360 : for (Int ich=0; ich<nChanMat(); ich++) {
1554 :
1555 1589760 : oneJones.reference(Jiter.array());
1556 1589760 : oneJOK.reference(JOKiter.array());
1557 :
1558 4769280 : for (Int ipar=0;ipar<nPar();ipar+=4) {
1559 3179520 : if (onePOK(ipar)) {
1560 967680 : Float freq = currFreq()(ich);
1561 967680 : Float k_disp = 1e-9*KDISPSCALE*C::_2pi*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
1562 :
1563 967680 : phase=onePar(ipar);
1564 967680 : phase+=2.0*C::pi*onePar(ipar+1)*
1565 967680 : (currFreq()(ich)-KrefFreqs_(currSpw()));
1566 967680 : phase+=2.0*C::pi*onePar(ipar+2)*KrefFreqs_(currSpw())*1e9*
1567 967680 : (currTime() - refTime());
1568 967680 : Float dph_d = onePar(ipar+3) * k_disp;
1569 : if (DEVDEBUG && 0) {
1570 : cerr << "fmin_ " << fmin_ << " fmax " << fmax << " k_disp " << k_disp
1571 : << " param " << onePar(ipar+3) << " dph_d " << dph_d << endl;
1572 : }
1573 967680 : phase+=dph_d;
1574 967680 : oneJones(ipar/4)=Complex(cos(phase),sin(phase));
1575 967680 : oneJOK(ipar/4)=True;
1576 : } else {
1577 2211840 : oneJOK(ipar/4)=False;
1578 : }
1579 : }
1580 : // Advance iterators
1581 1589760 : Jiter.next();
1582 1589760 : JOKiter.next();
1583 : }
1584 : // Step to next antenns's pars
1585 48600 : Piter.next();
1586 48600 : POKiter.next();
1587 : }
1588 4680 : }
1589 :
1590 :
1591 : void
1592 219 : FringeJones::calculateSNR(Int nCorr, DelayRateFFT& drf) {
1593 219 : Matrix<Float> sRP(solveRPar().nonDegenerate(1));
1594 219 : Matrix<Bool> sPok(solveParOK().nonDegenerate(1));
1595 219 : Matrix<Float> sSNR(solveParSNR().nonDegenerate(1));
1596 :
1597 575 : for (size_t icor=0; icor != (size_t) nCorr; icor++) {
1598 356 : const set<Int>& activeAntennas = drf.getActiveAntennasCorrelation(icor);
1599 3848 : for (Int iant=0; iant != nAnt(); iant++) {
1600 3492 : if (iant == refant()) {
1601 356 : Double maxsnr = 999.0;
1602 356 : sSNR(4*icor + 0, iant) = maxsnr;
1603 356 : sSNR(4*icor + 1, iant) = maxsnr;
1604 356 : sSNR(4*icor + 2, iant) = maxsnr;
1605 356 : sSNR(4*icor + 3, iant) = maxsnr;
1606 : }
1607 3136 : else if (activeAntennas.find(iant) != activeAntennas.end()) {
1608 1800 : Double delay = sRP(4*icor + 1, iant);
1609 1800 : Double rate = sRP(4*icor + 2, iant);
1610 : // Note that DelayRateFFT::snr is also used to calculate SNRs for the least square values!
1611 1800 : Float snrval = drf.snr(icor, iant, delay, rate);
1612 1800 : sSNR(4*icor + 0, iant) = snrval;
1613 1800 : sSNR(4*icor + 1, iant) = snrval;
1614 1800 : sSNR(4*icor + 2, iant) = snrval;
1615 1800 : sSNR(4*icor + 3, iant) = snrval;
1616 : } else {
1617 1336 : sPok(4*icor + 0, iant) = false;
1618 1336 : sPok(4*icor + 1, iant) = false;
1619 1336 : sPok(4*icor + 2, iant) = false;
1620 1336 : sPok(4*icor + 3, iant) = false;
1621 : }
1622 : }
1623 : }
1624 219 : }
1625 :
1626 :
1627 :
1628 :
1629 : void
1630 219 : FringeJones::selfSolveOne(SDBList& sdbs) {
1631 219 : solveRPar()=0.0;
1632 219 : solveParOK()=false;
1633 219 : solveParErr()=1.0; // Does nothing?
1634 : // Maybe we put refFreq, refTime stuff in here?
1635 219 : Vector<Double> myRefFreqs;
1636 : // Cannot assume we have a calibration table (ct_) in this method.
1637 : // MSSpectralWindow msSpw(ct_->spectralWindow());
1638 : /// MSSpWindowColumns spwCol(msSpw);
1639 : // spwCol.refFrequency().getColumn(myRefFreqs, true);
1640 : //Double ref_freq = myRefFreqs(currSpw());
1641 : //Double ref_freq = sdbs.freqs()(0);
1642 : if (DEVDEBUG) {
1643 : cerr << "Solving for fringes for spw=" << currSpw() << endl
1644 : << "Valid spws are " << freqMetaData_.validSpws() << endl;
1645 : }
1646 219 : Int spw = currSpw();
1647 219 : Double reffreq = freqMetaData_.freq(spw)(0);
1648 219 : Double t0 = sdbs(0).time()(0);
1649 219 : Double dt0 = refTime() - t0;
1650 : //Double df0 = ref_freq - sdbs.freqs()(0);
1651 : //Double df0 = 0;
1652 219 : Double df0 = reffreq - sdbs(0).freqs()(0); // global center to global edge
1653 :
1654 219 : logSink() << "Solving for fringes for spw=" << currSpw() << " at t="
1655 219 : << MVTime(refTime()/C::day).string(MVTime::YMD,7) << LogIO::POST;
1656 :
1657 219 : std::map<Int, Double> aggregateTime;
1658 : // Set the refant to the first choice that has data!
1659 219 : refant() = findRefAntWithData(sdbs, refantlist(), prtlev());
1660 219 : if (refant()<0) {
1661 0 : logSink() << "Antennas " << refantlist() << LogIO::POST;
1662 : logSink() << "dt0 " << dt0
1663 : << " nSDB() " << sdbs.nSDB()
1664 0 : << LogIO::POST;
1665 0 : throw(AipsError("No valid reference antenna supplied."));
1666 : }
1667 : else {
1668 219 : logSink() << "Using reference antenna " << refant() << LogIO::POST;
1669 : }
1670 219 : aggregateTimeCentroid(sdbs, refant(), aggregateTime);
1671 :
1672 : if (DEVDEBUG) {
1673 : std::cerr << "Weighted time centroids" << endl;
1674 : for (auto it=aggregateTime.begin(); it!=aggregateTime.end(); ++it)
1675 : std::cerr << it->first << " => " << it->second - t0 << std::endl;
1676 : }
1677 :
1678 : // We arrange that we can use either the DelayRateFFT based on concatenation
1679 : // or the new one that combines spectral windows after the FFT
1680 : if (DEVDEBUG) {
1681 : std::cerr << "Making a DelayRateFFT" << endl;
1682 : }
1683 219 : DelayRateFFT *drfp = DelayRateFFT::makeAChild(concatSPWs(), sdbs, refant(), delayWindow(), rateWindow());
1684 : // DelayRateFFT *drfp = new DelayRateFFTConcat(sdbs, refant(), delayWindow(), rateWindow());
1685 : // DelayRateFFT *drfp = new DelayRateFFTConcat(sdbs, refant(), delayWindow(), rateWindow());
1686 : if (DEVDEBUG) {
1687 : cerr << "Made a DelayRateFFT of some kind!" << endl;
1688 : }
1689 219 : drfp->FFT();
1690 219 : drfp->searchPeak();
1691 219 : Matrix<Float> sRP(solveRPar().nonDegenerate(1));
1692 219 : Matrix<Bool> sPok(solveParOK().nonDegenerate(1));
1693 219 : Matrix<Float> sSNR(solveParSNR().nonDegenerate(1));
1694 219 : logSink() << "sPok " << sPok.shape() << LogIO::POST;
1695 :
1696 : // Map from MS antenna number to index
1697 : // transcribe fft results to sRP
1698 219 : Int ncol = drfp->param().ncolumn();
1699 219 : Int nrow = drfp->param().nrow();
1700 : if (DEVDEBUG) {
1701 : std::cerr << "nrow " << nrow << ", ncol " << ncol << endl;
1702 : std::cerr << "drfp->flag() " << drfp->flag() << endl;
1703 : }
1704 :
1705 1709 : for (Int i=0; i!=ncol; i++) { // i==iant
1706 8690 : for (Int j=0; j!=nrow; j++) { // j is parameter number
1707 7200 : Int oj = (j>=3) ? j+1 : j;
1708 7200 : sRP(IPosition(2, oj, i)) = drfp->param()(IPosition(2, j, i));
1709 7200 : sPok(IPosition(2, oj, i)) = !(drfp->flag()(IPosition(2, j, i)));
1710 : }
1711 : // Our estimate for dispersion is zero, unconditionally, and we stand by it.
1712 1490 : sPok(IPosition(2, 3, i)) = true;
1713 1490 : if (nrow > 3)
1714 910 : sPok(IPosition(2, 7, i)) = true;
1715 : }
1716 219 : size_t nCorrOrig(sdbs(0).nCorrelations());
1717 219 : size_t nCorr = (nCorrOrig> 1 ? 2 : 1); // number of p-hands
1718 219 : calculateSNR(nCorr, *drfp);
1719 219 : set<Int> belowThreshold;
1720 219 : Float threshold = minSNR();
1721 :
1722 575 : for (size_t icor=0; icor != nCorr; icor++) {
1723 356 : const set<Int>& activeAntennas = drfp->getActiveAntennasCorrelation(icor);
1724 3848 : for (Int iant=0; iant != nAnt(); iant++) {
1725 3492 : if (iant != refant() && (activeAntennas.find(iant) != activeAntennas.end())) {
1726 1800 : Float s = sSNR(4*icor + 0, iant);
1727 : // Start the log message; finished below
1728 1800 : logSink() << "Antenna " << iant << " correlation " << icor << " has (FFT) SNR of " << s;
1729 1800 : if (s < threshold) {
1730 0 : belowThreshold.insert(iant);
1731 0 : logSink() << " below threshold (" << threshold << ")";
1732 : // Don't assume these will be flagged later; do it right away.
1733 : // (The least squares routine will eventually become optional.)
1734 0 : sPok(4*icor + 0, iant) = false;
1735 0 : sPok(4*icor + 1, iant) = false;
1736 0 : sPok(4*icor + 2, iant) = false;
1737 0 : sPok(4*icor + 3, iant) = false;
1738 : }
1739 : // Finish the log message
1740 1800 : logSink() << "." << LogIO::POST;
1741 : }
1742 : }
1743 : // We currently remove the antennas below SNR threshold from
1744 : // the object used to handle the FFT fringe search.
1745 356 : drfp->removeAntennasCorrelation(icor, belowThreshold);
1746 : if (DEVDEBUG) {
1747 : drfp->printActive();
1748 : }
1749 : }
1750 219 : if (globalSolve()) {
1751 219 : logSink() << "Starting least squares optimization." << LogIO::POST;
1752 : // Note that least_squares_driver is *not* a method of
1753 : // FringeJones so we pass everything in, including the logSink
1754 : // reference. Note also that sRP is passed by reference and
1755 : // altered in place.
1756 438 : least_squares_driver(sdbs, sRP, sPok, sSNR, refant(), drfp->getActiveAntennas(), maxits(),
1757 219 : paramActive(), logSink());
1758 : }
1759 : else {
1760 0 : logSink() << "Skipping least squares optimisation." << LogIO::POST;
1761 : }
1762 :
1763 : if (DEVDEBUG) {
1764 : cerr << "Ref time " << MVTime(refTime()/C::day).string(MVTime::YMD,7) << endl;
1765 : cerr << "df0 " << df0 << " dt0 " << dt0 << " reffreq*dt0 " << reffreq*dt0 << endl;
1766 : cerr << "reffreq " << reffreq << endl;
1767 : cerr << "df0 " << df0 << " dt0 " << dt0 << " reffreq*dt0 " << reffreq*dt0 << endl;
1768 : cerr << "sRP " << sRP << endl;
1769 : }
1770 :
1771 2395 : for (Int iant=0; iant != nAnt(); iant++) {
1772 5668 : for (size_t icor=0; icor != nCorr; icor++) {
1773 3492 : const set<Int>& activeAntennas = drfp->getActiveAntennasCorrelation(icor);
1774 3492 : if (activeAntennas.find(iant) == activeAntennas.end()) {
1775 1336 : continue;
1776 : }
1777 2156 : Double phi0 = sRP(4*icor + 0, iant);
1778 2156 : Double delay = sRP(4*icor + 1, iant);
1779 2156 : Double rate = sRP(4*icor + 2, iant);
1780 2156 : Double k_disp = sRP(4*icor + 3, iant);
1781 2156 : aggregateTime.find(iant);
1782 : // We assume the reference frequency for fringe fitting
1783 : // (which is NOT the one stored in the SPECTRAL_WINDOW
1784 : // table) is the left-hand edge of the frequency grid.
1785 2156 : Double delta1 = df0*delay/1e9;
1786 2156 : Double delta2 = reffreq*dt0*rate;
1787 2156 : Double delta3 = C::_2pi*(delta1+delta2);
1788 : Double dt;
1789 2156 : auto p = aggregateTime.find(iant);
1790 2156 : if (zeroRates() && p != aggregateTime.end()) {
1791 72 : dt = p->second - t0;
1792 : } else {
1793 2084 : dt = refTime() - t0;
1794 : }
1795 : if (DEVDEBUG) {
1796 : cerr << "Antenna " << iant << ": phi0 " << phi0 << " delay " << delay << " rate " << rate << " dt " << dt << endl
1797 : << "k_disp " << k_disp << endl
1798 : << "reffreq "<< reffreq << " Adding corrections for frequency (" << 360*delta1 << ")"
1799 : << " and time (" << 360*delta2 << ") degrees." << endl;
1800 : }
1801 2156 : sRP(4*icor + 0, iant) += delta3;
1802 : }
1803 : }
1804 :
1805 : // We can zero the rates here (if needed) whether we did least squares or not.
1806 219 : if (zeroRates()) {
1807 12 : logSink() << "Zeroing delay rates in calibration table." << LogIO::POST;
1808 :
1809 36 : for (size_t icor=0; icor != nCorr; icor++) {
1810 224 : for (Int iant=0; iant != nAnt(); iant++) {
1811 200 : sRP(4*icor + 2, iant) = 0.0;
1812 : }
1813 : }
1814 : }
1815 : // Copy the results to the second polarisation for combined pols
1816 219 : if (corrcomb()=="all") {
1817 42 : logSink() << "Correlations combined: Copying results to other correlation" << LogIO::POST;
1818 462 : for (Int iant=0; iant != nAnt(); iant++) {
1819 2100 : for (Int i=0; i !=4; i++) {
1820 1680 : sRP(4+i, iant) = sRP(i, iant);
1821 1680 : sPok(4+i, iant) = sPok(i, iant);
1822 1680 : sSNR(4+i, iant) = sSNR(i, iant);
1823 : }
1824 : }
1825 : }
1826 : if (DEVDEBUG) {
1827 : std::cerr << "sPok " << sPok << endl;
1828 : }
1829 219 : delete drfp;
1830 219 : }
1831 :
1832 : void
1833 0 : FringeJones::solveOneVB(const VisBuffer&) {
1834 0 : throw(AipsError("VisBuffer interface not supported!"));
1835 : }
1836 :
1837 :
1838 25 : void FringeJones::globalPostSolveTinker() {
1839 :
1840 : // Re-reference the phase, if requested
1841 25 : if (refantlist()(0)>-1) applyRefAnt();
1842 25 : }
1843 :
1844 26 : void FringeJones::applyRefAnt() {
1845 :
1846 : // TBD:
1847 : // 1. Synchronize refant changes on par axis
1848 : // 2. Implement minimum mean deviation algorithm
1849 :
1850 26 : if (refantlist()(0)<0)
1851 0 : throw(AipsError("No refant specified."));
1852 :
1853 26 : Int nUserRefant=refantlist().nelements();
1854 :
1855 : // Get the preferred refant names from the MS
1856 26 : String refantName(msmc().antennaName(refantlist()(0)));
1857 26 : if (nUserRefant>1) {
1858 0 : refantName+=" (";
1859 0 : for (Int i=1;i<nUserRefant;++i) {
1860 0 : refantName+=msmc().antennaName(refantlist()(i));
1861 0 : if (i<nUserRefant-1) refantName+=",";
1862 : }
1863 0 : refantName+=")";
1864 : }
1865 :
1866 : logSink() << "Applying refant: " << refantName
1867 26 : << " refantmode = " << refantmode();
1868 26 : if (refantmode()=="flex")
1869 26 : logSink() << " (hold alternate refants' phase constant) when refant flagged";
1870 26 : if (refantmode()=="strict")
1871 0 : logSink() << " (flag all antennas when refant flagged)";
1872 26 : logSink() << LogIO::POST;
1873 :
1874 : // Generate a prioritized refant choice list
1875 : // The first entry in this list is the user's primary refant,
1876 : // the second entry is the refant used on the previous interval,
1877 : // and the rest is a prioritized list of alternate refants,
1878 : // starting with the user's secondary (if provided) refants,
1879 : // followed by the rest of the array, in distance order. This
1880 : // makes the priorities correct all the time, and prevents
1881 : // a semi-stochastic alternation (by preferring the last-used
1882 : // alternate, even if nominally higher-priority refants become
1883 : // available)
1884 :
1885 :
1886 : // Extract antenna positions
1887 26 : Matrix<Double> xyz;
1888 26 : if (msName()!="<noms>") {
1889 26 : MeasurementSet ms(msName());
1890 26 : MSAntennaColumns msant(ms.antenna());
1891 26 : msant.position().getColumn(xyz);
1892 26 : }
1893 : else {
1894 : // TBD RO*
1895 0 : CTColumns ctcol(*ct_);
1896 0 : CTAntennaColumns& antcol(ctcol.antenna());
1897 0 : antcol.position().getColumn(xyz);
1898 0 : }
1899 :
1900 : // Calculate (squared) antenna distances, relative
1901 : // to last preferred antenna
1902 26 : Vector<Double> dist2(xyz.ncolumn(),0.0);
1903 104 : for (Int i=0;i<3;++i) {
1904 78 : Vector<Double> row=xyz.row(i);
1905 78 : row-=row(refantlist()(nUserRefant-1));
1906 78 : dist2+=square(row);
1907 78 : }
1908 : // Move preferred antennas to a large distance
1909 52 : for (Int i=0;i<nUserRefant;++i)
1910 26 : dist2(refantlist()(i))=DBL_MAX;
1911 :
1912 : // Generated sorted index
1913 26 : Vector<uInt> ord;
1914 26 : genSort(ord,dist2);
1915 :
1916 : // Assemble the whole choices list
1917 26 : Int nchoices=nUserRefant+1+ord.nelements();
1918 26 : Vector<Int> refantchoices(nchoices,0);
1919 52 : Vector<Int> r(refantchoices(IPosition(1,nUserRefant+1),IPosition(1,refantchoices.nelements()-1)));
1920 26 : convertArray(r,ord);
1921 :
1922 : // set first two to primary preferred refant
1923 26 : refantchoices(0)=refantchoices(1)=refantlist()(0);
1924 :
1925 : // set user's secondary refants (if any)
1926 26 : if (nUserRefant>1)
1927 0 : refantchoices(IPosition(1,2),IPosition(1,nUserRefant))=
1928 0 : refantlist()(IPosition(1,1),IPosition(1,nUserRefant-1));
1929 :
1930 : //cout << "refantchoices = " << refantchoices << endl;
1931 :
1932 :
1933 26 : if (refantmode()=="strict") {
1934 0 : nchoices=1;
1935 0 : refantchoices.resize(1,True);
1936 : }
1937 :
1938 26 : Vector<Int> nPol(nSpw(),2); // TBD:or 1, if data was single pol
1939 :
1940 26 : if (nPar()==8) {
1941 : // Verify that 2nd poln has unflagged solutions, PER SPW
1942 26 : ROCTMainColumns ctmc(*ct_);
1943 :
1944 26 : Block<String> cols(1);
1945 26 : cols[0]="SPECTRAL_WINDOW_ID";
1946 26 : CTIter ctiter(*ct_,cols);
1947 26 : Cube<Bool> fl;
1948 :
1949 116 : while (!ctiter.pastEnd()) {
1950 :
1951 90 : Int ispw=ctiter.thisSpw();
1952 90 : fl.assign(ctiter.flag());
1953 :
1954 90 : IPosition blc(3,0,0,0), trc(fl.shape());
1955 90 : trc-=1; blc(0)=nPar()/2;
1956 :
1957 : // cout << "ispw = " << ispw << " nfalse(fl(1,:,:)) = " << nfalse(fl(blc,trc)) << endl;
1958 :
1959 : // If there are no unflagged solutions in 2nd pol,
1960 : // avoid it in refant calculations
1961 90 : if (nfalse(fl(blc,trc))==0)
1962 40 : nPol(ispw)=1;
1963 :
1964 90 : ctiter.next();
1965 90 : }
1966 26 : }
1967 : // cout << "nPol = " << nPol << endl;
1968 :
1969 26 : Bool usedaltrefant(false);
1970 26 : Int currrefant(refantchoices(0)), lastrefant(-1);
1971 :
1972 26 : Block<String> cols(2);
1973 26 : cols[0]="SPECTRAL_WINDOW_ID";
1974 26 : cols[1]="TIME";
1975 26 : CTIter ctiter(*ct_,cols);
1976 :
1977 : // Arrays to hold per-timestamp solutions
1978 26 : Cube<Float> solA, solB;
1979 26 : Cube<Bool> flA, flB;
1980 26 : Vector<Int> ant1A, ant1B, ant2B;
1981 26 : Matrix<Complex> refPhsr; // the reference phasor [npol,nchan]
1982 26 : Int lastspw(-1);
1983 26 : Bool first(true);
1984 253 : while (!ctiter.pastEnd()) {
1985 227 : Int ispw=ctiter.thisSpw();
1986 227 : if (ispw!=lastspw) first=true; // spw changed, start over
1987 :
1988 : // Read in the current sol, fl, ant1:
1989 227 : solB.assign(ctiter.fparam());
1990 227 : flB.assign(ctiter.flag());
1991 227 : ant1B.assign(ctiter.antenna1());
1992 227 : ant2B.assign(ctiter.antenna2());
1993 :
1994 : // First time thru, 'previous' solution same as 'current'
1995 227 : if (first) {
1996 90 : solA.reference(solB);
1997 90 : flA.reference(flB);
1998 90 : ant1A.reference(ant1B);
1999 : }
2000 227 : IPosition shB(solB.shape());
2001 227 : IPosition shA(solA.shape());
2002 :
2003 : // Find a good refant at this time
2004 : // A good refant is one that is unflagged in all polarizations
2005 : // in the current(B) and previous(A) intervals (so they can be connected)
2006 227 : Int irefA(0),irefB(0); // index on 3rd axis of solution arrays
2007 227 : Int ichoice(0); // index in refantchoicelist
2008 227 : Bool found(false);
2009 227 : IPosition blcA(3,0,0,0),trcA(shA),blcB(3,0,0,0),trcB(shB);
2010 227 : trcA-=1; trcA(0)=trcA(2)=0;
2011 227 : trcB-=1; trcB(0)=trcB(2)=0;
2012 227 : ichoice=0;
2013 454 : while (!found && ichoice<nchoices) {
2014 : // Find index of current refant choice
2015 227 : irefA=irefB=0;
2016 255 : while (ant1A(irefA)!=refantchoices(ichoice) && irefA<shA(2)) ++irefA;
2017 255 : while (ant1B(irefB)!=refantchoices(ichoice) && irefB<shB(2)) ++irefB;
2018 :
2019 227 : if (irefA<shA(2) && irefB<shB(2)) {
2020 :
2021 : // cout << " Trial irefA,irefB: " << irefA << "," << irefB
2022 : // << " Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;
2023 :
2024 227 : blcA(2)=trcA(2)=irefA;
2025 227 : blcB(2)=trcB(2)=irefB;
2026 227 : found=true; // maybe
2027 641 : for (Int ipol=0;ipol<nPol(ispw);++ipol) {
2028 414 : blcA(0)=blcB(0)=ipol*nPar()/2;
2029 414 : trcA(0)=trcB(0)=(ipol+1)*nPar()/2-1;
2030 414 : found &= (nfalse(flA(blcA,trcA))>0); // previous interval
2031 414 : found &= (nfalse(flB(blcB,trcB))>0); // current interval
2032 : }
2033 : }
2034 : else
2035 : // irefA or irefB out-of-range
2036 0 : found=false; // Just to be sure
2037 :
2038 227 : if (!found) ++ichoice; // try next choice next round
2039 :
2040 : }
2041 :
2042 227 : if (found) {
2043 : // at this point, irefA/irefB point to a good refant
2044 :
2045 : // Keep track
2046 227 : usedaltrefant|=(ichoice>0);
2047 227 : currrefant=refantchoices(ichoice);
2048 227 : refantchoices(1)=currrefant; // 2nd priorty next time
2049 :
2050 : // cout << " currrefant = " << currrefant << " (" << ichoice << ")" << endl;
2051 :
2052 : // cout << " Final irefA,irefB: " << irefA << "," << irefB
2053 : // << " Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;
2054 :
2055 :
2056 : // Only report if using an alternate refant
2057 227 : if (currrefant!=lastrefant && ichoice>0) {
2058 : logSink()
2059 : << "At "
2060 0 : << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2061 : << " ("
2062 : << "Spw=" << ctiter.thisSpw()
2063 : << ", Fld=" << ctiter.thisField()
2064 : << ")"
2065 0 : << ", using refant " << msmc().antennaName(currrefant)
2066 : << " (id=" << currrefant
2067 : << ")" << " (alternate)"
2068 0 : << LogIO::POST;
2069 : }
2070 :
2071 : // Form reference phasor [nPar,nChan]
2072 227 : Matrix<Float> rA,rB;
2073 227 : Matrix<Bool> rflA,rflB;
2074 227 : rB.assign(solB.xyPlane(irefB));
2075 227 : rflB.assign(flB.xyPlane(irefB));
2076 :
2077 227 : if (!first) {
2078 : // Get and condition previous phasor for the current refant
2079 137 : rA.assign(solA.xyPlane(irefA));
2080 137 : rflA.assign(flA.xyPlane(irefA));
2081 137 : rB-=rA;
2082 :
2083 : // Accumulate flags
2084 137 : rflB&=rflA;
2085 : }
2086 :
2087 : // cout << " rB = " << rB << endl;
2088 : // cout << boolalpha << " rflB = " << rflB << endl;
2089 : // TBD: fillChanGaps?
2090 :
2091 : // Now apply reference phasor to all antennas
2092 227 : Matrix<Float> thissol;
2093 2491 : for (Int iant=0;iant<shB(2);++iant) {
2094 2264 : thissol.reference(solB.xyPlane(iant));
2095 2264 : thissol-=rB;
2096 : }
2097 :
2098 : // Set refant, so we can put it back
2099 227 : ant2B=currrefant;
2100 :
2101 : // put back referenced solutions
2102 227 : ctiter.setfparam(solB);
2103 227 : ctiter.setantenna2(ant2B);
2104 :
2105 : // Next time thru, solB is previous
2106 227 : solA.reference(solB);
2107 227 : flA.reference(flB);
2108 227 : ant1A.reference(ant1B);
2109 227 : solB.resize(); // (break references)
2110 227 : flB.resize();
2111 227 : ant1B.resize();
2112 :
2113 227 : lastrefant=currrefant;
2114 227 : first=false; // avoid first-pass stuff from now on
2115 :
2116 227 : } // found
2117 : else {
2118 : logSink()
2119 : << "At "
2120 0 : << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2121 : << " ("
2122 : << "Spw=" << ctiter.thisSpw()
2123 : << ", Fld=" << ctiter.thisField()
2124 : << ")"
2125 : << ", refant (id=" << currrefant
2126 : << ") was flagged; flagging all antennas strictly."
2127 0 : << LogIO::POST;
2128 : // Flag all solutions in this interval
2129 0 : flB.set(True);
2130 0 : ctiter.setflag(flB);
2131 0 : if (ichoice == nchoices){
2132 : logSink() << LogIO::WARN
2133 : << "From time: "
2134 0 : << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2135 : << "in Spw: "
2136 : << ctiter.thisSpw()
2137 : << " refant list exhausted, flagging all solutions"
2138 0 : << LogIO::POST;
2139 : }
2140 : }
2141 :
2142 : // advance to the next interval
2143 227 : lastspw=ispw;
2144 227 : ctiter.next();
2145 227 : }
2146 :
2147 26 : if (usedaltrefant)
2148 : logSink() << LogIO::NORMAL
2149 : << " NB: An alternate refant was used at least once to maintain" << endl
2150 : << " phase continuity where the user's preferred refant drops out." << endl
2151 : << " Alternate refants are held constant in phase (_not_ zeroed)" << endl
2152 : << " during these periods, and the preferred refant may return at" << endl
2153 : << " a non-zero phase. This is generally harmless."
2154 0 : << LogIO::POST;
2155 :
2156 52 : return;
2157 :
2158 26 : }
2159 :
2160 : } //# NAMESPACE CASA - END
2161 :
2162 :
2163 : /*
2164 : Example of use in at Python level:
2165 :
2166 : fringefit(vis="n14c2.ms", caltable="fail.fj", field="",spw="1",intent="",
2167 : selectdata=True, timerange="", antenna="", scan="5", observation="",
2168 : msselect="", solint="inf", refant="EF", minsnr=3.0, append=False,
2169 : gaintable=['n14c2.gcal'], parang=False)
2170 : */
2171 :
|
