Line data Source code
1 : //# VBContinuumSubtractor.cc: Subtract continuum from spectral line data
2 : //# Copyright (C) 2004
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 : //# $Id$
27 : //#
28 : //#include <casa/Arrays/ArrayLogical.h>
29 : //#include <casa/Arrays/ArrayMath.h>
30 : //#include <casa/Arrays/ArrayUtil.h>
31 : #include <casacore/casa/Arrays/Cube.h>
32 : //#include <casa/Arrays/MaskedArray.h>
33 : //#include <casa/Arrays/MaskArrMath.h>
34 : //#include <casa/Containers/Record.h>
35 : //#include <casa/Containers/RecordFieldId.h>
36 : //#include <casa/Exceptions/Error.h>
37 : #include <casacore/casa/Logging/LogIO.h>
38 : //#include <casa/Quanta/MVTime.h>
39 : //#include <casa/Quanta/QuantumHolder.h>
40 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
41 : #include <casacore/ms/MeasurementSets/MSColumns.h>
42 : #include <msvis/MSVis/VBContinuumSubtractor.h>
43 : #include <msvis/MSVis/VisBuffGroupAcc.h>
44 : #include <msvis/MSVis/VisBuffer.h>
45 : #include <casacore/scimath/Fitting/LinearFitSVD.h>
46 : #include <casacore/scimath/Functionals/Polynomial.h>
47 :
48 : #include <msvis/MSVis/VisBuffer2.h>
49 :
50 : //#include <algorithm>
51 :
52 : using namespace casacore;
53 : using namespace casa::vi;
54 :
55 : namespace casa { //# NAMESPACE CASA - BEGIN
56 :
57 0 : VBContinuumSubtractor::VBContinuumSubtractor():
58 0 : fitorder_p(-1),
59 0 : lofreq_p(-1.0),
60 0 : hifreq_p(-1.0),
61 0 : midfreq_p(-1.0),
62 0 : freqscale_p(1.0),
63 0 : maxAnt_p(-1),
64 0 : nHashes_p(0),
65 0 : ncorr_p(0),
66 0 : totnumchan_p(0),
67 0 : tvi_debug(False)
68 : {
69 0 : }
70 :
71 0 : VBContinuumSubtractor::VBContinuumSubtractor(const Double lofreq,
72 0 : const Double hifreq):
73 0 : fitorder_p(-1),
74 0 : lofreq_p(lofreq),
75 0 : hifreq_p(hifreq),
76 0 : midfreq_p(-1.0),
77 0 : freqscale_p(1.0),
78 0 : maxAnt_p(-1),
79 0 : nHashes_p(0),
80 0 : ncorr_p(0),
81 0 : totnumchan_p(0),
82 0 : tvi_debug(False)
83 : {
84 0 : midfreq_p = 0.5 * (lofreq + hifreq);
85 0 : freqscale_p = calcFreqScale();
86 0 : }
87 :
88 0 : void VBContinuumSubtractor::resize(Cube<Complex>& coeffs,
89 : Cube<Bool>& coeffsOK) const
90 : {
91 0 : if(maxAnt_p < 0 || fitorder_p < 0 || ncorr_p < 1)
92 0 : throw(AipsError("The fit order, # of corrs, and/or max antenna # must be set."));
93 :
94 : // An nth order polynomial has n + 1 coefficients.
95 0 : coeffs.resize(ncorr_p, fitorder_p + 1, nHashes_p);
96 : // Calibrater wants coeffsOK to be a Cube, even though a Matrix would do for
97 : // VBContinuumSubtractor. Unfortunately problems arise (worse, quietly) in
98 : // SolvableVisCal::keep() and store() if one tries to get away with
99 : //coeffsOK.resize(ncorr_p, 1, nHashes_p);
100 0 : coeffsOK.resize(ncorr_p, fitorder_p + 1, nHashes_p);
101 0 : }
102 :
103 0 : void VBContinuumSubtractor::init(const IPosition& shp, const uInt maxAnt,
104 : const uInt totnumchan,
105 : const Double lof, const Double hif)
106 : {
107 0 : ncorr_p = shp[0];
108 0 : fitorder_p = shp[1] - 1;
109 :
110 : //// Going from the number of baselines to the number of antennas is a little
111 : //// backwards, but so is this function.
112 : // uInt nAnt = round((-1 + sqrt(1 + 8 * shp[2])) / 2);
113 0 : setNAnt(maxAnt + 1);
114 :
115 0 : totnumchan_p = totnumchan;
116 0 : setScalingFreqs(lof, hif);
117 0 : }
118 :
119 0 : Bool VBContinuumSubtractor::initFromVBGA(VisBuffGroupAcc& vbga)
120 : {
121 0 : Bool retval = true;
122 :
123 0 : if(vbga.nBuf() > 0){
124 0 : ncorr_p = vbga(0).nCorr();
125 0 : setNAnt(vbga.nAnt());
126 :
127 : // Count the total number of channels, and get the minimum and maximum
128 : // frequencies for scaling.
129 0 : totnumchan_p = 0;
130 0 : hifreq_p = -1.0;
131 0 : lofreq_p = DBL_MAX;
132 0 : for(Int ibuf = 0; ibuf < vbga.nBuf(); ++ibuf){
133 0 : VisBuffer& vb(vbga(ibuf));
134 :
135 0 : totnumchan_p += vb.nChannel();
136 0 : getMinMaxFreq(vb, lofreq_p, hifreq_p, false);
137 : }
138 0 : midfreq_p = 0.5 * (lofreq_p + hifreq_p);
139 0 : freqscale_p = calcFreqScale();
140 : }
141 : else
142 0 : retval = false;
143 0 : return retval;
144 : }
145 :
146 0 : VBContinuumSubtractor::~VBContinuumSubtractor()
147 0 : {}
148 :
149 0 : void VBContinuumSubtractor::fit(VisBuffGroupAcc& vbga, const Int fitorder,
150 : MS::PredefinedColumns whichcol,
151 : Cube<Complex>& coeffs,
152 : Cube<Bool>& coeffsOK, const Bool doInit,
153 : const Bool doResize,
154 : const Bool squawk)
155 : {
156 0 : LogIO os(LogOrigin("VBContinuumSubtractor", "fit()", WHERE));
157 :
158 : // jagonzal: Debug code
159 : /*
160 : uInt row_debug = 0;
161 : uInt corr_debug = 0;
162 : */
163 :
164 0 : fitorder_p = fitorder;
165 :
166 0 : if(!(whichcol == MS::DATA || whichcol == MS::MODEL_DATA ||
167 : whichcol == MS::CORRECTED_DATA)){
168 0 : if(squawk)
169 : os << LogIO::SEVERE
170 : << MS::columnName(whichcol) << " is not supported.\n"
171 : << MS::columnName(MS::DATA) << " will be used instead."
172 0 : << LogIO::POST;
173 0 : whichcol = MS::DATA;
174 : }
175 :
176 0 : if(doInit)
177 0 : initFromVBGA(vbga);
178 :
179 0 : if(maxAnt_p < 0 || fitorder_p < 0 || ncorr_p < 1 || totnumchan_p < 1
180 0 : || lofreq_p < 0.0 || hifreq_p < 0.0)
181 0 : throw(AipsError("The continuum fitter must first be initialized."));
182 :
183 0 : if(doResize)
184 0 : resize(coeffs, coeffsOK);
185 :
186 0 : if(!checkSize(coeffs, coeffsOK))
187 0 : throw(AipsError("Shape mismatch in the coefficient storage cubes."));
188 :
189 : // Make the estimate
190 0 : LinearFitSVD<Float> fitter;
191 0 : fitter.asWeight(true); // Makes the "sigma" arg = w = 1/sig**2
192 :
193 0 : coeffsOK.set(false);
194 :
195 : // Translate vbga to arrays for use by LinearFitSVD.
196 :
197 : // The fitorder will actually be clamped on a baseline-by-baseline basis
198 : // because of flagging, but a summary note is in order here.
199 0 : if(static_cast<Int>(totnumchan_p) < fitorder_p)
200 : os << LogIO::WARN
201 : << "fitorder = " << fitorder_p
202 : << ", but only " << totnumchan_p << " channels were selected.\n"
203 : << "The polynomial order will be lowered accordingly."
204 0 : << LogIO::POST;
205 : // Scale frequencies to [-1, 1].
206 0 : midfreq_p = 0.5 * (lofreq_p + hifreq_p);
207 0 : freqscale_p = calcFreqScale();
208 0 : Vector<Float> freqs(totnumchan_p);
209 0 : uInt totchan = 0;
210 0 : for(Int ibuf = 0; ibuf < vbga.nBuf(); ++ibuf){
211 0 : VisBuffer& vb(vbga(ibuf));
212 0 : Vector<Double> freq(vb.frequency());
213 0 : uInt nchan = vb.nChannel();
214 :
215 0 : for(uInt c = 0; c < nchan; ++c){
216 0 : freqs[totchan] = freqscale_p * (freq[c] - midfreq_p);
217 0 : ++totchan;
218 : }
219 0 : }
220 :
221 0 : Vector<Float> wt(totnumchan_p);
222 0 : Vector<Float> unflaggedfreqs(totnumchan_p);
223 0 : Vector<Complex> vizzes(totnumchan_p);
224 0 : Vector<Float> floatvs(totnumchan_p);
225 0 : Vector<Float> realsolution(fitorder_p + 1);
226 0 : Vector<Float> imagsolution(fitorder_p + 1);
227 :
228 0 : for(uInt corrind = 0; corrind < ncorr_p; ++corrind){
229 0 : for(uInt blind = 0; blind < nHashes_p; ++blind){
230 0 : uInt totchan = 0;
231 0 : uInt totunflaggedchan = 0;
232 :
233 : // Fill wt, unflaggedfreqs, and vizzes with the baseline's values for
234 : // all channels being used in the fit.
235 0 : wt.resize(totnumchan_p);
236 0 : vizzes.resize(totnumchan_p);
237 0 : unflaggedfreqs.resize(totnumchan_p);
238 :
239 0 : for(Int ibuf = 0; ibuf < vbga.nBuf(); ++ibuf){
240 0 : VisBuffer& vb(vbga(ibuf));
241 0 : uInt nchan = vb.nChannel();
242 : //Int vbrow = vbga.outToInRow(ibuf, false)[blind];
243 :
244 0 : if(!vb.flagRow()[blind]){
245 0 : Cube<Complex>& viscube(vb.dataCube(whichcol));
246 : Float w;
247 :
248 : // 2/24/2011: VisBuffer doesn't (yet) have sigmaSpectrum, and I have
249 : // never seen it in an MS anyway. Settle for 1/sqrt(weightSpectrum)
250 : // if it is available or sigmaMat otherwise.
251 : //const Bool haveWS = vb.existsWeightSpectrum();
252 : // 5/13/2011: Sigh. VisBuffAccumulator doesn't even handle
253 : // WeightSpectrum, let alone sigma.
254 : //const Bool haveWS = false;
255 :
256 : //if(!haveWS) // w is needed either way, in case ws == 0.0.
257 0 : w = vb.weightMat()(corrind, blind);
258 :
259 : // w needs a sanity check, because a VisBuffer from vbga is not
260 : // necessarily still attached to the MS and sigmaMat() is not one
261 : // of the accumulated quantities. This caused problems for
262 : // the last integration in CAS-3135. checkVisIter() didn't do the
263 : // trick in that case. Fortunately w isn't all that important; if
264 : // all the channels have the same weight the only consequence of
265 : // setting w to 1 is that the estimated errors (which we don't yet
266 : // use) will be wrong.
267 : //
268 : // 5e-45 ended up getting squared in the fitter and producing a NaN.
269 0 : if(isnan(w) || w < 1.0e-20 || w > 1.0e20)
270 0 : w = 1.0; // Emit a warning?
271 :
272 : Bool chanFlag;
273 0 : for(uInt c = 0; c < nchan; ++c){
274 : // AAARRGGGHHH!! With Calibrater you have to use vb.flag(), not
275 : // flagCube(), to get the channel selection!
276 : //if(!vb.flagCube()(corrind, c, vbrow)){0
277 0 : if (!tvi_debug)
278 : {
279 0 : chanFlag = vb.flag()(c, blind);
280 : }
281 : else
282 : {
283 : // jagonzal: Use flagCube, meaning that each correlation is fit
284 : // independently even if the other correlations are fully flagged
285 0 : chanFlag = vb.flagCube()(corrind, c, blind);
286 : // jagonzal: Apparently line-free chan mask is contained in vb.flag()
287 0 : chanFlag = vb.flagCube()(corrind, c, blind) | vb.flag()(c, blind);
288 : }
289 : //if(!vb.flag()(c, blind)){
290 : // if(!vb.flagCube()(corrind, c, blind)){
291 0 : if (!chanFlag) {
292 0 : unflaggedfreqs[totunflaggedchan] = freqs[totchan];
293 : // if(haveWS){
294 : // Double ws = vb.weightSpectrum()(corrind, c, vbrow);
295 :
296 : // wt[totunflaggedchan] = ws;
297 : // }
298 : // else
299 0 : wt[totunflaggedchan] = w / nchan;
300 0 : vizzes[totunflaggedchan] = viscube(corrind, c, blind);
301 0 : ++totunflaggedchan;
302 : }
303 0 : ++totchan;
304 : }
305 : }
306 : else
307 0 : totchan += nchan;
308 : }
309 :
310 : // jagonzal: Debug code
311 : /*
312 : if (tvi_debug)
313 : {
314 : VisBuffer& vb(vbga(0));
315 : uInt rowId = vb.rowIds()[blind];
316 : if (rowId == row_debug and corrind == corr_debug and totunflaggedchan <= 0)
317 : {
318 : os << LogIO::WARN << "All channels flagged" << LogIO::POST;
319 : }
320 : }
321 : */
322 :
323 0 : if(totunflaggedchan > 0){ // OK, try a fit.
324 : // Truncate the Vectors.
325 0 : wt.resize(totunflaggedchan, true);
326 : //vizzes.resize(totunflaggedchan, true);
327 0 : floatvs.resize(totunflaggedchan);
328 0 : unflaggedfreqs.resize(totunflaggedchan, true);
329 :
330 : // perform least-squares fit of a polynomial.
331 : // Don't try to solve for more coefficients than valid channels.
332 0 : Int locFitOrd = min(fitorder_p, static_cast<Int>(totunflaggedchan) - 1);
333 :
334 : // if(locFitOrd < 1)
335 : // os << LogIO::DEBUG1
336 : // << "locFitOrd = " << locFitOrd
337 : // << LogIO::POST;
338 :
339 0 : Polynomial<AutoDiff<Float> > pnom(locFitOrd);
340 :
341 : // The way LinearFit is templated, "y" can be Complex, but at the cost
342 : // of "x" being Complex as well, and worse, wt too. It is better to
343 : // separately fit the reals and imags.
344 : // Do reals.
345 0 : for(Int ordind = 0; ordind <= locFitOrd; ++ordind) // Note <=.
346 0 : pnom.setCoefficient(ordind, 1.0);
347 :
348 0 : for(uInt c = 0; c < totunflaggedchan; ++c)
349 0 : floatvs[c] = vizzes[c].real();
350 :
351 0 : fitter.setFunction(pnom);
352 0 : realsolution(Slice(0,locFitOrd+1,1)) = fitter.fit(unflaggedfreqs, floatvs, wt);
353 :
354 : // jagonzal: Debug code
355 : /*
356 : if (tvi_debug)
357 : {
358 : VisBuffer& vb(vbga(0));
359 : uInt rowId = vb.rowIds()[blind];
360 : if (rowId == row_debug and corrind == corr_debug)
361 : {
362 : os << "locFitOrd=" << locFitOrd << LogIO::POST;
363 : os << "flag=" << vb.flag().column(blind) << LogIO::POST;
364 : os << "flagCube=" << vb.flagCube().xyPlane(blind).row(corrind) << LogIO::POST;
365 : os << "unflaggedfreqs=" << unflaggedfreqs << LogIO::POST;
366 : os << "wt=" << wt<< LogIO::POST;
367 : os << "real floatvs=" << floatvs << LogIO::POST;
368 : os << "real realsolution=" << realsolution << LogIO::POST;
369 : }
370 : }
371 : */
372 :
373 : // if(isnan(realsolution[0])){
374 : // os << LogIO::DEBUG1 << "NaN found." << LogIO::POST;
375 : // for(uInt c = 0; c < totunflaggedchan; ++c){
376 : // if(isnan(unflaggedfreqs[c]))
377 : // os << LogIO::DEBUG1
378 : // << "unflaggedfreqs[" << c << "] is a NaN."
379 : // << LogIO::POST;
380 : // if(isnan(floatvs[c]))
381 : // os << LogIO::DEBUG1
382 : // << "floatvs[" << c << "] is a NaN."
383 : // << LogIO::POST;
384 : // if(isnan(wt[c]))
385 : // os << LogIO::DEBUG1
386 : // << "wt[" << c << "] is a NaN."
387 : // << LogIO::POST;
388 : // else if(wt[c] <= 0.0)
389 : // os << LogIO::DEBUG1
390 : // << "wt[" << c << "] = " << wt[c]
391 : // << LogIO::POST;
392 : // }
393 : // }
394 :
395 : // Do imags.
396 0 : for(Int ordind = 0; ordind <= locFitOrd; ++ordind) // Note <=.
397 0 : pnom.setCoefficient(ordind, 1.0);
398 :
399 0 : for(uInt c = 0; c < totunflaggedchan; ++c)
400 0 : floatvs[c] = vizzes[c].imag();
401 :
402 0 : fitter.setFunction(pnom);
403 0 : imagsolution(Slice(0,locFitOrd+1,1)) = fitter.fit(unflaggedfreqs, floatvs, wt);
404 :
405 : // jagonzal: Debug code
406 : /*
407 : if (tvi_debug)
408 : {
409 : VisBuffer& vb(vbga(0));
410 : uInt rowId = vb.rowIds()[blind];
411 : if (rowId == row_debug and corrind == corr_debug)
412 : {
413 : os << "imag floatvs=" << floatvs << LogIO::POST;
414 : os << "imag solution=" << imagsolution << LogIO::POST;
415 : }
416 : }
417 : */
418 :
419 :
420 0 : for(Int ordind = 0; ordind <= locFitOrd; ++ordind){ // Note <=.
421 0 : coeffs(corrind, ordind, blind) = Complex(realsolution[ordind],
422 0 : imagsolution[ordind]);
423 0 : coeffsOK(corrind, ordind, blind) = true;
424 : }
425 :
426 : // Pad remaining orders (if any) with 0.0. Note <=.
427 0 : for(Int ordind = locFitOrd + 1; ordind <= fitorder_p; ++ordind){
428 0 : coeffs(corrind, ordind, blind) = 0.0;
429 :
430 : // Since coeffs(corrind, ordind, blind) == 0, it isn't necessary to
431 : // pay attention to coeffsOK(corrind, ordind, blind) (especially?) if
432 : // ordind > 0. But Calibrater's SolvableVisCal::keep() and store()
433 : // quietly go awry if you try coeffsOK.resize(ncorr_p, 1, nHashes_p);
434 0 : coeffsOK(corrind, ordind, blind) = false;
435 : }
436 :
437 : // TODO: store uncertainties
438 0 : }
439 : }
440 : }
441 0 : }
442 :
443 0 : void VBContinuumSubtractor::getMinMaxFreq(VisBuffer& vb,
444 : Double& minfreq,
445 : Double& maxfreq,
446 : const Bool initialize)
447 : {
448 0 : const Vector<Double>& freq(vb.frequency());
449 0 : Int hichan = vb.nChannel() - 1;
450 0 : Int lochan = 0;
451 :
452 0 : if(initialize){
453 0 : maxfreq = -1.0;
454 0 : minfreq = DBL_MAX;
455 : }
456 :
457 0 : if(freq[hichan] < freq[lochan]){
458 0 : lochan = hichan;
459 0 : hichan = 0;
460 : }
461 0 : if(freq[hichan] > maxfreq)
462 0 : maxfreq = freq[hichan];
463 0 : if(freq[lochan] < minfreq)
464 0 : minfreq = freq[lochan];
465 0 : }
466 :
467 0 : Bool VBContinuumSubtractor::areFreqsInBounds(VisBuffer& vb,
468 : const Bool squawk) const
469 : {
470 : Double maxfreq, minfreq;
471 :
472 0 : getMinMaxFreq(vb, minfreq, maxfreq);
473 0 : Bool result = minfreq >= lofreq_p && maxfreq <= hifreq_p;
474 :
475 0 : if(squawk && !result){
476 : // The truth was considered too alarming (CAS-1968).
477 : // LogIO os(LogOrigin("VBContinuumSubtractor", "areFreqsInBounds"));
478 0 : LogIO os(LogOrigin("VBContinuumSubtractor", "apply"));
479 :
480 : os << LogIO::WARN
481 : << "Extrapolating to cover [" << 1.0e-9 * minfreq << ", "
482 : << 1.0e-9 * maxfreq << "] (GHz).\n"
483 : << "The frequency range used for the continuum fit was ["
484 0 : << 1.0e-9 * lofreq_p << ", "
485 0 : << 1.0e-9 * hifreq_p << "] (GHz)."
486 0 : << LogIO::POST;
487 0 : }
488 0 : return result;
489 : }
490 :
491 0 : Bool VBContinuumSubtractor::doShapesMatch(VisBuffer& vb,
492 : LogIO& os, const Bool squawk) const
493 : {
494 0 : Bool theydo = true;
495 :
496 0 : if(vb.nCorr() != static_cast<Int>(ncorr_p)){
497 0 : theydo = false;
498 0 : if(squawk)
499 : os << LogIO::SEVERE
500 0 : << "The supplied number of correlations, " << vb.nCorr()
501 0 : << ", does not match the expected " << ncorr_p
502 0 : << LogIO::POST;
503 : }
504 : // It's no longer the # of rows that matter but the maximum antenna #.
505 : // if(vb.nRow() != nrow_p){
506 0 : if(max(vb.antenna2()) > maxAnt_p){
507 0 : theydo = false; // Should it just flag unknown baselines?
508 0 : if(squawk)
509 : os << LogIO::SEVERE
510 0 : << "The fit is only valid for antennas with indices <= " << maxAnt_p
511 0 : << LogIO::POST;
512 : }
513 0 : return theydo;
514 : }
515 :
516 : // Do the subtraction
517 0 : Bool VBContinuumSubtractor::apply(VisBuffer& vb,
518 : const MS::PredefinedColumns whichcol,
519 : const Cube<Complex>& coeffs,
520 : const Cube<Bool>& coeffsOK,
521 : const Bool doSubtraction,
522 : const Bool squawk)
523 : {
524 : using casacore::operator*;
525 :
526 0 : LogIO os(LogOrigin("VBContinuumSubtractor", "apply"));
527 :
528 0 : if(!doShapesMatch(vb, os, squawk))
529 0 : return false;
530 :
531 0 : Bool ok = areFreqsInBounds(vb, squawk); // A Bool might be too Boolean here.
532 0 : ok = true; // Yep, returning false for a slight
533 : // extrapolation is too harsh.
534 :
535 0 : if(!(whichcol == MS::DATA || whichcol == MS::MODEL_DATA ||
536 : whichcol == MS::CORRECTED_DATA)){
537 0 : if(squawk)
538 : os << LogIO::SEVERE
539 : << MS::columnName(whichcol) << " is not supported."
540 0 : << LogIO::POST;
541 0 : return false;
542 : }
543 0 : Cube<Complex>& viscube(vb.dataCube(whichcol));
544 :
545 0 : uInt nchan = vb.nChannel();
546 0 : uInt nvbrow = vb.nRow();
547 :
548 : // DEBUGGING
549 : // os << LogIO::DEBUG1
550 : // << "nvbrow: " << nvbrow << ", nchan: " << nchan
551 : // << LogIO::POST;
552 : // // Check coeffs.
553 : // for(uInt vbrow = 0; vbrow < nvbrow; ++vbrow){
554 : // uInt blind = hashFunction(vb.antenna1()[vbrow],
555 : // vb.antenna2()[vbrow]);
556 :
557 : // for(uInt corrind = 0; corrind < ncorr_p; ++corrind){
558 : // if(coeffsOK(corrind, 0, blind)){
559 : // Complex cont = coeffs(corrind, 0, blind);
560 :
561 : // if(fabs(cont) < 0.001)
562 : // os << LogIO::WARN
563 : // << "cont(" << corrind << ", 0, " << blind << ") = "
564 : // << cont
565 : // << LogIO::POST;
566 : // }
567 : // }
568 : // }
569 : // END DEBUGGING
570 :
571 0 : Vector<Double> freqpow(fitorder_p + 1); // sf**ordind
572 0 : freqpow[0] = 1.0;
573 0 : Vector<Double>& freq(vb.frequency());
574 :
575 0 : for(uInt c = 0; c < nchan; ++c){
576 0 : Double sf = freqscale_p * (freq[c] - midfreq_p); // scaled frequency
577 :
578 0 : for(Int ordind = 1; ordind <= fitorder_p; ++ordind)
579 0 : freqpow[ordind] = sf * freqpow[ordind - 1];
580 :
581 0 : for(uInt vbrow = 0; vbrow < nvbrow; ++vbrow){
582 0 : uInt blind = hashFunction(vb.antenna1()[vbrow],
583 0 : vb.antenna2()[vbrow]);
584 :
585 0 : for(uInt corrind = 0; corrind < ncorr_p; ++corrind){
586 0 : if(coeffsOK(corrind, 0, blind)){
587 0 : Complex cont = coeffs(corrind, 0, blind);
588 :
589 0 : for(Int ordind = 1; ordind <= fitorder_p; ++ordind)
590 0 : cont += coeffs(corrind, ordind, blind) * freqpow[ordind];
591 0 : if(doSubtraction)
592 0 : viscube(corrind, c, vbrow) -= cont;
593 : else
594 0 : viscube(corrind, c, vbrow) = cont;
595 :
596 : // TODO: Adjust WEIGHT_SPECTRUM (create if necessary?), WEIGHT, and
597 : // SIGMA.
598 : }
599 : else
600 : {
601 : // jagonzal: Do not flag data in case of not successful
602 : // sfit because it might be cause by the line-free mask
603 : // and we may end up flagging 'good' line channels
604 0 : if (!tvi_debug) vb.flagCube()(corrind, c, vbrow) = true;
605 : // jagonzal: When returning the continuum, in case of not
606 : // successful fit we should return 0, not the input data
607 0 : if (tvi_debug and !doSubtraction) viscube(corrind, c, vbrow) = 0;
608 : }
609 : //vb.flag()(c, vbrow) = true;
610 : }
611 : }
612 : }
613 0 : return ok;
614 0 : }
615 :
616 :
617 : // Do the subtraction VB2 version
618 0 : Bool VBContinuumSubtractor::apply2(VisBuffer2& vb,
619 : Cube<Complex>& Vout,
620 : //const MS::PredefinedColumns whichcol,
621 : const Cube<Complex>& coeffs,
622 : const Cube<Bool>& coeffsOK,
623 : const Bool doSubtraction,
624 : const Bool squawk)
625 : {
626 : using casacore::operator*;
627 :
628 0 : LogIO os(LogOrigin("VBContinuumSubtractor", "apply2"));
629 :
630 0 : if(!doShapesMatch(vb, os, squawk))
631 0 : return false;
632 :
633 0 : Bool ok = areFreqsInBounds(vb, squawk); // A Bool might be too Boolean here.
634 0 : ok = true; // Yep, returning false for a slight
635 : // extrapolation is too harsh.
636 :
637 0 : Cube<Complex>& viscube(Vout);
638 :
639 0 : Cube<Bool> fC;
640 0 : fC.reference(vb.flagCube());
641 :
642 0 : uInt nchan = vb.nChannels();
643 0 : uInt nvbrow = vb.nRows();
644 :
645 0 : Vector<Double> freqpow(fitorder_p + 1); // sf**ordind
646 0 : freqpow[0] = 1.0;
647 0 : const Vector<Double>& freq(vb.getFrequencies(0));
648 :
649 0 : for(uInt c = 0; c < nchan; ++c){
650 0 : Double sf = freqscale_p * (freq[c] - midfreq_p); // scaled frequency
651 :
652 0 : for(Int ordind = 1; ordind <= fitorder_p; ++ordind)
653 0 : freqpow[ordind] = sf * freqpow[ordind - 1];
654 :
655 0 : for(uInt vbrow = 0; vbrow < nvbrow; ++vbrow){
656 0 : uInt blind = hashFunction(vb.antenna1()(vbrow),
657 0 : vb.antenna2()(vbrow));
658 :
659 0 : for(uInt corrind = 0; corrind < ncorr_p; ++corrind){
660 0 : if(coeffsOK(corrind, 0, blind)){
661 0 : Complex cont = coeffs(corrind, 0, blind);
662 :
663 0 : for(Int ordind = 1; ordind <= fitorder_p; ++ordind)
664 0 : cont += coeffs(corrind, ordind, blind) * freqpow[ordind];
665 0 : if(doSubtraction)
666 0 : viscube(corrind, c, vbrow) -= cont;
667 : else
668 0 : viscube(corrind, c, vbrow) = cont;
669 :
670 : // TODO: Adjust WEIGHT_SPECTRUM (create if necessary?), WEIGHT, and
671 : // SIGMA.
672 : }
673 : else
674 : {
675 : // jagonzal: Do not flag data in case of not successful
676 : // sfit because it might be cause by the line-free mask
677 : // and we may end up flagging 'good' line channels
678 0 : if (!tvi_debug) fC(corrind, c, vbrow) = true;
679 : // jagonzal: When returning the continuum, in case of not
680 : // successful fit we should return 0, not the input data
681 0 : if (tvi_debug and !doSubtraction) viscube(corrind, c, vbrow) = 0;
682 : }
683 : //vb.flag()(c, vbrow) = true;
684 : }
685 : }
686 : }
687 0 : return ok;
688 0 : }
689 :
690 0 : Bool VBContinuumSubtractor::doShapesMatch(VisBuffer2& vb,
691 : LogIO& os, const Bool squawk) const
692 : {
693 0 : Bool theydo = true;
694 :
695 0 : if(vb.nCorrelations() != static_cast<Int>(ncorr_p)){
696 0 : theydo = false;
697 0 : if(squawk)
698 : os << LogIO::SEVERE
699 0 : << "The supplied number of correlations, " << vb.nCorrelations()
700 0 : << ", does not match the expected " << ncorr_p
701 0 : << LogIO::POST;
702 : }
703 0 : if(max(vb.antenna2()) > maxAnt_p){
704 0 : theydo = false; // Should it just flag unknown baselines?
705 0 : if(squawk)
706 : os << LogIO::SEVERE
707 0 : << "The fit is only valid for antennas with indices <= " << maxAnt_p
708 0 : << LogIO::POST;
709 : }
710 0 : return theydo;
711 : }
712 :
713 0 : void VBContinuumSubtractor::getMinMaxFreq(VisBuffer2& vb,
714 : Double& minfreq,
715 : Double& maxfreq,
716 : const Bool initialize)
717 : {
718 0 : const Vector<Double>& freq(vb.getFrequencies(0)); // row=0, assumed uniform in VB2
719 0 : Int hichan = vb.nChannels() - 1;
720 0 : Int lochan = 0;
721 :
722 0 : if(initialize){
723 0 : maxfreq = -1.0;
724 0 : minfreq = DBL_MAX;
725 : }
726 :
727 0 : if(freq[hichan] < freq[lochan]){
728 0 : lochan = hichan;
729 0 : hichan = 0;
730 : }
731 0 : if(freq[hichan] > maxfreq)
732 0 : maxfreq = freq[hichan];
733 0 : if(freq[lochan] < minfreq)
734 0 : minfreq = freq[lochan];
735 0 : }
736 :
737 :
738 0 : Bool VBContinuumSubtractor::areFreqsInBounds(VisBuffer2& vb,
739 : const Bool squawk) const
740 : {
741 : Double maxfreq, minfreq;
742 :
743 0 : getMinMaxFreq(vb, minfreq, maxfreq);
744 0 : Bool result = minfreq >= lofreq_p && maxfreq <= hifreq_p;
745 :
746 0 : if(squawk && !result){
747 : // The truth was considered too alarming (CAS-1968).
748 : // LogIO os(LogOrigin("VBContinuumSubtractor", "areFreqsInBounds"));
749 0 : LogIO os(LogOrigin("VBContinuumSubtractor", "apply"));
750 :
751 : os << LogIO::WARN
752 : << "Extrapolating to cover [" << 1.0e-9 * minfreq << ", "
753 : << 1.0e-9 * maxfreq << "] (GHz).\n"
754 : << "The frequency range used for the continuum fit was ["
755 0 : << 1.0e-9 * lofreq_p << ", "
756 0 : << 1.0e-9 * hifreq_p << "] (GHz)."
757 0 : << LogIO::POST;
758 0 : }
759 0 : return result;
760 : }
761 :
762 :
763 :
764 : } //# NAMESPACE CASA - END
765 :
|
