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 225 : AuxParamBundle(SDBList& sdbs_, size_t refant, const std::map< Int, std::set<Int> >& activeAntennas_, Vector<Bool> paramActive) :
107 225 : sdbs(sdbs_),
108 225 : nCalls(0),
109 225 : refant(refant),
110 225 : nCorrelations(sdbs.nCorrelations() > 1 ? 2 : 1),
111 225 : corrStep(sdbs.nCorrelations() > 2 ? 3 : 1),
112 225 : activeAntennas(activeAntennas_),
113 225 : parameterFlags(),
114 225 : parameterMap(),
115 225 : activeCorr(-1)
116 : // corrStep(3)
117 : {
118 225 : Int last_index = sdbs.nSDB() - 1 ;
119 225 : t0 = sdbs(0).time()(0);
120 225 : Double tlast = sdbs(last_index).time()(0);
121 225 : reftime = 0.5*(t0 + tlast);
122 :
123 225 : parameterFlags = paramActive.tovector();
124 225 : Int j = 0; // The CASA parameter index (0=peculiar phase, 1=delay, 2=rate, 3=dispersive)
125 225 : Int i = 0; // the Least Squares parameter vector index, depending on what's being solved for
126 1125 : for (auto p=parameterFlags.begin(); p!=parameterFlags.end(); p++) {
127 900 : if (*p) {
128 671 : parameterMap.insert(std::pair<Int, Int>(j, i));
129 671 : i++;
130 : }
131 900 : j++;
132 : }
133 225 : if (i==0) {
134 0 : throw(AipsError("No parameters specified!"));
135 : }
136 225 : nParams = i; // There's always at least one parameter!
137 : // cerr << "AuxParamBundle reftime " << reftime << " t0 " << t0 <<" dt " << tlast - t0 << endl;
138 225 : }
139 :
140 7862003326 : Int nParameters() {
141 7862003326 : return nParams;
142 : }
143 6572 : Double get_t0() {
144 6572 : return t0;
145 : }
146 : Double
147 : get_ref_time() {
148 : return reftime;
149 : }
150 : size_t
151 587 : get_num_corrs() {
152 : //return sdbs.nCorrelations() > 1 ? 2 : 1;
153 587 : return nCorrelations;
154 : }
155 : size_t
156 1086 : get_num_antennas() {
157 1086 : if (activeCorr < 0) {
158 0 : throw(AipsError("Correlation out of range."));
159 : }
160 1086 : std::set< Int > ants = activeAntennas.find(activeCorr)->second;
161 2172 : return (size_t) ants.size();
162 1086 : }
163 : size_t
164 5740 : get_max_antenna_index() {
165 5740 : if (activeCorr < 0) {
166 0 : throw(AipsError("Correlation out of range."));
167 : }
168 5740 : 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 nTotalChans = 0;
174 : Int nTotalRows = 0;
175 : for (Int i = 0; i != sdbs.nSDB(); i++) {
176 : SolveDataBuffer& s (sdbs(i));
177 : nTotalChans += s.nRows()*s.nChannels();
178 : }
179 : return nTotalChans;
180 : }
181 : Int
182 362 : get_actual_num_data_points() {
183 362 : Int nTotalRows = 0;
184 362 : Int nTotalChans = 0;
185 43358 : for (Int i = 0; i != sdbs.nSDB(); i++) {
186 42996 : SolveDataBuffer& s (sdbs(i));
187 854846 : for (Int irow=0; irow!=s.nRows(); irow++) {
188 811850 : if (s.flagRow()(irow)) continue;
189 809876 : nTotalChans += s.nChannels();
190 : }
191 : }
192 362 : return nTotalChans;
193 : }
194 : size_t
195 18042308 : get_data_corr_index(size_t icorr) {
196 18042308 : if (icorr > nCorrelations) {
197 0 : throw(AipsError("Correlation out of range."));
198 : }
199 18042308 : size_t dcorr = icorr * corrStep;
200 18042308 : return dcorr;
201 : }
202 : bool
203 37961824 : isActive(size_t iant) {
204 37961824 : std::set<Int> ants = activeAntennas.find(activeCorr)->second;
205 37961824 : if (iant == refant) return true;
206 33317408 : else return (ants.find(iant) != ants.end());
207 37961824 : }
208 : Int
209 11916533911 : get_param_corr_param_index(size_t iant0, size_t ipar) {
210 11916533911 : if (iant0 == refant) return -1;
211 10487183449 : int iant1 = antennaIndexMap[iant0];
212 10487183449 : if (iant1 > antennaIndexMap[refant]) {
213 10487183449 : iant1 -= 1;
214 : }
215 : int ipar1;
216 10487183449 : auto p = parameterMap.find(ipar);
217 10487183449 : if (p==parameterMap.end()) {
218 2625182386 : ipar1 = -1;
219 : }
220 : else {
221 7862001063 : ipar1 = (iant1 * nParameters()) + p->second;
222 : }
223 10487183449 : return ipar1;
224 : }
225 : size_t
226 18042308 : get_active_corr() {
227 18042308 : return activeCorr;
228 : }
229 : void
230 362 : set_active_corr(size_t icorr) {
231 362 : activeCorr = icorr;
232 362 : antennaIndexMap.clear();
233 362 : Int i = 0;
234 362 : std::set<Int>::iterator it;
235 362 : std::set<Int> ants = activeAntennas.find(activeCorr)->second;
236 2625 : for (it = ants.begin(); it != ants.end(); it++) {
237 2263 : antennaIndexMap[*it] = i++;
238 : }
239 362 : }
240 : };
241 :
242 :
243 :
244 : void
245 0 : print_baselines(std::set<std::pair< Int, Int > > baselines) {
246 0 : cerr << "Baselines encountered ";
247 0 : std::set<std::pair< Int, Int > >::iterator it;
248 0 : for (it=baselines.begin(); it != baselines.end(); ++it) {
249 0 : cerr << "(" << it->first << ", " << it->second << ") ";
250 : }
251 0 : cerr << endl;
252 0 : }
253 :
254 :
255 : int
256 3105 : expb_f(const gsl_vector *param, void *d, gsl_vector *f)
257 : {
258 3105 : AuxParamBundle *bundle = (AuxParamBundle *)d;
259 3105 : SDBList& sdbs = bundle->sdbs;
260 3105 : Double refTime = bundle->get_t0();
261 :
262 3105 : gsl_vector_set_zero(f);
263 : // Vector<Double> freqs = sdbs.freqs();
264 :
265 3105 : const Double reffreq0=sdbs(0).freqs()(0); // First freq in first SDB
266 :
267 3105 : size_t count = 0; // This is the master index.
268 :
269 3105 : Double sumwt = 0.0;
270 3105 : Double xi_squared = 0.0;
271 :
272 358095 : for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++) {
273 354990 : SolveDataBuffer& s (sdbs(ibuf));
274 354990 : if (!s.Ok()) continue;
275 :
276 354709 : const Vector<Double> freqs(s.freqs()); // This ibuf's freqs
277 354709 : Float fmin_ = min(freqs);
278 354709 : Float fmax = max(freqs);
279 :
280 354709 : const Cube<Complex>& v(s.visCubeCorrected());
281 354709 : const Cube<Bool>& fl(s.flagCube());
282 354709 : const Cube<Float>& weights(s.weightSpectrum());
283 :
284 9725805 : for (Int irow=0; irow!=s.nRows(); irow++) {
285 9371096 : if (s.flagRow()(irow)) continue;
286 :
287 9363793 : Int ant1(s.antenna1()(irow));
288 9363793 : Int ant2(s.antenna2()(irow));
289 9363793 : if (!bundle->isActive(ant1) || !bundle->isActive(ant2))
290 223832 : continue;
291 9139961 : if (ant1==ant2) continue;
292 :
293 : // VisBuffer.h seems to suggest that a vb.visCube may have shape
294 : // (nCorr(), nChannel(), nRow())
295 8661055 : size_t icorr0 = bundle->get_active_corr();
296 8661055 : size_t dcorr = bundle->get_data_corr_index(icorr0);
297 : // We also need to get the right parameters for this,
298 : // polarization (icorr is an encoding of the
299 : // polarization of the correlation products).
300 :
301 : Double phi0, tau, r, disp;
302 : {
303 : Int i;
304 : Double phi0_1, tau1, r1, disp1;
305 : Double phi0_2, tau2, r2, disp2;
306 :
307 8661055 : phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
308 8661055 : tau1 = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
309 8661055 : r1 = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
310 8661055 : disp1 = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
311 :
312 8661055 : phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
313 8661055 : tau2 = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
314 8661055 : r2 = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
315 8661055 : disp2 = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
316 :
317 8661055 : phi0 = phi0_2 - phi0_1;
318 8661055 : tau = tau2 - tau1;
319 8661055 : r = r2-r1;
320 8661055 : disp = disp2-disp1;
321 : }
322 :
323 180949719 : for (size_t ichan = 0; ichan != v.ncolumn(); ichan++) {
324 172288664 : if (fl(dcorr, ichan, irow)) continue;
325 :
326 171721550 : Float freq = freqs(ichan);
327 171721550 : Float k_disp = KDISPSCALE*(2.0*M_PI)*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
328 :
329 171721550 : Complex vis = v(dcorr, ichan, irow);
330 171721550 : Double w0 = weights(dcorr, ichan, irow);
331 : // FIXME: what should we use to scale the weights?
332 : // Double weightScale = norm(vis);
333 : // Double weightScale = abs(vis);
334 : // Double weightScale = 1;
335 : // Double weightScale = 1/sqrt(w0); // Actually AIPS 0, not AIPS 1!
336 : // Double weightScale = pow(w0, -0.75); // AIPS 2
337 : // Double weightScale = 1/w0; // AIPS 3
338 : // Double weightScale = norm(vis); // Casa 1, I guess
339 171721550 : Double w = sqrt(w0);
340 171721550 : sumwt += w*w;
341 171721550 : if (fabs(w) < FLT_EPSILON) continue;
342 : // We have to turn the delay back into seconds from nanoseconds.
343 : // Freq difference is in Hz, which comes out typically as 1e6 bands
344 : //Double wDf = (2.0*M_PI)*(freqs(ichan) - freqs(0))*1e-9;
345 171721550 : Double wDf = (2.0*M_PI)*(freqs(ichan) - reffreq0)*1e-9;
346 : //
347 171721550 : Double t1 = s.time()(0);
348 : // FIXME: Remind me why we *do* scale wDf with 1e-9
349 : // but do *not* do that with ref_freq?
350 : // I have a theory which is mine:
351 : // this is because tau is in nanoseconds.
352 : //Double ref_freq = freqs(0);
353 : //Double wDt = (2.0*M_PI)*(t1 - refTime) * ref_freq;
354 171721550 : Double wDt = (2.0*M_PI)*(t1 - refTime) * reffreq0;
355 :
356 171721550 : Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp);
357 171721550 : Double vtheta = arg(vis);
358 :
359 171721550 : Double c_r = w*(cos(mtheta) - cos(vtheta));
360 171721550 : Double c_i = w*(sin(mtheta) - sin(vtheta));
361 171721550 : gsl_vector_set(f, count, c_r);
362 171721550 : gsl_vector_set(f, count+1, c_i);
363 :
364 171721550 : count += 2;
365 171721550 : xi_squared += c_r*c_r + c_i*c_i;
366 : }
367 : }
368 354709 : }
369 : // cerr << "Residual xi-squared = " << xi_squared << endl;
370 3105 : return GSL_SUCCESS;
371 : }
372 :
373 :
374 : int
375 3105 : expb_df(CBLAS_TRANSPOSE_t TransJ, const gsl_vector* param, const gsl_vector *u, void *bundle_, gsl_vector *v, gsl_matrix *JTJ)
376 : {
377 :
378 : // param is the current vector for which we're finding the jacobian.
379 : // if TransJ is true, evaluate J^T u and store in v.
380 : // Also store J^T . J in lower half of JTJ.
381 3105 : std::set <std::pair < Int, Int> > baselines;
382 3105 : AuxParamBundle *bundle = (AuxParamBundle *)bundle_;
383 :
384 3105 : SDBList& sdbs = bundle->sdbs;
385 3105 : const Double reffreq0=sdbs(0).freqs()(0); // First freq in first SDB
386 :
387 3105 : size_t count = 0; // This is the master index.
388 :
389 3105 : gsl_vector_set_zero(v);
390 3105 : gsl_matrix_set_zero(JTJ);
391 :
392 3105 : Double refTime = bundle->get_t0();
393 :
394 358095 : for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++) {
395 354990 : SolveDataBuffer& s (sdbs(ibuf));
396 354990 : if (!s.Ok()) continue;
397 :
398 354709 : const Vector<Double>& freqs(s.freqs()); // This ibuf's freqs
399 354709 : Float fmin_ = min(freqs);
400 354709 : Float fmax = max(freqs);
401 354709 : const Cube<Complex>& vis(s.visCubeCorrected());
402 354709 : const Cube<Bool>& fl(s.flagCube());
403 354709 : const Cube<Float>& weights(s.weightSpectrum());
404 :
405 354709 : Double t1 = s.time()(0);
406 9725805 : for (Int irow=0; irow!=s.nRows(); irow++) {
407 10073834 : if (s.flagRow()(irow)) continue;
408 :
409 9363793 : Int ant1(s.antenna1()(irow));
410 9363793 : Int ant2(s.antenna2()(irow));
411 :
412 9363793 : if (ant1==ant2) continue;
413 8882023 : if (!bundle->isActive(ant1) || !bundle->isActive(ant2)) {
414 220968 : continue;
415 : }
416 :
417 : // VisBuffer.h seems to suggest that a vb.visCube may have shape
418 : // (nCorr(), nChannel(), nRow())
419 :
420 8661055 : size_t icorr0 = bundle->get_active_corr();
421 8661055 : size_t dcorr = bundle->get_data_corr_index(icorr0);
422 : // We also need to get the right parameters for this
423 : // polarization (icorr is an encoding of the
424 : // polarization of the correlation products).
425 :
426 : Double phi0, tau, r, disp;
427 : {
428 : Int i;
429 : Double phi0_1, tau1, r1, disp1;
430 : Double phi0_2, tau2, r2, disp2;
431 :
432 8661055 : phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
433 8661055 : tau1 = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
434 8661055 : r1 = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
435 8661055 : disp1 = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
436 :
437 8661055 : phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
438 8661055 : tau2 = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
439 8661055 : r2 = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
440 8661055 : disp2 = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
441 :
442 8661055 : phi0 = phi0_2 - phi0_1;
443 8661055 : tau = tau2 - tau1;
444 8661055 : r = r2-r1;
445 8661055 : disp = disp2-disp1;
446 : }
447 :
448 :
449 : //Double ref_freq = freqs(0);
450 : //Double wDt = (2.0*M_PI)*(t1 - refTime) * ref_freq;
451 8661055 : Double wDt = (2.0*M_PI)*(t1 - refTime) * reffreq0;
452 8661055 : bool found_data = false;
453 :
454 180949719 : for (size_t ichan = 0; ichan != vis.ncolumn(); ichan++) {
455 172288664 : if (fl(dcorr, ichan, irow)) continue;
456 171721550 : Double w0 = weights(dcorr, ichan, irow);
457 171721550 : Double w = sqrt(w0);
458 171721550 : if (fabs(w) < FLT_EPSILON) continue;
459 171721550 : found_data = true;
460 :
461 171721550 : Float freq = freqs(ichan);
462 171721550 : Float k_disp = KDISPSCALE*(2.0*M_PI)*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
463 :
464 : // Add a 1e-9 factor because tau parameter is in nanoseconds.
465 : //Double wDf = (2.0*M_PI)*(freqs(ichan) - freqs(0))*1e-9;
466 171721550 : Double wDf = (2.0*M_PI)*(freqs(ichan) - reffreq0)*1e-9;
467 : //
468 171721550 : Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp);
469 171721550 : Double ws = sin(mtheta);
470 171721550 : Double wc = cos(mtheta);
471 :
472 171721550 : Double p0 = 1.0;
473 171721550 : Double p1 = wDf;
474 171721550 : Double p2 = wDt;
475 171721550 : Double p3 = k_disp;
476 :
477 171721550 : Vector<Double> dterm2(4);
478 171721550 : dterm2(0) = -p0;
479 171721550 : dterm2(1) = -p1;
480 171721550 : dterm2(2) = -p2;
481 171721550 : dterm2(3) = -p3;
482 :
483 171721550 : Vector<Double> dterm1(4);
484 171721550 : dterm1(0) = p0;
485 171721550 : dterm1(1) = p1;
486 171721550 : dterm1(2) = p2;
487 171721550 : dterm1(3) = p3;
488 :
489 : /*
490 : What we want to express is just:
491 : J[count + 0, iparam2 + 0] = w*-ws*-1.0;
492 : J[count + 1, iparam2 + 0] = w*+wc*-1.0;
493 : J[count + 0, iparam2 + 1] = w*-ws*-wDf;
494 : J[count + 1, iparam2 + 1] = w*+wc*-wDf;
495 : J[count + 0, iparam2 + 2] = w*-ws*-wDt;
496 : J[count + 1, iparam2 + 2] = w*+wc*-wDt;
497 :
498 : But in the GSL multilarge framework we have to
499 : be ready to calculate either J*u for a given u
500 : or J^T*u, depending on the flag TransJ, and we also have to fill in the
501 :
502 : v[iparam + ...] = J[count + ..., iparam + ...] * u[iparam + ...]
503 :
504 : or
505 :
506 : v[iparam + ...] = J^T[iparam + ..., count + ...] * u[count + ...]
507 :
508 : <https://www.gnu.org/software/gsl/doc/html/nls.html#c.gsl_multifit_nlinear_default_parameters>
509 :
510 : "Additionally, the normal equations matrix J^T J should be stored in the lower half of JTJ."
511 :
512 : So we should also use
513 : JTJ[iparam + ..., iparam + ...] += J^T[iparam + ..., count + ...] J[count + ..., iparam + ...]
514 :
515 : */
516 171721550 : if (TransJ==CblasNoTrans) {
517 0 : for (Int di=0; di<4; di++) {
518 : Int i;
519 0 : if ((i = bundle->get_param_corr_param_index(ant2, di))>=0) {
520 0 : (*gsl_vector_ptr(v, count + 0)) += (w*-ws*dterm2(di)) * gsl_vector_get(u, i);
521 0 : (*gsl_vector_ptr(v, count + 1)) += (w*+wc*dterm2(di)) * gsl_vector_get(u, i);
522 : }
523 0 : if ((i = bundle->get_param_corr_param_index(ant1, di))>=0) {
524 0 : (*gsl_vector_ptr(v, count + 0)) += gsl_vector_get(u, i) * (w*-ws*dterm1(di));
525 0 : (*gsl_vector_ptr(v, count + 1)) += gsl_vector_get(u, i) * (w*+wc*dterm1(di));
526 : }
527 : }
528 : } else {
529 858607750 : for (Int di=0; di<4; di++) {
530 : Int i;
531 686886200 : if ((i = bundle->get_param_corr_param_index(ant2, di))>=0) {
532 514357436 : (*gsl_vector_ptr(v, i)) += (w*-ws*dterm2(di)) * gsl_vector_get(u, count + 0);
533 514357436 : (*gsl_vector_ptr(v, i)) += (w*+wc*dterm2(di)) * gsl_vector_get(u, count + 1);
534 : }
535 686886200 : if ((i = bundle->get_param_corr_param_index(ant1, di))>=0) {
536 385205432 : (*gsl_vector_ptr(v, i)) += gsl_vector_get(u, count + 0) * (w*-ws*dterm1(di));
537 385205432 : (*gsl_vector_ptr(v, i)) += gsl_vector_get(u, count + 1) * (w*+wc*dterm1(di));
538 : }
539 : }
540 : }
541 171721550 : if (JTJ) {
542 : Int i, j;
543 171721550 : Double wterm = (-ws) * (-ws) + (+wc) * (+wc);
544 171721550 : if (fabs(1-wterm) > 1e-15)
545 0 : throw AipsError("Insufficiently at one");
546 858607750 : for (Int di=0; di<4; di++) {
547 2404101700 : for (Int dj=0; dj<=di; dj++) {
548 2748972808 : if (((i = bundle->get_param_corr_param_index(ant2, di))>=0) &&
549 1031757308 : ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
550 1029231260 : (*gsl_matrix_ptr(JTJ, i, j)) += w0*dterm2(di)*dterm2(dj);
551 : }
552 2488667300 : if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
553 771451800 : ((j = bundle->get_param_corr_param_index(ant1, dj))>=0)) {
554 770609784 : (*gsl_matrix_ptr(JTJ, i, j)) += w0*dterm1(di)*dterm1(dj);
555 : }
556 : }
557 : }
558 : // iant1 != iant2, so we don't have to worry about collisions
559 858607750 : for (Int di=0; di<4; di++) {
560 3434431000 : for (Int dj=0; dj<4; dj++) {
561 : Int i0, j0;
562 4288366528 : if (((i0 = bundle->get_param_corr_param_index(ant1, di))>=0) &&
563 1540821728 : ((j0 = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
564 1156014136 : Int i1 = max(i0, j0);
565 1156014136 : Int j1 = min(i0, j0);
566 1156014136 : (*gsl_matrix_ptr(JTJ, i1, j1)) += w0*dterm2(di)*dterm1(dj);
567 : }
568 : }
569 : }
570 171721550 : count += 2;
571 : } // if JTJ
572 171721550 : } // loop over channels
573 8661055 : if (found_data) {
574 8661055 : std::pair<Int, Int> antpair = std::make_pair(ant1, ant2);
575 8661055 : bool newBaseline = (baselines.find(antpair) == baselines.end());
576 8661055 : if (newBaseline) {
577 87542 : baselines.insert(antpair);
578 : // cerr << "paramFlagging for antenna "<< ant1 << ": ";
579 : // for (size_t di=0; di<4; di++) {
580 : // cerr << (bundle->get_param_corr_param_index(ant1, di)>=0) << " ";
581 : // }
582 : // cerr << endl;
583 : // cerr << "indices for antenna "<< ant1 << ": ";
584 : // if (bundle->get_param_corr_param_index(ant1, 0) >= 0) {
585 : // for (size_t di=0; di<4; di++) {
586 : // cerr << bundle->get_param_corr_param_index(ant1, di) << " ";
587 : // }
588 : // cerr << endl;
589 : // }
590 : // cerr << "phi0 " << phi0 << " tau " << tau << " r " << r << endl;
591 : }
592 : }
593 : }
594 : }
595 : if (DEVDEBUG && 0) {
596 : print_baselines(baselines);
597 : cerr << "count " << count << endl;
598 : cerr << "v = ";
599 : for (size_t i=0; i!=v->size; i++) {
600 : cerr << gsl_vector_get(v, i) << " ";
601 : }
602 : cerr << endl;
603 : // if (JTJ) {
604 : // cerr <<"JTJ " << std::scientific << endl;
605 : // for (size_t i=0; i!=JTJ->size1; i++) {
606 : // for (size_t j=0; j!=JTJ->size2; j++) {
607 : // cerr << gsl_matrix_get(JTJ, i, j) << " ";
608 : // }
609 : // cerr << endl;
610 : // }
611 : // cerr << endl;
612 : // }
613 : }
614 3105 : return GSL_SUCCESS;
615 3105 : }
616 :
617 :
618 :
619 : int
620 362 : expb_hess(gsl_vector *param, AuxParamBundle *bundle, gsl_matrix *hess, Double xi_squared, gsl_vector *snr_vector, LogIO& logSink)
621 : {
622 : // We calculate the diagonal for the hessian as used by AIPS for
623 : // the signal to noise. The AIPS formulation, by Fred Schwab, is a
624 : // hand-rolled routine that solves a different problem to ours: by
625 : // using a triangular matrix for the Jacobian (requiring antenna i<
626 : // antenna j) the J^T * J term vanishes throughout and the Hessian
627 : // of the Xi^2 functional *only* includes the second-order
628 : // derivative terms, which are usually neglected.
629 : //
630 : // This is very clever, but it also means different covariance and
631 : // information matrices, and therefore a different SNR. Here we
632 : // use a generic least squares solver but we cheat slightly and use
633 : // the AIPS form for the Hessian and SNR calculations.
634 : //
635 : // FIXME: Is there any compelling reason to use gsl_vectors for
636 : // this, given that we're not really hooked in to the gsl
637 : // least squares framework by the time we do this?
638 :
639 362 : SDBList& sdbs = bundle->sdbs;
640 362 : Double refTime = bundle->get_t0();
641 :
642 : // Dimensions of (num_antennas); is the same dimension as
643 : // param vector here.
644 362 : gsl_matrix_set_zero(hess);
645 :
646 362 : const Double reffreq0=sdbs(0).freqs()(0); // First freq in first SDB
647 :
648 362 : size_t nobs = 0;
649 362 : Double sumwt = 0;
650 362 : size_t numpar = param->size;
651 :
652 43358 : for (Int ibuf=0; ibuf < sdbs.nSDB(); ibuf++)
653 : {
654 42996 : SolveDataBuffer& s (sdbs(ibuf));
655 42996 : if (!s.Ok()) continue;
656 :
657 42936 : const Vector<Double> freqs(s.freqs()); // This ibuf's freqs
658 42936 : Float fmin_ = min(freqs);
659 42936 : Float fmax = max(freqs);
660 :
661 42936 : const Cube<Complex>& v(s.visCubeCorrected());
662 42936 : const Cube<Bool>& fl(s.flagCube());
663 42936 : const Cube<Float>& weights(s.weightSpectrum());
664 :
665 :
666 854426 : for (Int irow=0; irow!=s.nRows(); irow++) {
667 811490 : if (s.flagRow()(irow)) continue;
668 :
669 809816 : Int ant1(s.antenna1()(irow));
670 809816 : Int ant2(s.antenna2()(irow));
671 809816 : if (!bundle->isActive(ant1) || !bundle->isActive(ant2))
672 11408 : continue;
673 798408 : if (ant1==ant2) continue;
674 :
675 : // VisBuffer.h seems to suggest that a vb.visCube may have shape
676 : // (nCorr(), nChannel(), nRow())
677 720198 : size_t icorr0 = bundle->get_active_corr();
678 720198 : size_t dcorr = bundle->get_data_corr_index(icorr0);
679 : // We also need to get the right parameters for this,
680 : // polarization (icorr is an encoding of the
681 : // polarization of the correlation products).
682 :
683 : Double phi0, tau, r, disp;
684 : {
685 : Int i;
686 : Double phi0_1, tau1, r1, disp1;
687 : Double phi0_2, tau2, r2, disp2;
688 :
689 720198 : phi0_1 = ((i = bundle->get_param_corr_param_index(ant1, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
690 720198 : tau1 = ((i = bundle->get_param_corr_param_index(ant1, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
691 720198 : r1 = ((i = bundle->get_param_corr_param_index(ant1, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
692 720198 : disp1 = ((i = bundle->get_param_corr_param_index(ant1, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
693 :
694 :
695 720198 : phi0_2 = ((i = bundle->get_param_corr_param_index(ant2, 0))>=0) ? gsl_vector_get(param, i) : 0.0;
696 720198 : tau2 = ((i = bundle->get_param_corr_param_index(ant2, 1))>=0) ? gsl_vector_get(param, i) : 0.0;
697 720198 : r2 = ((i = bundle->get_param_corr_param_index(ant2, 2))>=0) ? gsl_vector_get(param, i) : 0.0;
698 720198 : disp2 = ((i = bundle->get_param_corr_param_index(ant2, 3))>=0) ? gsl_vector_get(param, i) : 0.0;
699 :
700 720198 : phi0 = phi0_2 - phi0_1;
701 720198 : tau = tau2 - tau1;
702 720198 : r = r2-r1;
703 720198 : disp = disp2-disp1;
704 : }
705 :
706 12199726 : for (size_t ichan = 0; ichan != v.ncolumn(); ichan++) {
707 11479528 : if (fl(dcorr, ichan, irow)) continue;
708 11428500 : Complex vis = v(dcorr, ichan, irow);
709 : // Fixme: this isn't a square root.
710 11428500 : Double w0 = weights(dcorr, ichan, irow);
711 11428500 : sumwt += w0;
712 11428500 : Double w = w0;
713 11428500 : nobs++;
714 11428500 : if (fabs(w) < FLT_EPSILON) continue;
715 :
716 : // We have to turn the delay back into seconds from nanoseconds.
717 : // Freq difference is in Hz, which comes out typically as 1e6 bands
718 : //Double wDf = (2.0*M_PI)*(freqs(ichan) - freqs(0))*1e-9;
719 11428500 : Double wDf = (2.0*M_PI)*(freqs(ichan) - reffreq0)*1e-9;
720 : //
721 11428500 : Double t1 = s.time()(0);
722 :
723 : //Double ref_freq = freqs(0);
724 : //Double wDt = (2.0*M_PI)*(t1 - refTime) * ref_freq;
725 11428500 : Double wDt = (2.0*M_PI)*(t1 - refTime) * reffreq0;
726 :
727 11428500 : Float freq = freqs(ichan);
728 11428500 : Float k_disp = KDISPSCALE*(2.0*M_PI)*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
729 :
730 11428500 : Double mtheta = -(phi0 + tau*wDf + r*wDt + disp*k_disp);
731 11428500 : Double vtheta = arg(vis);
732 :
733 : // Hold on a minute though!
734 11428500 : Double cx = w*cos(vtheta - mtheta);
735 :
736 11428500 : Matrix<Double> dterm(4,4);
737 11428500 : dterm(0, 0) = cx;
738 11428500 : dterm(0, 1) = wDf*cx;
739 11428500 : dterm(0, 2) = wDt*cx;
740 11428500 : dterm(0, 3) = k_disp*dterm(0, 1);
741 11428500 : dterm(1, 1) = wDf*dterm(0, 1);
742 11428500 : dterm(1, 2) = wDt*dterm(0, 1);
743 11428500 : dterm(1, 3) = wDf*dterm(0, 3);
744 11428500 : dterm(2, 2) = wDt*dterm(1, 2);
745 11428500 : dterm(2, 3) = wDt*dterm(1, 3);
746 11428500 : dterm(3, 3) = k_disp*dterm(2, 3);
747 :
748 : // Symmetry terms:
749 11428500 : dterm(1, 0) = dterm(0, 1);
750 11428500 : dterm(2, 0) = dterm(0, 2);
751 11428500 : dterm(3, 0) = dterm(0, 3);
752 11428500 : dterm(2, 1) = dterm(1, 2);
753 11428500 : dterm(3, 1) = dterm(1, 3);
754 11428500 : dterm(3, 2) = dterm(2, 3);
755 :
756 :
757 :
758 57142500 : for (Int di=0; di<4; di++) {
759 228570000 : for (Int dj=0; dj<4; dj++) {
760 : Int i, j;
761 276510944 : if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
762 93654944 : ((j = bundle->get_param_corr_param_index(ant1, dj))>=0)) {
763 70316600 : *gsl_matrix_ptr(hess, i, j) += dterm(di, dj);
764 : }
765 : // Exactly the same logic, but with antenna2
766 319368000 : if (((i = bundle->get_param_corr_param_index(ant2, di))>=0) &&
767 136512000 : ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
768 102578832 : *gsl_matrix_ptr(hess, i, j) += dterm(di, dj);
769 : }
770 : }
771 : }
772 : // FIXME: Not just diagonal terms any more!
773 : // Note that some of these are not in the lower
774 : // triangular part, even though they are copied
775 : // faithfully from AIPS which thinks it is filling
776 : // a triangular matrix and handles symmetry
777 : // later. Unless I've missed something (again).
778 57142500 : for (Int di=0; di<4; di++) {
779 228570000 : for (Int dj=0; dj<4; dj++) {
780 : Int i, j;
781 276510944 : if (((i = bundle->get_param_corr_param_index(ant1, di))>=0) &&
782 93654944 : ((j = bundle->get_param_corr_param_index(ant2, dj))>=0)) {
783 70316600 : *gsl_matrix_ptr(hess, i, j) -= dterm(di, dj);
784 70316600 : *gsl_matrix_ptr(hess, j, i) -= dterm(dj, di);
785 : } // if
786 : } // dj
787 : } // di
788 11428500 : } // ichan
789 : } // irow
790 42936 : } // ibuff
791 : // s is more tricky: it is the xi^2 term from exp_f
792 :
793 362 : xi_squared = max(xi_squared, 1e-25);
794 :
795 : if (DEVDEBUG && 0) {
796 : cerr << "The matrix is" << endl;
797 : cerr << setprecision(3) << scientific;
798 : for (size_t i = 0; i != hess->size1; i++) {
799 : for (size_t j = 0; j < hess->size2; j++) {
800 : cerr << gsl_matrix_get(hess,i,j) << " ";
801 : }
802 : cerr << endl;
803 : }
804 : cerr << endl;
805 : // str.unsetf(cerr:floatfield);
806 : cerr << std::defaultfloat;
807 : }
808 : //
809 362 : size_t hsize = hess->size1;
810 : int signum;
811 362 : gsl_permutation *perm = gsl_permutation_alloc (hsize);
812 362 : gsl_matrix *lu = gsl_matrix_alloc(hsize, hsize);
813 362 : gsl_matrix *inv_hess = gsl_matrix_alloc(hsize, hsize);
814 :
815 362 : gsl_linalg_LU_decomp(hess, perm, &signum);
816 362 : Double det = gsl_linalg_LU_det(hess, signum);
817 362 : if ((fabs(det) < GSL_DBL_EPSILON) || std::isnan(det)) {
818 0 : logSink << "Hessian matrix singular (determinant=" << det << "); setting signal-to-noise ratio to zero." << LogIO::POST;
819 : // Singular matrix; fill snrs with zero.
820 0 : for (size_t i=0; i < hess->size1; i+=bundle->nParameters()) {
821 0 : Double snr = 0;
822 0 : gsl_vector_set(snr_vector, i, snr);
823 : }
824 : }
825 : else {
826 362 : gsl_linalg_LU_invert(hess, perm, inv_hess);
827 :
828 362 : Double sigma2 = xi_squared / (nobs - numpar) * nobs / sumwt;
829 : // cerr << "xi_squared " << xi_squared << " Nobs " << nobs << " sumwt " << sumwt << " sigma2 " << sigma2 << endl;
830 2263 : for (size_t i=0; i < hess->size1; i+=bundle->nParameters()) {
831 1901 : Double h = gsl_matrix_get(inv_hess, i, i);
832 1901 : Double snr0 = sqrt(sigma2*h*0.5);
833 1901 : snr0 = min(snr0, 9999.999);
834 1901 : Double snr = (snr0 > 1e-20) ? 1.0/snr0 : snr0;
835 1901 : gsl_vector_set(snr_vector, i, snr);
836 : }
837 : }
838 362 : gsl_matrix_free(lu);
839 362 : gsl_matrix_free(inv_hess);
840 362 : gsl_permutation_free(perm);
841 : // SNR[i], according to aips, is 1/sqrt(sigma2*hess(i1,i1)*0.5);
842 : // Note that in AIPS i1 ranges over 1..NANT
843 : // We use 1 as a success code.
844 362 : return 1;
845 : }
846 :
847 :
848 : Int
849 225 : findRefAntWithData(SDBList& sdbs, Vector<Int>& refAntList, Int prtlev) {
850 225 : std::set<Int> activeAntennas;
851 27744 : for (Int ibuf=0; ibuf != sdbs.nSDB(); ibuf++) {
852 27519 : SolveDataBuffer& s(sdbs(ibuf));
853 27519 : if (!s.Ok())
854 30 : continue;
855 27489 : Cube<Bool> fl = s.flagCube();
856 563779 : for (Int irow=0; irow!=s.nRows(); irow++) {
857 536290 : if (s.flagRow()(irow))
858 1653 : continue;
859 534637 : Int a1(s.antenna1()(irow));
860 534637 : Int a2(s.antenna2()(irow));
861 : // Not using irow
862 534637 : Matrix<Bool> flr = fl.xyPlane(irow);
863 534637 : if (!allTrue(flr)) {
864 504152 : activeAntennas.insert(a1);
865 504152 : activeAntennas.insert(a2);
866 : }
867 534637 : }
868 27489 : }
869 225 : if (prtlev > 2) {
870 0 : cout << "[FringeJones.cc::findRefAntWithData] refantlist " << refAntList << endl;
871 0 : cout << "[FringeJones.cc::findRefAntWithData] activeAntennas: ";
872 0 : std::copy(
873 : activeAntennas.begin(),
874 : activeAntennas.end(),
875 0 : std::ostream_iterator<Int>(std::cout, " ")
876 : );
877 0 : cout << endl;
878 : }
879 225 : Int refAnt = -1;
880 450 : for (Vector<Int>::ConstIteratorSTL a = refAntList.begin(); a != refAntList.end(); a++) {
881 225 : if (activeAntennas.find(*a) != activeAntennas.end()) {
882 225 : if (prtlev > 2)
883 0 : cout << "[FringeJones.cc::findRefAntWithData] We are choosing refant " << *a << endl;
884 225 : refAnt = *a;
885 225 : break;
886 : } else {
887 0 : if (prtlev > 2)
888 0 : cout << "[FringeJones.cc::findRefAntWithData] No data for refant " << *a << endl;
889 : }
890 225 : }
891 225 : return refAnt;
892 225 : }
893 :
894 :
895 : // Stolen from SolveDataBuffer
896 : void
897 225 : aggregateTimeCentroid(SDBList& sdbs, Int refAnt, std::map<Int, Double>& aggregateTime) {
898 : // Weighted average of SDBs' timeCentroids
899 225 : std::map<Int, Double> aggregateWeight;
900 27744 : for (Int i=0; i < sdbs.nSDB(); ++i) {
901 27519 : SolveDataBuffer& s = sdbs(i);
902 27519 : Vector<Double> wtvD;
903 : // Sum over correlations and channels to get a vector of weights for each row
904 55038 : Vector<Float> wtv(partialSums(s.weightSpectrum(), IPosition(2, 0, 1)));
905 27519 : wtvD.resize(wtv.nelements());
906 27519 : convertArray(wtvD, wtv);
907 563989 : for (Int j=0; j<s.nRows(); j++) {
908 536470 : Int a1 = s.antenna1()(j);
909 536470 : Int a2 = s.antenna2()(j);
910 : Int ant;
911 536470 : if (a1 == refAnt) { ant = a2; }
912 393973 : else if (a2 == refAnt) { ant = a1; }
913 393253 : else continue;
914 143217 : Double w = wtv(j);
915 143217 : aggregateTime[ant] += w*s.timeCentroid()(j);
916 143217 : aggregateWeight[ant] += w;
917 : }
918 27519 : }
919 1548 : for (auto it=aggregateTime.begin(); it!=aggregateTime.end(); ++it) {
920 1323 : Int a = it->first;
921 1323 : it->second /= aggregateWeight[a];
922 : }
923 :
924 225 : }
925 :
926 :
927 : void
928 0 : print_gsl_vector(gsl_vector *v)
929 : {
930 0 : const size_t n = v->size;
931 0 : for (size_t i=0; i!=n; i++) {
932 0 : cerr << gsl_vector_get(v, i) << " ";
933 0 : if (i>0 && (i % 4)==3) cerr << endl;
934 : }
935 0 : cerr << endl;
936 0 : }
937 :
938 : void
939 0 : print_max_gsl3(gsl_vector *v)
940 : {
941 0 : double phi_max = 0.0;
942 0 : double del_max = 0.0;
943 0 : double rat_max = 0.0;
944 :
945 0 : const size_t n = v->size;
946 0 : for (size_t i=0; i!=n/3; i++) {
947 0 : if (fabs(gsl_vector_get(v, 3*i+0)) > fabs(phi_max)) phi_max = gsl_vector_get(v, 3*i+0);
948 0 : if (fabs(gsl_vector_get(v, 3*i+1)) > fabs(del_max)) del_max = gsl_vector_get(v, 3*i+1);
949 0 : if (fabs(gsl_vector_get(v, 3*i+2)) > fabs(rat_max)) rat_max = gsl_vector_get(v, 3*i+2);
950 : }
951 0 : cerr << "phi_max " << phi_max << " del_max " << del_max << " rat_max " << rat_max << endl;
952 0 : }
953 :
954 :
955 :
956 : /*
957 : gsl-2.4/multilarge_nlinear/fdf.c defines gsl_multilarge_nlinear_driver,
958 : which I have butchered for my purposes here into
959 : not_gsl_multilarge_nlinear_driver(). We still iterate the nonlinear
960 : least squares solver until completion, but we adopt a convergence
961 : criterion copied from AIPS.
962 :
963 : Inputs: maxiter - maximum iterations to allow
964 : w - workspace
965 :
966 : Additionally I've removed the info parameter, and I may yet regret it.
967 :
968 : Originally:
969 : info - (output) info flag on why iteration terminated
970 : 1 = stopped due to small step size ||dx|
971 : 2 = stopped due to small gradient
972 : 3 = stopped due to small change in f
973 : GSL_ETOLX = ||dx|| has converged to within machine
974 : precision (and xtol is too small)
975 : GSL_ETOLG = ||g||_inf is smaller than machine
976 : precision (gtol is too small)
977 : GSL_ETOLF = change in ||f|| is smaller than machine
978 : precision (ftol is too small)
979 :
980 : Return:
981 : GSL_SUCCESS if converged
982 : GSL_MAXITER if maxiter exceeded without converging
983 : GSL_ENOPROG if no accepted step found on first iteration
984 : */
985 :
986 : int
987 362 : least_squares_inner_driver (const size_t maxiter,
988 : gsl_multilarge_nlinear_workspace * w,
989 : AuxParamBundle *bundle)
990 : {
991 : int status;
992 362 : size_t iter = 0;
993 : /* call user callback function prior to any iterations
994 : * with initial system state */
995 : Double s;
996 362 : Double last_s = 1.0e30;
997 362 : Bool converged = false;
998 : do {
999 2743 : status = gsl_multilarge_nlinear_iterate (w);
1000 : /*
1001 : * If the solver reports no progress on the first iteration,
1002 : * then it didn't find a single step to reduce the
1003 : * cost function and more iterations won't help so return.
1004 : *
1005 : * If we get a no progress flag on subsequent iterations,
1006 : * it means we did find a good step in a previous iteration,
1007 : * so continue iterating since the solver has now reset
1008 : * mu to its initial value.
1009 : */
1010 2743 : if (status == GSL_ENOPROG && iter == 0) {
1011 0 : return GSL_EMAXITER;
1012 : }
1013 :
1014 2743 : Double fnorm = gsl_blas_dnrm2(w->f);
1015 2743 : s = 0.5 * fnorm * fnorm;
1016 : if (DEVDEBUG) {
1017 : cerr << "Iter: " << iter << " ";
1018 : cerr << "Parameters: " << endl;
1019 : print_gsl_vector(w->x);
1020 : print_max_gsl3(w->dx);
1021 : cerr << "1 - s/last_s=" << 1 - s/last_s << endl;
1022 : }
1023 2743 : ++iter;
1024 2743 : if ((1 - s/last_s < 5e-6) && (iter > 1)) converged = true;
1025 2743 : last_s = s;
1026 : /* old test for convergence:
1027 : status = not_gsl_multilarge_nlinear_test(xtol, gtol, ftol, info, w); */
1028 2743 : } while (!converged && iter < maxiter);
1029 : /*
1030 : * the following error codes mean that the solution has converged
1031 : * to within machine precision, so record the error code in info
1032 : * and return success
1033 : */
1034 362 : if (status == GSL_ETOLF || status == GSL_ETOLX || status == GSL_ETOLG)
1035 : {
1036 0 : status = GSL_SUCCESS;
1037 : }
1038 : /* check if max iterations reached */
1039 362 : if (iter >= maxiter && status != GSL_SUCCESS)
1040 362 : status = GSL_EMAXITER; return status;
1041 : } /* gsl_multilarge_nlinear_driver() */
1042 :
1043 :
1044 :
1045 :
1046 : void
1047 225 : least_squares_driver(SDBList& sdbs, Matrix<Float>& casa_param, Matrix<Bool>& casa_flags, Matrix<Float>& casa_snr, Int refant,
1048 : const std::map< Int, std::set<Int> >& activeAntennas, Int maxits, Vector<Bool> paramActive, LogIO& logSink) {
1049 : // The variable casa_param is the Casa calibration framework's RParam matrix; we transcribe our results into it only at the end.
1050 : // n below is number of variables,
1051 : // p is number of parameters
1052 :
1053 : // We could pass in an AuxParamBundle instead I guess?
1054 225 : AuxParamBundle bundle(sdbs, refant, activeAntennas, paramActive);
1055 587 : for (size_t icor=0; icor != bundle.get_num_corrs(); icor++) {
1056 362 : bundle.set_active_corr(icor);
1057 362 : if (bundle.get_num_antennas() == 0) {
1058 0 : logSink << "No antennas for correlation " << icor << "; not running least-squares solver." << LogIO::POST;
1059 0 : continue;
1060 : }
1061 362 : if (bundle.get_num_antennas() == 1) {
1062 0 : logSink << "No baselines for correlation " << icor << "; not running least-squares solver." << LogIO::POST;
1063 0 : continue;
1064 : }
1065 : // Four parameters for every antenna, with dispersion
1066 362 : size_t p = bundle.nParameters() * (bundle.get_num_antennas() - 1);
1067 : // We need to store complex visibilities in a real matrix so we
1068 : // just store real and imaginary components separately.
1069 362 : size_t n = 2 * bundle.get_actual_num_data_points();
1070 :
1071 : if (DEVDEBUG) {
1072 : cerr << "bundle.nParameters() " << bundle.nParameters()
1073 : << " bundle.get_num_antennas() " <<bundle.get_num_antennas() << endl;
1074 : cerr << "p " << p << " n " << n << endl;
1075 : }
1076 : // Parameters for the least-squares solver.
1077 : // param_tol sets roughly the number of decimal places accuracy you want in the answer;
1078 : // I feel that 3 is probably plenty for fringe fitting.
1079 : // param_tol is not used
1080 : //const double param_tol = 1.0e-3;
1081 : // gtol is not used
1082 : // const double gtol = pow(GSL_DBL_EPSILON, 1.0/3.0);
1083 : // ftol is not used
1084 : // const double ftol = 1.0e-20;
1085 362 : const size_t max_iter = maxits ;
1086 :
1087 362 : const gsl_multilarge_nlinear_type *T = gsl_multilarge_nlinear_trust;
1088 362 : gsl_multilarge_nlinear_parameters params = gsl_multilarge_nlinear_default_parameters();
1089 : // the Moré scaling is the best equipped to handle very different scales;
1090 : //it should be the best choice to accommodate dispersive terms of O(f)
1091 362 : params.scale = gsl_multilarge_nlinear_scale_more;
1092 362 : params.trs = gsl_multilarge_nlinear_trs_lm;
1093 362 : params.solver = gsl_multilarge_nlinear_solver_cholesky;
1094 362 : gsl_multilarge_nlinear_workspace *w = gsl_multilarge_nlinear_alloc(T, ¶ms, n, p);
1095 : gsl_multilarge_nlinear_fdf f;
1096 :
1097 362 : f.f = &expb_f;
1098 : /* Can't set to NULL for finite-difference Jacobian in multilarge case. */
1099 362 : f.df = &expb_df;
1100 362 : f.n = n; /* number of data points */
1101 362 : f.p = p; /* number of parameters */
1102 362 : f.params = &bundle;
1103 :
1104 :
1105 : // Our original param is a matrix of (3*nCorr, nElem).
1106 : // We have to transcribe it to a vector.
1107 :
1108 362 : gsl_vector *gp = gsl_vector_alloc(p);
1109 362 : gsl_vector_set_zero(gp);
1110 :
1111 : // We transcribe Casa parameters into gsl vector format, as required by the solver.
1112 2870 : for (size_t iant=0; iant != bundle.get_max_antenna_index()+1; iant++) {
1113 2508 : if (!bundle.isActive(iant)) {
1114 245 : continue;
1115 : }
1116 11315 : for (int di=0; di<4; di++) {
1117 9052 : Int param_ind = bundle.get_param_corr_param_index(iant, di);
1118 9052 : if (param_ind < 0) continue;
1119 5719 : gsl_vector_set(gp, param_ind, casa_param(4*icor + di, iant));
1120 : }
1121 : }
1122 362 : gsl_vector *gp_orig = gsl_vector_alloc(p);
1123 : // Keep a copy of original parameters
1124 362 : gsl_vector_memcpy (gp_orig, gp);
1125 : // initialise workspace
1126 362 : gsl_multilarge_nlinear_init(gp, &f, w);
1127 :
1128 : // compute initial residual norm */
1129 362 : gsl_vector *res_f = gsl_multilarge_nlinear_residual(w);
1130 :
1131 : int info;
1132 362 : int status = least_squares_inner_driver(max_iter, w, &bundle);
1133 362 : double chi1 = gsl_blas_dnrm2(res_f);
1134 :
1135 362 : gsl_vector_sub(gp_orig, w->x);
1136 362 : gsl_vector *res = gsl_multilarge_nlinear_position(w);
1137 :
1138 : // We transcribe values back from gsl_vector to the param matrix
1139 :
1140 362 : gsl_matrix *hess = gsl_matrix_alloc(p,p);
1141 362 : gsl_vector *snr_vector = gsl_vector_alloc(p);
1142 362 : gsl_matrix_set_zero(hess);
1143 362 : gsl_vector_set_zero(snr_vector);
1144 362 : expb_hess(gp, &bundle, hess, chi1*chi1, snr_vector, logSink);
1145 :
1146 : // Transcribe parameters back into CASA arrays
1147 2870 : for (size_t iant=0; iant != bundle.get_max_antenna_index()+1; iant++) {
1148 2508 : if (!bundle.isActive(iant)) continue;
1149 2263 : Int iparam = bundle.get_param_corr_param_index(iant, 0);
1150 2263 : if (iparam<0) {
1151 362 : continue;
1152 : }
1153 : if (DEVDEBUG) {
1154 : logSink << "iparam " << iparam << endl;
1155 : logSink << "Old values for ant " << iant << " correlation " << icor
1156 : << " delay " << casa_param(4*icor + 1, iant) << " ns "
1157 : << " rate " << casa_param(4*icor + 2, iant)
1158 : << " angle " << casa_param(4*icor + 0, iant)
1159 : << endl;
1160 : logSink << "New values for ant " << iant << " correlation " << icor << ":";
1161 : int i;
1162 : if ((i=bundle.get_param_corr_param_index(iant, 0))>=0) {
1163 : logSink << " Angle " << gsl_vector_get(res, i);
1164 : }
1165 : if ((i=bundle.get_param_corr_param_index(iant, 1))>=0) {
1166 : logSink << " delay " << gsl_vector_get(res, i) << " ns ";
1167 : }
1168 : if ((i=bundle.get_param_corr_param_index(iant, 2))>=0) {
1169 : logSink << " rate " << gsl_vector_get(res, i);
1170 : }
1171 : logSink << "." << LogIO::POST;
1172 : }
1173 1901 : if (status==GSL_SUCCESS || status==GSL_EMAXITER) {
1174 : // Current policy is to assume that exceeding max
1175 : // number of iterations is not a deal-breaker, leave it
1176 : // to SNR calculation to decide if the results are
1177 : // useful.
1178 9505 : for (size_t di=0; di<4; di++) {
1179 7604 : int i=bundle.get_param_corr_param_index(iant, di);
1180 7604 : int i0 = bundle.get_param_corr_param_index(iant, 0);
1181 7604 : if (i>=0) {
1182 5719 : casa_param(4*icor + di, iant) = gsl_vector_get(res, i);
1183 5719 : casa_snr(4*icor + di, iant) = gsl_vector_get(snr_vector, i0);
1184 : } else {
1185 1885 : casa_param(4*icor + di, iant) = 0.0;
1186 1885 : casa_snr(4*icor + di, iant) = 0.0;
1187 : }
1188 : }
1189 1901 : } else { // gsl solver failed; flag data
1190 0 : logSink << "Least-squares solver failed to converge; flagging" << endl;
1191 0 : for (size_t di=0; di<4; di++) {
1192 0 : casa_flags(4*icor + di, iant) = false;
1193 : }
1194 : }
1195 : }
1196 :
1197 : logSink << "Least squares complete for correlation " << icor
1198 362 : << " after " << gsl_multilarge_nlinear_niter(w) << " iterations." << LogIO::POST;
1199 :
1200 : // << "reason for stopping: " << ((info == 1) ? "small step size" : "small gradient") << endl
1201 : // << "initial |f(x)| = " << chi0 << endl
1202 : // << "final |f(x)| = " << chi1 << endl
1203 : // << "final step taken = " << diffsize
1204 :
1205 : if (DEVDEBUG) {
1206 : cerr << "casa_param " << casa_param << endl;
1207 :
1208 : switch (info) {
1209 : case 1:
1210 : logSink << "Small step size." << endl;
1211 : break;
1212 : case 2:
1213 : logSink << "Flatness." << endl;
1214 : }
1215 : logSink << LogIO::POST;
1216 : }
1217 362 : gsl_vector_free(gp);
1218 362 : gsl_vector_free(gp_orig);
1219 362 : gsl_vector_free(snr_vector);
1220 362 : gsl_matrix_free(hess);
1221 362 : gsl_multilarge_nlinear_free(w);
1222 : }
1223 225 : }
1224 :
1225 : // **********************************************************
1226 : // CTRateAwareTimeInterp1 Implementations
1227 : //
1228 :
1229 276 : CTRateAwareTimeInterp1::CTRateAwareTimeInterp1(NewCalTable& ct,
1230 : const casacore::String& timetype,
1231 : casacore::Array<Float>& result,
1232 276 : casacore::Array<Bool>& rflag) :
1233 276 : CTTimeInterp1(ct,timetype,result,rflag)
1234 276 : {}
1235 :
1236 : // Destructor (nothing to do locally)
1237 552 : 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*M_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*M_PI*rate(0)*centroidFreq*dtime(0);
1295 39336 : phase(1)+=2.0*M_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 39 : FringeJones::FringeJones(const MSMetaInfoForCal& msmc) :
1329 : VisCal(msmc), // virtual base
1330 : VisMueller(msmc), // virtual base
1331 39 : GJones(msmc) // immediate parent
1332 : {
1333 39 : if (prtlev()>2) cout << "FringeJones::FringeJones(msmc)" << endl;
1334 39 : }
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 78 : FringeJones::~FringeJones() {
1345 39 : if (prtlev()>2) cout << "FringeJones::~FringeJones()" << endl;
1346 78 : }
1347 :
1348 5 : void FringeJones::setApply(const Record& apply) {
1349 : // Call parent to do conventional things
1350 5 : GJones::setApply(apply);
1351 :
1352 5 : 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 5 : 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 5 : MSSpectralWindow msSpw(ct_->spectralWindow());
1379 5 : MSSpWindowColumns msCol(msSpw);
1380 5 : Vector<Double> tmpfreqs;
1381 5 : Vector<Double> chanfreq;
1382 5 : tmpfreqs.resize(msSpw.nrow());
1383 25 : for (rownr_t ispw=0;ispw<msSpw.nrow();++ispw) {
1384 20 : if (ispw < msSpw.nrow()) {
1385 20 : msCol.chanFreq().get(ispw,chanfreq,true); // reshape, if nec.
1386 20 : Int nch=chanfreq.nelements();
1387 20 : tmpfreqs(ispw)=chanfreq(nch/2);
1388 : }
1389 : }
1390 :
1391 5 : KrefFreqs_.resize(nSpw()); KrefFreqs_.set(0.0);
1392 25 : for (rownr_t ispw=0;ispw<(rownr_t)nSpw();++ispw) {
1393 20 : if (ispw < tmpfreqs.nelements())
1394 20 : KrefFreqs_(ispw)=tmpfreqs(ispw);
1395 : }
1396 :
1397 : /// Re-assign KrefFreq_ according spwmap (if any)
1398 5 : if (spwMap().nelements()>0) {
1399 10 : for (uInt ispw=0;ispw<spwMap().nelements();++ispw)
1400 5 : if (spwMap()(ispw)>-1 && ispw < tmpfreqs.nelements())
1401 0 : KrefFreqs_(ispw)=tmpfreqs(spwMap()(ispw));
1402 : }
1403 :
1404 5 : KrefFreqs_/=1.0e9; // in GHz
1405 5 : }
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 32 : void FringeJones::setSolve(const Record& solve) {
1473 :
1474 : // Call parent to do conventional things
1475 32 : GJones::setSolve(solve);
1476 32 : 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 32 : refant() = -1;
1481 32 : 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 32 : if (solve.isDefined("zerorates")) {
1488 32 : zeroRates() = solve.asBool("zerorates");
1489 : }
1490 32 : if (solve.isDefined("globalsolve")) {
1491 32 : globalSolve() = solve.asBool("globalsolve");
1492 : }
1493 32 : if (solve.isDefined("delaywindow")) {
1494 32 : Array<Double> dw = solve.asArrayDouble("delaywindow");
1495 32 : delayWindow() = dw;
1496 32 : } else {
1497 0 : cerr << "No delay window!" << endl;
1498 : }
1499 32 : if (solve.isDefined("ratewindow")) {
1500 32 : rateWindow() = solve.asArrayDouble("ratewindow");
1501 : } else {
1502 0 : cerr << "No rate window!" << endl;
1503 : }
1504 32 : if (solve.isDefined("niter")) {
1505 32 : maxits() = solve.asInt("niter");
1506 : }
1507 32 : if (solve.isDefined("paramactive")) {
1508 32 : paramActive() = solve.asArrayBool("paramactive");
1509 : }
1510 32 : if (solve.isDefined("concatspws")) {
1511 32 : concatSPWs() = solve.asBool("concatspws");
1512 : }
1513 32 : if (solve.isDefined("corrcomb")) {
1514 : //cerr << "FringeJones::setsolve() Corrcomb is set! To:"
1515 : // << solve.asString("corrcomb")
1516 : // << endl;
1517 32 : corrcomb() = solve.asString("corrcomb");
1518 : }
1519 32 : }
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*(2.0*M_PI)*(1.0/freq + (freq-fmin_-fmax)/(fmin_*fmax));
1562 :
1563 967680 : phase=onePar(ipar);
1564 967680 : phase+=2.0*M_PI*onePar(ipar+1)*
1565 967680 : (currFreq()(ich)-KrefFreqs_(currSpw()));
1566 967680 : phase+=2.0*M_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 225 : FringeJones::calculateSNR(Int nCorr, DelayRateFFT& drf) {
1593 225 : Matrix<Float> sRP(solveRPar().nonDegenerate(1));
1594 225 : Matrix<Bool> sPok(solveParOK().nonDegenerate(1));
1595 225 : Matrix<Float> sSNR(solveParSNR().nonDegenerate(1));
1596 :
1597 587 : for (size_t icor=0; icor != (size_t) nCorr; icor++) {
1598 362 : const set<Int>& activeAntennas = drf.getActiveAntennasCorrelation(icor);
1599 3906 : for (Int iant=0; iant != nAnt(); iant++) {
1600 3544 : if (iant == refant()) {
1601 362 : Double maxsnr = 999.0;
1602 362 : sSNR(4*icor + 0, iant) = maxsnr;
1603 362 : sSNR(4*icor + 1, iant) = maxsnr;
1604 362 : sSNR(4*icor + 2, iant) = maxsnr;
1605 362 : sSNR(4*icor + 3, iant) = maxsnr;
1606 : }
1607 3182 : else if (activeAntennas.find(iant) != activeAntennas.end()) {
1608 1906 : Double delay = sRP(4*icor + 1, iant);
1609 1906 : Double rate = sRP(4*icor + 2, iant);
1610 : // Note that DelayRateFFT::snr is also used to calculate SNRs for the least square values!
1611 1906 : Float snrval = drf.snr(icor, iant, delay, rate);
1612 1906 : sSNR(4*icor + 0, iant) = snrval;
1613 1906 : sSNR(4*icor + 1, iant) = snrval;
1614 1906 : sSNR(4*icor + 2, iant) = snrval;
1615 1906 : sSNR(4*icor + 3, iant) = snrval;
1616 : } else {
1617 1276 : sPok(4*icor + 0, iant) = false;
1618 1276 : sPok(4*icor + 1, iant) = false;
1619 1276 : sPok(4*icor + 2, iant) = false;
1620 1276 : sPok(4*icor + 3, iant) = false;
1621 : }
1622 : }
1623 : }
1624 225 : }
1625 :
1626 :
1627 :
1628 :
1629 : void
1630 225 : FringeJones::selfSolveOne(SDBList& sdbs) {
1631 225 : solveRPar()=0.0;
1632 225 : solveParOK()=false;
1633 225 : solveParErr()=1.0; // Does nothing?
1634 : // Maybe we put refFreq, refTime stuff in here?
1635 225 : 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 225 : Int spw = currSpw();
1647 225 : Double reffreq = freqMetaData_.freq(spw)(0);
1648 225 : Double t0 = sdbs(0).time()(0);
1649 225 : Double dt0 = refTime() - t0;
1650 : //Double df0 = ref_freq - sdbs.freqs()(0);
1651 : //Double df0 = 0;
1652 225 : Double df0 = reffreq - sdbs(0).freqs()(0); // global center to global edge
1653 :
1654 225 : logSink() << "Solving for fringes for spw=" << currSpw() << " at t="
1655 225 : << MVTime(refTime()/C::day).string(MVTime::YMD,7) << LogIO::POST;
1656 :
1657 225 : std::map<Int, Double> aggregateTime;
1658 : // Set the refant to the first choice that has data!
1659 225 : refant() = findRefAntWithData(sdbs, refantlist(), prtlev());
1660 225 : 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 225 : logSink() << "Using reference antenna " << refant() << LogIO::POST;
1669 : }
1670 225 : 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 225 : 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 225 : drfp->FFT();
1690 225 : drfp->searchPeak();
1691 225 : Matrix<Float> sRP(solveRPar().nonDegenerate(1));
1692 225 : Matrix<Bool> sPok(solveParOK().nonDegenerate(1));
1693 225 : Matrix<Float> sSNR(solveParSNR().nonDegenerate(1));
1694 225 : logSink() << "sPok " << sPok.shape() << LogIO::POST;
1695 :
1696 : // Map from MS antenna number to index
1697 : // transcribe fft results to sRP
1698 225 : Int ncol = drfp->param().ncolumn();
1699 225 : 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 1799 : for (Int i=0; i!=ncol; i++) { // i==iant
1706 9098 : for (Int j=0; j!=nrow; j++) { // j is parameter number
1707 7524 : Int oj = (j>=3) ? j+1 : j;
1708 7524 : sRP(IPosition(2, oj, i)) = drfp->param()(IPosition(2, j, i));
1709 7524 : 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 1574 : sPok(IPosition(2, 3, i)) = true;
1713 1574 : if (nrow > 3)
1714 934 : sPok(IPosition(2, 7, i)) = true;
1715 : }
1716 225 : size_t nCorrOrig(sdbs(0).nCorrelations());
1717 225 : size_t nCorr = (nCorrOrig> 1 ? 2 : 1); // number of p-hands
1718 225 : calculateSNR(nCorr, *drfp);
1719 225 : set<Int> belowThreshold;
1720 225 : Float threshold = minSNR();
1721 :
1722 587 : for (size_t icor=0; icor != nCorr; icor++) {
1723 362 : const set<Int>& activeAntennas = drfp->getActiveAntennasCorrelation(icor);
1724 3906 : for (Int iant=0; iant != nAnt(); iant++) {
1725 3544 : if (iant != refant() && (activeAntennas.find(iant) != activeAntennas.end())) {
1726 1906 : Float s = sSNR(4*icor + 0, iant);
1727 : // Start the log message; finished below
1728 1906 : logSink() << "Antenna " << iant << " correlation " << icor << " has (FFT) SNR of " << s;
1729 1906 : if (s < threshold) {
1730 5 : belowThreshold.insert(iant);
1731 5 : 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 5 : sPok(4*icor + 0, iant) = false;
1735 5 : sPok(4*icor + 1, iant) = false;
1736 5 : sPok(4*icor + 2, iant) = false;
1737 5 : sPok(4*icor + 3, iant) = false;
1738 : }
1739 : // Finish the log message
1740 1906 : 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 362 : drfp->removeAntennasCorrelation(icor, belowThreshold);
1746 : if (DEVDEBUG) {
1747 : drfp->printActive();
1748 : }
1749 : }
1750 225 : if (globalSolve()) {
1751 225 : 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 450 : least_squares_driver(sdbs, sRP, sPok, sSNR, refant(), drfp->getActiveAntennas(), maxits(),
1757 225 : 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 2457 : for (Int iant=0; iant != nAnt(); iant++) {
1772 5776 : for (size_t icor=0; icor != nCorr; icor++) {
1773 3544 : const set<Int>& activeAntennas = drfp->getActiveAntennasCorrelation(icor);
1774 3544 : if (activeAntennas.find(iant) == activeAntennas.end()) {
1775 1281 : continue;
1776 : }
1777 2263 : Double phi0 = sRP(4*icor + 0, iant);
1778 2263 : Double delay = sRP(4*icor + 1, iant);
1779 2263 : Double rate = sRP(4*icor + 2, iant);
1780 2263 : Double k_disp = sRP(4*icor + 3, iant);
1781 2263 : 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 2263 : Double delta1 = df0*delay/1e9;
1786 2263 : Double delta2 = reffreq*dt0*rate;
1787 2263 : Double delta3 = (2.0*M_PI)*(delta1+delta2);
1788 : Double dt;
1789 2263 : auto p = aggregateTime.find(iant);
1790 2263 : if (zeroRates() && p != aggregateTime.end()) {
1791 72 : dt = p->second - t0;
1792 : } else {
1793 2191 : 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 2263 : 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 225 : 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 405 : if (corrcomb().contains("stokes") ||
1817 180 : corrcomb().contains("parallel")) {
1818 48 : logSink() << "Correlations combined: Copying results to other correlation" << LogIO::POST;
1819 528 : for (Int iant=0; iant != nAnt(); iant++) {
1820 2400 : for (Int i=0; i !=4; i++) {
1821 1920 : sRP(4+i, iant) = sRP(i, iant);
1822 1920 : sPok(4+i, iant) = sPok(i, iant);
1823 1920 : sSNR(4+i, iant) = sSNR(i, iant);
1824 : }
1825 : }
1826 : }
1827 : if (DEVDEBUG) {
1828 : std::cerr << "sPok " << sPok << endl;
1829 : }
1830 225 : delete drfp;
1831 225 : }
1832 :
1833 : void
1834 0 : FringeJones::solveOneVB(const VisBuffer&) {
1835 0 : throw(AipsError("VisBuffer interface not supported!"));
1836 : }
1837 :
1838 :
1839 31 : void FringeJones::globalPostSolveTinker() {
1840 :
1841 : // Re-reference the phase, if requested
1842 31 : if (refantlist()(0)>-1) applyRefAnt();
1843 31 : }
1844 :
1845 31 : void FringeJones::applyRefAnt() {
1846 :
1847 : // TBD:
1848 : // 1. Synchronize refant changes on par axis
1849 : // 2. Implement minimum mean deviation algorithm
1850 :
1851 31 : if (refantlist()(0)<0)
1852 0 : throw(AipsError("No refant specified."));
1853 :
1854 31 : Int nUserRefant=refantlist().nelements();
1855 :
1856 : // Get the preferred refant names from the MS
1857 31 : String refantName(msmc().antennaName(refantlist()(0)));
1858 31 : if (nUserRefant>1) {
1859 0 : refantName+=" (";
1860 0 : for (Int i=1;i<nUserRefant;++i) {
1861 0 : refantName+=msmc().antennaName(refantlist()(i));
1862 0 : if (i<nUserRefant-1) refantName+=",";
1863 : }
1864 0 : refantName+=")";
1865 : }
1866 :
1867 : logSink() << "Applying refant: " << refantName
1868 31 : << " refantmode = " << refantmode();
1869 31 : if (refantmode()=="flex")
1870 31 : logSink() << " (hold alternate refants' phase constant) when refant flagged";
1871 31 : if (refantmode()=="strict")
1872 0 : logSink() << " (flag all antennas when refant flagged)";
1873 31 : logSink() << LogIO::POST;
1874 :
1875 : // Generate a prioritized refant choice list
1876 : // The first entry in this list is the user's primary refant,
1877 : // the second entry is the refant used on the previous interval,
1878 : // and the rest is a prioritized list of alternate refants,
1879 : // starting with the user's secondary (if provided) refants,
1880 : // followed by the rest of the array, in distance order. This
1881 : // makes the priorities correct all the time, and prevents
1882 : // a semi-stochastic alternation (by preferring the last-used
1883 : // alternate, even if nominally higher-priority refants become
1884 : // available)
1885 :
1886 :
1887 : // Extract antenna positions
1888 31 : Matrix<Double> xyz;
1889 31 : if (msName()!="<noms>") {
1890 31 : MeasurementSet ms(msName());
1891 31 : MSAntennaColumns msant(ms.antenna());
1892 31 : msant.position().getColumn(xyz);
1893 31 : }
1894 : else {
1895 : // TBD RO*
1896 0 : CTColumns ctcol(*ct_);
1897 0 : CTAntennaColumns& antcol(ctcol.antenna());
1898 0 : antcol.position().getColumn(xyz);
1899 0 : }
1900 :
1901 : // Calculate (squared) antenna distances, relative
1902 : // to last preferred antenna
1903 31 : Vector<Double> dist2(xyz.ncolumn(),0.0);
1904 124 : for (Int i=0;i<3;++i) {
1905 93 : Vector<Double> row=xyz.row(i);
1906 93 : row-=row(refantlist()(nUserRefant-1));
1907 93 : dist2+=square(row);
1908 93 : }
1909 : // Move preferred antennas to a large distance
1910 62 : for (Int i=0;i<nUserRefant;++i)
1911 31 : dist2(refantlist()(i))=DBL_MAX;
1912 :
1913 : // Generated sorted index
1914 31 : Vector<uInt> ord;
1915 31 : genSort(ord,dist2);
1916 :
1917 : // Assemble the whole choices list
1918 31 : Int nchoices=nUserRefant+1+ord.nelements();
1919 31 : Vector<Int> refantchoices(nchoices,0);
1920 62 : Vector<Int> r(refantchoices(IPosition(1,nUserRefant+1),IPosition(1,refantchoices.nelements()-1)));
1921 31 : convertArray(r,ord);
1922 :
1923 : // set first two to primary preferred refant
1924 31 : refantchoices(0)=refantchoices(1)=refantlist()(0);
1925 :
1926 : // set user's secondary refants (if any)
1927 31 : if (nUserRefant>1)
1928 0 : refantchoices(IPosition(1,2),IPosition(1,nUserRefant))=
1929 0 : refantlist()(IPosition(1,1),IPosition(1,nUserRefant-1));
1930 :
1931 : //cout << "refantchoices = " << refantchoices << endl;
1932 :
1933 :
1934 31 : if (refantmode()=="strict") {
1935 0 : nchoices=1;
1936 0 : refantchoices.resize(1,True);
1937 : }
1938 :
1939 31 : Vector<Int> nPol(nSpw(),2); // TBD:or 1, if data was single pol
1940 :
1941 31 : if (nPar()==8) {
1942 : // Verify that 2nd poln has unflagged solutions, PER SPW
1943 31 : ROCTMainColumns ctmc(*ct_);
1944 :
1945 31 : Block<String> cols(1);
1946 31 : cols[0]="SPECTRAL_WINDOW_ID";
1947 31 : CTIter ctiter(*ct_,cols);
1948 31 : Cube<Bool> fl;
1949 :
1950 125 : while (!ctiter.pastEnd()) {
1951 :
1952 94 : Int ispw=ctiter.thisSpw();
1953 94 : fl.assign(ctiter.flag());
1954 :
1955 94 : IPosition blc(3,0,0,0), trc(fl.shape());
1956 94 : trc-=1; blc(0)=nPar()/2;
1957 :
1958 : // cout << "ispw = " << ispw << " nfalse(fl(1,:,:)) = " << nfalse(fl(blc,trc)) << endl;
1959 :
1960 : // If there are no unflagged solutions in 2nd pol,
1961 : // avoid it in refant calculations
1962 94 : if (nfalse(fl(blc,trc))==0)
1963 40 : nPol(ispw)=1;
1964 :
1965 94 : ctiter.next();
1966 94 : }
1967 31 : }
1968 : // cout << "nPol = " << nPol << endl;
1969 :
1970 31 : Bool usedaltrefant(false);
1971 31 : Int currrefant(refantchoices(0)), lastrefant(-1);
1972 :
1973 31 : Block<String> cols(2);
1974 31 : cols[0]="SPECTRAL_WINDOW_ID";
1975 31 : cols[1]="TIME";
1976 31 : CTIter ctiter(*ct_,cols);
1977 :
1978 : // Arrays to hold per-timestamp solutions
1979 31 : Cube<Float> solA, solB;
1980 31 : Cube<Bool> flA, flB;
1981 31 : Vector<Int> ant1A, ant1B, ant2B;
1982 31 : Matrix<Complex> refPhsr; // the reference phasor [npol,nchan]
1983 31 : Int lastspw(-1);
1984 31 : Bool first(true);
1985 256 : while (!ctiter.pastEnd()) {
1986 225 : Int ispw=ctiter.thisSpw();
1987 225 : if (ispw!=lastspw) first=true; // spw changed, start over
1988 :
1989 : // Read in the current sol, fl, ant1:
1990 225 : solB.assign(ctiter.fparam());
1991 225 : flB.assign(ctiter.flag());
1992 225 : ant1B.assign(ctiter.antenna1());
1993 225 : ant2B.assign(ctiter.antenna2());
1994 :
1995 : // First time thru, 'previous' solution same as 'current'
1996 225 : if (first) {
1997 94 : solA.reference(solB);
1998 94 : flA.reference(flB);
1999 94 : ant1A.reference(ant1B);
2000 : }
2001 225 : IPosition shB(solB.shape());
2002 225 : IPosition shA(solA.shape());
2003 :
2004 : // Find a good refant at this time
2005 : // A good refant is one that is unflagged in all polarizations
2006 : // in the current(B) and previous(A) intervals (so they can be connected)
2007 225 : Int irefA(0),irefB(0); // index on 3rd axis of solution arrays
2008 225 : Int ichoice(0); // index in refantchoicelist
2009 225 : Bool found(false);
2010 225 : IPosition blcA(3,0,0,0),trcA(shA),blcB(3,0,0,0),trcB(shB);
2011 225 : trcA-=1; trcA(0)=trcA(2)=0;
2012 225 : trcB-=1; trcB(0)=trcB(2)=0;
2013 225 : ichoice=0;
2014 450 : while (!found && ichoice<nchoices) {
2015 : // Find index of current refant choice
2016 225 : irefA=irefB=0;
2017 229 : while (ant1A(irefA)!=refantchoices(ichoice) && irefA<shA(2)) ++irefA;
2018 229 : while (ant1B(irefB)!=refantchoices(ichoice) && irefB<shB(2)) ++irefB;
2019 :
2020 225 : if (irefA<shA(2) && irefB<shB(2)) {
2021 :
2022 : // cout << " Trial irefA,irefB: " << irefA << "," << irefB
2023 : // << " Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;
2024 :
2025 225 : blcA(2)=trcA(2)=irefA;
2026 225 : blcB(2)=trcB(2)=irefB;
2027 225 : found=true; // maybe
2028 635 : for (Int ipol=0;ipol<nPol(ispw);++ipol) {
2029 410 : blcA(0)=blcB(0)=ipol*nPar()/2;
2030 410 : trcA(0)=trcB(0)=(ipol+1)*nPar()/2-1;
2031 410 : found &= (nfalse(flA(blcA,trcA))>0); // previous interval
2032 410 : found &= (nfalse(flB(blcB,trcB))>0); // current interval
2033 : }
2034 : }
2035 : else
2036 : // irefA or irefB out-of-range
2037 0 : found=false; // Just to be sure
2038 :
2039 225 : if (!found) ++ichoice; // try next choice next round
2040 :
2041 : }
2042 :
2043 225 : if (found) {
2044 : // at this point, irefA/irefB point to a good refant
2045 :
2046 : // Keep track
2047 225 : usedaltrefant|=(ichoice>0);
2048 225 : currrefant=refantchoices(ichoice);
2049 225 : refantchoices(1)=currrefant; // 2nd priorty next time
2050 :
2051 : // cout << " currrefant = " << currrefant << " (" << ichoice << ")" << endl;
2052 :
2053 : // cout << " Final irefA,irefB: " << irefA << "," << irefB
2054 : // << " Ants=" << ant1A(irefA) << "," << ant1B(irefB) << endl;
2055 :
2056 :
2057 : // Only report if using an alternate refant
2058 225 : if (currrefant!=lastrefant && ichoice>0) {
2059 : logSink()
2060 : << "At "
2061 0 : << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2062 : << " ("
2063 : << "Spw=" << ctiter.thisSpw()
2064 : << ", Fld=" << ctiter.thisField()
2065 : << ")"
2066 0 : << ", using refant " << msmc().antennaName(currrefant)
2067 : << " (id=" << currrefant
2068 : << ")" << " (alternate)"
2069 0 : << LogIO::POST;
2070 : }
2071 :
2072 : // Form reference phasor [nPar,nChan]
2073 225 : Matrix<Float> rA,rB;
2074 225 : Matrix<Bool> rflA,rflB;
2075 225 : rB.assign(solB.xyPlane(irefB));
2076 225 : rflB.assign(flB.xyPlane(irefB));
2077 :
2078 225 : if (!first) {
2079 : // Get and condition previous phasor for the current refant
2080 131 : rA.assign(solA.xyPlane(irefA));
2081 131 : rflA.assign(flA.xyPlane(irefA));
2082 131 : rB-=rA;
2083 :
2084 : // Accumulate flags
2085 131 : rflB&=rflA;
2086 : }
2087 :
2088 : // cout << " rB = " << rB << endl;
2089 : // cout << boolalpha << " rflB = " << rflB << endl;
2090 : // TBD: fillChanGaps?
2091 :
2092 : // Now apply reference phasor to all antennas
2093 225 : Matrix<Float> thissol;
2094 2457 : for (Int iant=0;iant<shB(2);++iant) {
2095 2232 : thissol.reference(solB.xyPlane(iant));
2096 2232 : thissol-=rB;
2097 : }
2098 :
2099 : // Set refant, so we can put it back
2100 225 : ant2B=currrefant;
2101 :
2102 : // put back referenced solutions
2103 225 : ctiter.setfparam(solB);
2104 225 : ctiter.setantenna2(ant2B);
2105 :
2106 : // Next time thru, solB is previous
2107 225 : solA.reference(solB);
2108 225 : flA.reference(flB);
2109 225 : ant1A.reference(ant1B);
2110 225 : solB.resize(); // (break references)
2111 225 : flB.resize();
2112 225 : ant1B.resize();
2113 :
2114 225 : lastrefant=currrefant;
2115 225 : first=false; // avoid first-pass stuff from now on
2116 :
2117 225 : } // found
2118 : else {
2119 : logSink()
2120 : << "At "
2121 0 : << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2122 : << " ("
2123 : << "Spw=" << ctiter.thisSpw()
2124 : << ", Fld=" << ctiter.thisField()
2125 : << ")"
2126 : << ", refant (id=" << currrefant
2127 : << ") was flagged; flagging all antennas strictly."
2128 0 : << LogIO::POST;
2129 : // Flag all solutions in this interval
2130 0 : flB.set(True);
2131 0 : ctiter.setflag(flB);
2132 0 : if (ichoice == nchoices){
2133 : logSink() << LogIO::WARN
2134 : << "From time: "
2135 0 : << MVTime(ctiter.thisTime()/C::day).string(MVTime::YMD,7)
2136 : << "in Spw: "
2137 : << ctiter.thisSpw()
2138 : << " refant list exhausted, flagging all solutions"
2139 0 : << LogIO::POST;
2140 : }
2141 : }
2142 :
2143 : // advance to the next interval
2144 225 : lastspw=ispw;
2145 225 : ctiter.next();
2146 225 : }
2147 :
2148 31 : if (usedaltrefant)
2149 : logSink() << LogIO::NORMAL
2150 : << " NB: An alternate refant was used at least once to maintain" << endl
2151 : << " phase continuity where the user's preferred refant drops out." << endl
2152 : << " Alternate refants are held constant in phase (_not_ zeroed)" << endl
2153 : << " during these periods, and the preferred refant may return at" << endl
2154 : << " a non-zero phase. This is generally harmless."
2155 0 : << LogIO::POST;
2156 :
2157 62 : return;
2158 :
2159 31 : }
2160 :
2161 : // CT Smooth for fringefit table
2162 : // Some duplication from CTglobals Smooth
2163 :
2164 0 : void FringeJones::smooth(Vector<Int>& fields,
2165 : const String& smtype,
2166 : const Double& smtime,
2167 : const bool ratesmooth) {
2168 0 : NewCalTable ct = *ct_;
2169 :
2170 : // half-width
2171 0 : Double thw(smtime/2.0);
2172 :
2173 : // Workspace
2174 0 : Vector<Double> times;
2175 0 : Vector<Float> p,newp,pRate, newRate;
2176 0 : Vector<Float> d;
2177 0 : Vector<Bool> pOK, newpOK;
2178 : // Unwrapped
2179 0 : vector<vector<float>> temp;
2180 :
2181 0 : Cube<Float> fpar;
2182 0 : Cube<Bool> fparok,newfparok;
2183 :
2184 0 : Vector<Bool> mask;
2185 :
2186 : // Set up spw col for ref freq
2187 0 : MSSpectralWindow msSpw(ct.spectralWindow());
2188 0 : MSSpWindowColumns msCol(msSpw);
2189 :
2190 0 : IPosition blc(3,0,0,0), fblc(3,0,0,0);
2191 0 : IPosition trc(3,0,0,0), ftrc(3,0,0,0);
2192 0 : IPosition vec(1,0);
2193 :
2194 0 : Block<String> cols(4);
2195 0 : cols[0]="SPECTRAL_WINDOW_ID";
2196 0 : cols[1]="FIELD_ID";
2197 0 : cols[2]="ANTENNA1";
2198 0 : cols[3]="ANTENNA2";
2199 0 : CTIter ctiter(ct,cols);
2200 0 : int counter = 0;
2201 :
2202 0 : while (!ctiter.pastEnd()) {
2203 :
2204 : //MSSpectralWindow msSpw(ct.spectralWindow());
2205 : //MSSpWindowColumns msCol(msSpw);
2206 :
2207 0 : Int nSlot=ctiter.nrow();
2208 0 : Int ifld=ctiter.thisField();
2209 0 : Int ispw=ctiter.thisSpw();
2210 :
2211 : // Only if more than one slot in this spw _AND_
2212 : // field is among those requested (if any)
2213 0 : if (nSlot>1 &&
2214 0 : (fields.nelements()<1 || anyEQ(fields,ifld))) {
2215 0 : vec(0)=nSlot;
2216 0 : trc(2)=ftrc(2)=nSlot-1;
2217 :
2218 0 : times.assign(ctiter.time());
2219 :
2220 0 : fpar.assign(ctiter.fparam());
2221 0 : fparok.assign(!ctiter.flag());
2222 0 : newfparok.assign(fparok);
2223 0 : IPosition fsh(fpar.shape());
2224 :
2225 0 : blc(1)=trc(1)=fblc(1)=ftrc(1)=0;
2226 :
2227 : // get chan Freqs
2228 0 : Vector<Double> freqs;
2229 0 : msCol.chanFreq().get(ispw,freqs,True);
2230 0 : Double refFreq = freqs(0);
2231 :
2232 : // For each param (pol)
2233 0 : counter = 0;
2234 0 : temp.clear();
2235 0 : vector<vector<float>> unwrap(2);
2236 : //int polId = 0;
2237 0 : bool polReset = false;
2238 :
2239 : // Smoothing beforehand for rates
2240 :
2241 :
2242 :
2243 : // Need a seperate iter to construct unwrapped phase estimates
2244 : // Iterate over polId rather than ipar
2245 0 : for (Int polId=0; polId<2;++polId) {
2246 0 : temp.clear();
2247 0 : counter = 0;
2248 :
2249 0 : blc(0)=trc(0)=polId * 4;
2250 0 : fblc(0)=ftrc(0)=(polId * 4) + 2;
2251 :
2252 : // Reference slices of par twice. Once for pol and again for delay rates
2253 0 : p.reference(fpar(blc,trc).reform(vec));
2254 0 : newp.assign(p);
2255 0 : pRate.reference(fpar(fblc,ftrc).reform(vec));
2256 0 : newRate.assign(pRate);
2257 0 : pOK.reference(fparok(fblc,ftrc).reform(vec));
2258 0 : newpOK.reference(newfparok(fblc,ftrc).reform(vec));
2259 :
2260 0 : Vector<Bool> mask;
2261 0 : int cycles = 0;
2262 :
2263 0 : for (Int i=0;i<nSlot;++i) {
2264 : // holder for phase, delay, and time
2265 0 : vector<float> holder {0.0, 0.0, 0.0};
2266 : // mask for rate smoothing
2267 : // Make mask
2268 0 : mask = pOK;
2269 0 : mask = (mask && ( (times > (times(i)-thw)) &&
2270 0 : (times <= (times(i)+thw)) ) );
2271 :
2272 :
2273 : // Save the phase rate and time to use for the estimates
2274 :
2275 0 : holder[0] = newp(i);
2276 0 : holder[1] = pRate(i);
2277 0 : holder[2] = times(i);
2278 :
2279 : // Smooth the rates
2280 0 : if (ntrue(mask)>0) {
2281 0 : if (smtype=="mean") {
2282 0 : pRate(i) = mean(newRate(mask));
2283 : }
2284 0 : else if (smtype=="median") {
2285 0 : pRate(i) = median(newRate(mask), false);
2286 : }
2287 : }
2288 :
2289 : // array of phases delays and times to be used in the cycle estimations
2290 0 : temp.push_back(holder);
2291 :
2292 : // Estimate the number of Phase cycles
2293 0 : if (counter == 0) {
2294 : // If we are at 0 we cant interpolate backwards. Just instert as starting value
2295 0 : unwrap[polId].push_back(temp[counter][0]);
2296 : }
2297 : else {
2298 : // Get the time difference between two points
2299 0 : float timeStep = temp[counter][2] - temp[counter-1][2];
2300 : // Get Forwards and backwards predictions (in cycles)
2301 0 : float predictFWDiff = ((temp[counter-1][0]/(2*M_PI)) + (temp[counter-1][1] * refFreq * timeStep * 2) * int(ratesmooth)) - (temp[counter][0]/(2*M_PI));
2302 0 : float predictBWDiff = ((temp[counter][0]/(2*M_PI)) - (temp[counter][1] * refFreq * timeStep * 2) * int(ratesmooth)) - (temp[counter-1][0]/(2*M_PI));
2303 : // Take the average prediction of cycles
2304 0 : float cycleDiff = ((predictFWDiff-predictBWDiff)/2);
2305 : // Adjust total cycle estimate
2306 0 : if (cycleDiff > 0.5) {
2307 0 : cycles += (int)(floor(cycleDiff + 0.5));
2308 : }
2309 0 : if (cycleDiff < 0.5) {
2310 0 : cycles -= (int)(floor(abs(cycleDiff) + 0.5));
2311 : }
2312 :
2313 : // Add unwrapped phase values to array
2314 0 : unwrap[polId].push_back(temp[counter][0] + 2 * M_PI * cycles);
2315 : }
2316 :
2317 : // increment counter
2318 0 : counter++;
2319 0 : }
2320 0 : }
2321 :
2322 : // Convert unwrap to casa Vector so we can use the same mean and masking functions
2323 0 : Vector<Float> unwrapPhasesPol1(unwrap[0]);
2324 0 : Vector<Float> unwrapPhasesPol2(unwrap[1]);
2325 :
2326 : // Regular ipar interation
2327 0 : for (Int ipar=0;ipar<fsh(0);++ipar) {
2328 0 : blc(0)=trc(0)=ipar;
2329 0 : fblc(0)=ftrc(0)=ipar;
2330 :
2331 : // Reference slices of par/parOK
2332 0 : p.reference(fpar(blc,trc).reform(vec));
2333 0 : newp.assign(p);
2334 0 : pOK.reference(fparok(fblc,ftrc).reform(vec));
2335 0 : newpOK.reference(newfparok(fblc,ftrc).reform(vec));
2336 :
2337 0 : Vector<Bool> mask;
2338 :
2339 : //cout << "IPAR: " << ipar << "\n"
2340 : //<< "VAL: " << newp << "\n" << endl;
2341 :
2342 0 : for (Int i=0;i<nSlot;++i) {
2343 : // Make mask
2344 0 : mask = pOK;
2345 0 : mask = (mask && ( (times > (times(i)-thw)) &&
2346 0 : (times <= (times(i)+thw)) ) );
2347 :
2348 0 : if (ntrue(mask)>0) {
2349 0 : if (smtype=="mean") {
2350 :
2351 : // If phases use our unwrapped vector
2352 0 : if (ipar==0) {newp(i)=mean(unwrapPhasesPol1(mask));};
2353 0 : if (ipar==4) {newp(i)=mean(unwrapPhasesPol2(mask));};
2354 :
2355 0 : if (ipar == 0 || ipar == 4){
2356 0 : while (newp(i) < -M_PI) {
2357 0 : newp(i) += 2*M_PI;
2358 : }
2359 0 : while (newp(i) > M_PI) {
2360 0 : newp(i) -= 2*M_PI;
2361 : }
2362 : }
2363 : else{
2364 0 : newp(i)=mean(p(mask));
2365 : }
2366 :
2367 : }
2368 0 : else if (smtype=="median") {
2369 0 : if (ipar==0) {newp(i)=median(unwrapPhasesPol1(mask),false);};
2370 0 : if (ipar==4) {newp(i)=median(unwrapPhasesPol2(mask),false);};
2371 :
2372 : if (ipar == 0 || 4) {
2373 0 : while (newp(i) < -M_PI) {
2374 0 : newp(i) += 2*M_PI;
2375 : }
2376 0 : while (newp(i) > M_PI) {
2377 0 : newp(i) -= 2*M_PI;
2378 : }
2379 : }
2380 : else {
2381 : newp(i)= median(p(mask),false);
2382 : }
2383 : }
2384 0 : newpOK(i)=true;
2385 : }
2386 : else
2387 0 : newpOK(i)=false;
2388 :
2389 : } // i
2390 :
2391 : // keep new ok info
2392 0 : p=newp;
2393 0 : } // ipar
2394 :
2395 : // Put info back
2396 0 : ctiter.setfparam(fpar);
2397 :
2398 0 : ctiter.setflag(!newfparok);
2399 :
2400 0 : } // nSlot>1
2401 :
2402 0 : ctiter.next();
2403 : } // ispw
2404 0 : }
2405 :
2406 :
2407 :
2408 : } //# NAMESPACE CASA - END
2409 :
2410 :
2411 : /*
2412 : Example of use in at Python level:
2413 :
2414 : fringefit(vis="n14c2.ms", caltable="fail.fj", field="",spw="1",intent="",
2415 : selectdata=True, timerange="", antenna="", scan="5", observation="",
2416 : msselect="", solint="inf", refant="EF", minsnr=3.0, append=False,
2417 : gaintable=['n14c2.gcal'], parang=False)
2418 : */
2419 :
|
