Line data Source code
1 : //# VBContinuumSubtractor.h: Fit a continuum model to a VisBuffer and provide
2 : //# the continuum and/or line estimates as VisBuffers.
3 : //# Copyright (C) 2011
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : //#
29 : #ifndef MSVIS_VBCONTINUUMSUBTRACTOR_H
30 : #define MSVIS_VBCONTINUUMSUBTRACTOR_H
31 :
32 : #include <casacore/casa/aips.h>
33 : #include <casacore/casa/Arrays/Cube.h>
34 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
35 : //#include <msvis/MSVis/VisBuffer.h>
36 : #include <msvis/MSVis/VisBuffer2.h>
37 :
38 :
39 : namespace casacore{
40 :
41 : class LogIO;
42 : }
43 :
44 : namespace casa { //# NAMESPACE CASA - BEGIN
45 :
46 : class VisBuffer;
47 : class VisBuffGroupAcc;
48 :
49 : // <summary>Fits and optionally subtracts the continuum in visibility spectra</summary>
50 : // <use visibility=export>
51 : //
52 : // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
53 : // </reviewed>
54 : //
55 : // <prerequisite>
56 : // <li> <linkto class=casacore::MeasurementSet>MeasurementSet</linkto>
57 : // </prerequisite>
58 : //
59 : // <etymology>
60 : // This class's main aim is to subtract the continuum from VisBuffers.
61 : // </etymology>
62 : //
63 : // <synopsis>
64 : // Spectral line observations often contain continuum emission which is
65 : // present in all channels (often with a small slope and curvature across the band). This
66 : // class fits a polynomial to this continuum, and can return the continuum
67 : // and/or line (continuum-subtracted) estimates as VisBuffers.
68 : // </synopsis>
69 : //
70 : // <example>
71 : // <srcBlock>
72 : // const casacore::Int fitorder = 4; // Careful! High orders might
73 : // // absorb power from the lines.
74 : // casacore::MS inMS(fileName);
75 : // casacore::Block<casacore::Int> sortcolumns;
76 : // ROVisIter vi(inMS, sortcolumns, 0.0);
77 : // VisBuffer cvb(vi); // Continuum estimate
78 : // VisBuffer lvb(vi); // Line estimate
79 : // VisBuffer vb(vi);
80 : // for(vi.originChunks(); vi.moreChunks(); vi.nextChunk()){
81 : // for(vi.origin(); vi.more(); ++vi){
82 : // VBContinuumSubtractor contestor(vb, fitorder); // Continuum Estimator
83 : //
84 : // contestor.cont_est(cvb); // Put continuum estimate into cvb.
85 : // contestor.cont_subtracted(lvb); // Put line estimate into lvb.
86 : //
87 : // // Do something with cvb and lvb...
88 : // }
89 : // }
90 : // </srcBlock>
91 : // </example>
92 : //
93 : // <motivation>
94 : // This class provides continuum fitting by polynomials with order >= 0, with
95 : // a convenient interface for use by other classes in synthesis and msvis
96 : // (i.e. calibration, imaging, and split).
97 : // </motivation>
98 : //
99 : // <todo asof="">
100 : // </todo>
101 : class VBContinuumSubtractor
102 : {
103 : public:
104 : VBContinuumSubtractor();
105 :
106 : // Construct it with the low and high frequencies used for scaling the
107 : // frequencies in the polynomial.
108 : VBContinuumSubtractor(const casacore::Double lofreq, const casacore::Double hifreq);
109 :
110 : ~VBContinuumSubtractor();
111 :
112 : // Set the # of correlations and fitorder from shp, the total number of input
113 : // channels to look at (including masked ones!), and the low and high scaling
114 : // frequencies.
115 : void init(const casacore::IPosition& shp, const casacore::uInt maxAnt, const casacore::uInt totnumchan,
116 : const casacore::Double lof, const casacore::Double hif);
117 :
118 : // Set the low and high frequencies, and #s of correlations, antennas, and
119 : // channels from vbga. Returns false if vbga is empty.
120 : casacore::Bool initFromVBGA(VisBuffGroupAcc& vbga);
121 :
122 : // Makes the continuum estimate by fitting a frequency polynomial of order
123 : // fitorder to the data in vbga. It sets the low and high frequencies used
124 : // for scaling the frequencies in the polynomial to the min and max
125 : // frequencies in vbga.
126 : // casacore::Input: vbga, The data
127 : // fitorder, e.g. 2 for a + bf + cf**2
128 : // doInit, if true call initFromVBGA(vbga)
129 : // doResize if true set coeffs and coeffsOK to the right shape.
130 : // Output (these will be resized):
131 : // coeffs: casacore::Polynomial coefficients for the continuum, indexed by (corr,
132 : // order, hash(ant1, ant2).
133 : // coeffsOK: and whether or not they're usable.
134 : void fit(VisBuffGroupAcc& vbga, const casacore::Int fitorder,
135 : casacore::MS::PredefinedColumns whichcol,
136 : casacore::Cube<casacore::Complex>& coeffs, casacore::Cube<casacore::Bool>& coeffsOK,
137 : const casacore::Bool doInit=false, const casacore::Bool doResize=false,
138 : const casacore::Bool squawk=true);
139 :
140 : // Apply the continuum estimate in coeffs (from fit) to vb. The affected
141 : // column of vb is chosen by whichcol, which must be exactly one of casacore::MS::DATA,
142 : // casacore::MS::MODEL_DATA, or casacore::MS::CORRECTED_DATA lest an exception will be thrown.
143 : // If doSubtraction is true, vb becomes the line estimate. Otherwise it is
144 : // replaced with the continuum estimate. Returns false if it detects an
145 : // error, true otherwise. squawk=true increases the potential number of
146 : // warnings sent to the logger.
147 : casacore::Bool apply(VisBuffer& vb,
148 : const casacore::MS::PredefinedColumns whichcol,
149 : const casacore::Cube<casacore::Complex>& coeffs,
150 : const casacore::Cube<casacore::Bool>& coeffsOK, const casacore::Bool doSubtraction=true,
151 : const casacore::Bool squawk=true);
152 :
153 : // Returns whether or not vb's frequencies are within the bounds used for the
154 : // continuum fit. If not, and squawk is true, a warning will be sent to the
155 : // logger. (Extrapolation is allowed, just not recommended.)
156 : casacore::Bool areFreqsInBounds(VisBuffer& vb, const casacore::Bool squawk) const;
157 : casacore::Bool areFreqsInBounds(vi::VisBuffer2& vb, const casacore::Bool squawk) const;
158 :
159 : // Fills minfreq and maxfreq with the minimum and maximum frequencies (in Hz,
160 : // acc. to the casacore::MS def'n v.2) of vb (not the continuum fit!). If you want to
161 : // get the extreme frequencies over a set of VisBuffers, set initialize to
162 : // false and initialize minfreq and maxfreq yourself to DBL_MAX and -1.0,
163 : // respectively.
164 : static void getMinMaxFreq(VisBuffer& vb, casacore::Double& minfreq, casacore::Double& maxfreq,
165 : const casacore::Bool initialize=true);
166 : static void getMinMaxFreq(vi::VisBuffer2& vb, casacore::Double& minfreq, casacore::Double& maxfreq,
167 : const casacore::Bool initialize=true);
168 :
169 : // These are provided for the calibration framework so that a
170 : // VBContinuumSubtractor can be c'ted from a VisBuffGroupAcc, make a
171 : // continuum fit, write its results to a caltable, destroyed, and then later
172 : // another one can be c'ted from the caltable and apply the fit to a
173 : // VisBuffer. Obviously it'd be safer and faster to use the same
174 : // VBContinuumSubtractor for fitting and application, but the rest of CASA
175 : // isn't ready for that yet (3/7/2011).
176 : casacore::Int getOrder() const {return fitorder_p;}
177 :
178 2 : casacore::Double getLowFreq() const {return lofreq_p;} // Lowest frequency used in the fit,
179 2 : casacore::Double getHighFreq() const {return hifreq_p;} // and highest, in Hz, acc. to
180 : // the casacore::MS def'n v.2.
181 : casacore::Int getMaxAntNum() const {return maxAnt_p;} // -1 if unready.
182 :
183 : // The total number of input channels that will be looked at (including
184 : // masked ones!)
185 2 : casacore::uInt getTotNumChan() const {return totnumchan_p;}
186 :
187 : // Low (lof) and high (hif) frequencies, in Hz, used for renormalizing
188 : // frequencies in the polynomials.
189 150 : void setScalingFreqs(casacore::Double lof, casacore::Double hif){
190 150 : lofreq_p = lof;
191 150 : hifreq_p = hif;
192 150 : midfreq_p = 0.5 * (lof + hif);
193 150 : freqscale_p = calcFreqScale();
194 150 : }
195 :
196 : // Set the maximum number of antennas (actually, 1 + the maximum antenna
197 : // number).
198 152 : void setNAnt(const casacore::uInt nAnt){
199 152 : maxAnt_p = nAnt - 1;
200 152 : nHashes_p = (nAnt * (nAnt + 1)) / 2; // Allows for autocorrs.
201 152 : }
202 :
203 : // Set the total number of input channels that will be looked at (including
204 : // masked ones!)
205 : void setTotNumChan(const casacore::uInt tnc) {totnumchan_p = tnc;}
206 :
207 : // A convenience function for prepping coeffs and coeffsOK to hold the
208 : // polynomial coefficients and their validities.
209 : void resize(casacore::Cube<casacore::Complex>& coeffs, casacore::Cube<casacore::Bool>& coeffsOK) const;
210 :
211 39 : casacore::Bool checkSize(casacore::Cube<casacore::Complex>& coeffs, casacore::Cube<casacore::Bool>& coeffsOK) const
212 : {
213 39 : return coeffs.shape()[0] == static_cast<casacore::Int>(ncorr_p) &&
214 39 : coeffs.shape()[1] == fitorder_p + 1 &&
215 39 : coeffs.shape()[2] == static_cast<casacore::Int>(nHashes_p) &&
216 39 : coeffsOK.shape()[0] == static_cast<casacore::Int>(ncorr_p) &&
217 117 : coeffsOK.shape()[1] == fitorder_p + 1 &&
218 78 : coeffsOK.shape()[2] == static_cast<casacore::Int>(nHashes_p);
219 : }
220 :
221 191 : casacore::Double calcFreqScale() const {
222 191 : return hifreq_p > midfreq_p ? 1.0 / (hifreq_p - midfreq_p) : 1.0;
223 : }
224 :
225 0 : void setTVIDebug(bool debug) {tvi_debug = debug;}
226 :
227 :
228 : casacore::Bool apply2(vi::VisBuffer2& vb,
229 : casacore::Cube<casacore::Complex>& Vout,
230 : //const casacore::MS::PredefinedColumns whichcol,
231 : const casacore::Cube<casacore::Complex>& coeffs,
232 : const casacore::Cube<casacore::Bool>& coeffsOK, const casacore::Bool doSubtraction=true,
233 : const casacore::Bool squawk=true);
234 :
235 :
236 :
237 : private:
238 : // Disable default copying, and assignment.
239 : VBContinuumSubtractor& operator=(VBContinuumSubtractor& other);
240 :
241 : // Can the fit be _applied_ to vb?
242 : // If not and squawk is true, send a severe message to os.
243 : casacore::Bool doShapesMatch(VisBuffer& vb, casacore::LogIO& os, const casacore::Bool squawk=true) const;
244 : casacore::Bool doShapesMatch(vi::VisBuffer2& vb, casacore::LogIO& os, const casacore::Bool squawk=true) const;
245 :
246 : // Compute baseline (row) index in coeffs_p for (ant1, ant2).
247 : // It ASSUMES that ant1 and ant2 are both <= maxAnt_p.
248 137280 : casacore::uInt hashFunction(const casacore::Int ant1, const casacore::Int ant2)
249 : {
250 137280 : return (maxAnt_p + 1) * ant1 - (ant1 * (ant1 - 1)) / 2 + ant2 - ant1;
251 : }
252 :
253 : //VisBuffGroupAcc& vbga_p; // Holds the VisBuffers
254 : casacore::Int fitorder_p; // Order of the fitting polynomial.
255 :
256 : casacore::Double lofreq_p; // Lowest frequency used in the continuum fit,
257 : casacore::Double hifreq_p; // and highest, in Hz, acc. to the casacore::MS def'n v.2.
258 : casacore::Double midfreq_p; // 0.5 * (lofreq_p + hifreq_p)
259 : casacore::Double freqscale_p;
260 : casacore::Int maxAnt_p; // Highest 0 based antenna number that can be
261 : // used with the fits. -1 if not ready.
262 : casacore::uInt nHashes_p; // Calculated and cached from maxAnt_p.
263 : casacore::uInt ncorr_p;
264 : casacore::uInt totnumchan_p;
265 :
266 : casacore::PtrBlock<casacore::Vector<casacore::Bool> * > chanmask_p;
267 :
268 : bool tvi_debug;
269 : };
270 :
271 : } //# NAMESPACE CASA - END
272 :
273 : #endif
|
