Line data Source code
1 : //# UVContSubTVI.h: This file contains the interface definition of the MSTransformManager class.
2 : //#
3 : //# CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
4 : //# Copyright (C) Associated Universities, Inc. Washington DC, USA 2011, All rights reserved.
5 : //# Copyright (C) European Southern Observatory, 2011, All rights reserved.
6 : //#
7 : //# This library is free software; you can redistribute it and/or
8 : //# modify it under the terms of the GNU Lesser General Public
9 : //# License as published by the Free software Foundation; either
10 : //# version 2.1 of the License, or (at your option) any later version.
11 : //#
12 : //# This library is distributed in the hope that it will be useful,
13 : //# but WITHOUT ANY WARRANTY, without even the implied warranty of
14 : //# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 : //# Lesser General Public License for more details.
16 : //#
17 : //# You should have received a copy of the GNU Lesser General Public
18 : //# License along with this library; if not, write to the Free Software
19 : //# Foundation, Inc., 59 Temple Place, Suite 330, Boston,
20 : //# MA 02111-1307 USA
21 : //# $Id: $
22 :
23 : #ifndef UVContSubTVI_H_
24 : #define UVContSubTVI_H_
25 :
26 : #include <unordered_map>
27 :
28 : // Base class
29 : #include <mstransform/TVI/FreqAxisTVI.h>
30 :
31 : // Fitting classes
32 : #include <casacore/scimath/Fitting/LinearFitSVD.h>
33 : #include <casacore/scimath/Functionals/Polynomial.h>
34 : #include <mstransform/TVI/DenoisingLib.h>
35 : #include <mstransform/TVI/UVContSubResult.h>
36 :
37 : using namespace casacore;
38 :
39 : namespace casa { //# NAMESPACE CASA - BEGIN
40 :
41 : namespace vi { //# NAMESPACE VI - BEGIN
42 :
43 : // One spw:chan (as str) + fit_order (as int) from the input fitspec record
44 : typedef std::pair<std::string, unsigned int> InFitSpec;
45 :
46 : // field -> fitspec spec with spw:chan string and fitorder
47 : typedef std::unordered_map<int, std::vector<InFitSpec>> InFitSpecMap;
48 :
49 : // Specification of a fit (given per field, per spw, or globally)
50 : struct FitSpec {
51 373 : FitSpec() = default; // should be =delete when C++17/insert_or_assign is available
52 : FitSpec(Vector<bool> mask, unsigned int order);
53 :
54 : Vector<bool> lineFreeChannelMask;
55 : unsigned int fitOrder = -1;
56 : };
57 :
58 : //////////////////////////////////////////////////////////////////////////
59 : // UVContSubTVI class
60 : //////////////////////////////////////////////////////////////////////////
61 :
62 : class UVContSubTVI : public FreqAxisTVI
63 : {
64 :
65 : public:
66 :
67 : UVContSubTVI( ViImplementation2 * inputVii,
68 : const Record &configuration);
69 :
70 : ~UVContSubTVI();
71 :
72 : // Report the the ViImplementation type
73 44 : virtual String ViiType() const { return String("UVContSub( ")+getVii()->ViiType()+" )"; };
74 :
75 : virtual void floatData (Cube<float> & vis) const;
76 : virtual void visibilityObserved (Cube<Complex> & vis) const;
77 : virtual void visibilityCorrected (Cube<Complex> & vis) const;
78 : void savePrecalcModel(const Cube<Complex> & origVis,
79 : const Cube<Complex> & contSubVis) const;
80 : virtual void visibilityModel (Cube<Complex> & vis) const;
81 : virtual void result(casacore::Record &res) const;
82 :
83 : protected:
84 :
85 : InFitSpecMap parseConfiguration(const Record &configuration);
86 : void initialize(const InFitSpecMap &fitspec);
87 :
88 : template<class T> void transformDataCube( const Cube<T> &inputVis,
89 : const Cube<float> &inputWeight,
90 : Cube<T> &outputVis) const;
91 :
92 : template<class T> Complex transformDataCore(denoising::GslPolynomialModel<Double>* model,
93 : Vector<bool> *lineFreeChannelMask,
94 : const Cube<T> &inputVis,
95 : const Cube<bool> &inputFlags,
96 : const Cube<float> &inputWeight,
97 : Cube<T> &outputVis,
98 : int parallelCorrAxis=-1) const;
99 :
100 : private:
101 : rownr_t getMaxMSFieldID() const;
102 : InFitSpecMap parseFitSpec(const Record &configuration) const;
103 : InFitSpecMap parseDictFitSpec(const Record &configuration) const;
104 : void parseFieldSubDict(const Record &fieldRec, const std::vector<int> &fieldIdxs,
105 : InFitSpecMap &fitspec) const;
106 : void printInputFitSpec(const InFitSpecMap &fitspec) const;
107 : void populatePerFieldSpec(int fieldID, const std::vector<InFitSpec> &fitSpecs);
108 : void insertToFieldSpecMap(const std::vector<int> &fieldIdxs, const InFitSpec &spec,
109 : InFitSpecMap &fitspec) const;
110 : void fitSpecToPerFieldMap(const InFitSpecMap &fitspec);
111 : void spwInputChecks(const MeasurementSet *msVii, MSSelection &spwChan) const;
112 : unordered_map<int, vector<int>> makeLineFreeChannelSelMap(std::string spwChanStr) const;
113 :
114 : mutable uint fitOrder_p;
115 : mutable bool want_cont_p;
116 : // If the user gives SPWs in fitspec that are not in this list, give a warning
117 : std::string allowedSpws_p;
118 : mutable bool withDenoisingLib_p;
119 : mutable uint nThreads_p;
120 : mutable uint niter_p;
121 :
122 : // Maps field->SPW->(order, channel_mask), 1s will be or-ed with flags to exclude chans.
123 : unordered_map<int, unordered_map<int, FitSpec>> perFieldSpecMap_p;
124 : mutable map<int, unique_ptr<denoising::GslPolynomialModel<double>>> inputFrequencyMap_p;
125 :
126 : mutable UVContSubResult result_p;
127 :
128 : // cache model when writemodel / MODEL_DATA has to be produced
129 : mutable Cube<casacore::Complex> precalcModelVis_p;
130 : bool precalcModel_p = false;
131 : };
132 :
133 : //////////////////////////////////////////////////////////////////////////
134 : // UVContSubTVIFactory class
135 : //////////////////////////////////////////////////////////////////////////
136 :
137 : class UVContSubTVIFactory : public ViFactory
138 : {
139 :
140 : public:
141 :
142 : UVContSubTVIFactory(Record &configuration,ViImplementation2 *inputVII);
143 :
144 : protected:
145 :
146 : vi::ViImplementation2 * createVi (VisibilityIterator2 *) const;
147 : vi::ViImplementation2 * createVi () const;
148 :
149 : Record configuration_p;
150 : ViImplementation2 *inputVii_p;
151 : };
152 :
153 : //////////////////////////////////////////////////////////////////////////
154 : // UVContSubTVILayerFactory class (for _recursive_ layering)
155 : //////////////////////////////////////////////////////////////////////////
156 :
157 : class UVContSubTVILayerFactory : public ViiLayerFactory
158 : {
159 :
160 : public:
161 :
162 : UVContSubTVILayerFactory(Record &configuration);
163 108 : virtual ~UVContSubTVILayerFactory() {};
164 :
165 : protected:
166 :
167 : ViImplementation2 * createInstance(ViImplementation2* vii0) const;
168 :
169 : Record configuration_p;
170 :
171 : };
172 :
173 : //////////////////////////////////////////////////////////////////////////
174 : // UVContSubTransformEngine class
175 : //////////////////////////////////////////////////////////////////////////
176 :
177 : template<class T> class UVContSubKernel; // Forward declaration
178 :
179 : template<class T> class UVContSubTransformEngine : public FreqAxisTransformEngine2<T>
180 : {
181 : using FreqAxisTransformEngine2<T>::inputData_p;
182 : using FreqAxisTransformEngine2<T>::outputData_p;
183 : using FreqAxisTransformEngine2<T>::debug_p;
184 :
185 : public:
186 :
187 : UVContSubTransformEngine( UVContSubKernel<T> *kernel,
188 : DataCubeMap *inputData,
189 : DataCubeMap *outputData );
190 :
191 : void transform();
192 :
193 : void transformCore(DataCubeMap *inputData,DataCubeMap *outputData);
194 :
195 : protected:
196 :
197 : // This member has to be a pointer, otherwise there
198 : // are compile time problems due to the fact that
199 : // it is a pure virtual class.
200 : UVContSubKernel<T> *uvContSubKernel_p;
201 : };
202 :
203 : //////////////////////////////////////////////////////////////////////////
204 : // UVContSubKernel class
205 : //////////////////////////////////////////////////////////////////////////
206 :
207 : template<class T> class UVContSubKernel
208 : {
209 :
210 : public:
211 :
212 : UVContSubKernel( denoising::GslPolynomialModel<Double>* model,
213 : Vector<bool> *lineFreeChannelMask);
214 :
215 : virtual void kernel(DataCubeMap *inputData,
216 : DataCubeMap *outputData);
217 :
218 : virtual void changeFitOrder(size_t order) = 0;
219 :
220 : virtual void defaultKernel(Vector<T> &inputVector,
221 : Vector<T> &outputVector) = 0;
222 :
223 : virtual Complex kernelCore(Vector<T> &inputVector,
224 : Vector<bool> &inputFlags,
225 : Vector<float> &inputWeights,
226 : Vector<T> &outputVector) = 0;
227 :
228 3422 : Complex getChiSquared() { return chisq_p; }
229 69280 : void setDebug(bool debug) { debug_p = debug; }
230 :
231 : protected:
232 :
233 : bool debug_p;
234 : size_t fitOrder_p;
235 : denoising::GslPolynomialModel<Double> *model_p;
236 : Matrix<Double> freqPows_p;
237 : Vector<float> frequencies_p;
238 : Vector<bool> *lineFreeChannelMask_p;
239 : Complex chisq_p = Complex(std::numeric_limits<float>::infinity(),
240 : std::numeric_limits<float>::infinity());
241 : };
242 :
243 : //////////////////////////////////////////////////////////////////////////
244 : // UVContSubtractionKernel class
245 : //////////////////////////////////////////////////////////////////////////
246 :
247 : template<class T> class UVContSubtractionKernel : public UVContSubKernel<T>
248 : {
249 : using UVContSubKernel<T>::fitOrder_p;
250 : using UVContSubKernel<T>::model_p;
251 : using UVContSubKernel<T>::freqPows_p;
252 : using UVContSubKernel<T>::frequencies_p;
253 : using UVContSubKernel<T>::lineFreeChannelMask_p;
254 : using UVContSubKernel<T>::debug_p;
255 :
256 :
257 : public:
258 :
259 : UVContSubtractionKernel( denoising::GslPolynomialModel<Double>* model,
260 : Vector<bool> *lineFreeChannelMask=nullptr);
261 :
262 : void changeFitOrder(size_t order);
263 :
264 : void defaultKernel(Vector<Complex> &inputVector,
265 : Vector<Complex> &outputVector);
266 :
267 : void defaultKernel(Vector<float> &inputVector,
268 : Vector<float> &outputVector);
269 :
270 : Complex kernelCore(Vector<Complex> &inputVector,
271 : Vector<bool> &inputFlags,
272 : Vector<float> &inputWeights,
273 : Vector<Complex> &outputVector);
274 :
275 : Complex kernelCore(Vector<float> &inputVector,
276 : Vector<bool> &inputFlags,
277 : Vector<float> &inputWeights,
278 : Vector<float> &outputVector);
279 :
280 : private:
281 :
282 : LinearFitSVD<float> fitter_p;
283 : };
284 :
285 : //////////////////////////////////////////////////////////////////////////
286 : // UVContEstimationKernel class
287 : //////////////////////////////////////////////////////////////////////////
288 :
289 : template<class T> class UVContEstimationKernel : public UVContSubKernel<T>
290 : {
291 :
292 : using UVContSubKernel<T>::fitOrder_p;
293 : using UVContSubKernel<T>::model_p;
294 : using UVContSubKernel<T>::freqPows_p;
295 : using UVContSubKernel<T>::frequencies_p;
296 : using UVContSubKernel<T>::lineFreeChannelMask_p;
297 : using UVContSubKernel<T>::debug_p;
298 :
299 : public:
300 :
301 : UVContEstimationKernel( denoising::GslPolynomialModel<Double>* model,
302 : Vector<bool> *lineFreeChannelMask=nullptr);
303 :
304 : void changeFitOrder(size_t order);
305 :
306 : void defaultKernel(Vector<Complex> &inputVector,
307 : Vector<Complex> &outputVector);
308 :
309 : void defaultKernel(Vector<float> &inputVector,
310 : Vector<float> &outputVector);
311 :
312 : Complex kernelCore(Vector<Complex> &inputVector,
313 : Vector<bool> &inputFlags,
314 : Vector<float> &inputWeights,
315 : Vector<Complex> &outputVector);
316 :
317 : Complex kernelCore(Vector<float> &inputVector,
318 : Vector<bool> &inputFlags,
319 : Vector<float> &inputWeights,
320 : Vector<float> &outputVector);
321 :
322 : private:
323 :
324 : LinearFitSVD<float> fitter_p;
325 : };
326 :
327 : //////////////////////////////////////////////////////////////////////////
328 : // UVContSubtractionDenoisingKernel class
329 : //////////////////////////////////////////////////////////////////////////
330 :
331 : template<class T> class UVContSubtractionDenoisingKernel : public UVContSubKernel<T>
332 : {
333 : using UVContSubKernel<T>::fitOrder_p;
334 : using UVContSubKernel<T>::model_p;
335 : using UVContSubKernel<T>::freqPows_p;
336 : using UVContSubKernel<T>::frequencies_p;
337 : using UVContSubKernel<T>::lineFreeChannelMask_p;
338 : using UVContSubKernel<T>::debug_p;
339 :
340 :
341 : public:
342 :
343 : UVContSubtractionDenoisingKernel( denoising::GslPolynomialModel<Double>* model,
344 : size_t nIter,
345 : Vector<bool> *lineFreeChannelMask=nullptr);
346 :
347 : void changeFitOrder(size_t order);
348 :
349 : void defaultKernel(Vector<T> &inputVector,
350 : Vector<T> &outputVector);
351 :
352 : Complex kernelCore(Vector<T> &inputVector,
353 : Vector<bool> &inputFlags,
354 : Vector<float> &inputWeights,
355 : Vector<T> &outputVector);
356 :
357 : private:
358 :
359 : denoising::IterativelyReweightedLeastSquares fitter_p;
360 : };
361 :
362 : //////////////////////////////////////////////////////////////////////////
363 : // UVContEstimationDenoisingKernel class
364 : //////////////////////////////////////////////////////////////////////////
365 :
366 : template<class T> class UVContEstimationDenoisingKernel : public UVContSubKernel<T>
367 : {
368 :
369 : using UVContSubKernel<T>::fitOrder_p;
370 : using UVContSubKernel<T>::model_p;
371 : using UVContSubKernel<T>::freqPows_p;
372 : using UVContSubKernel<T>::frequencies_p;
373 : using UVContSubKernel<T>::lineFreeChannelMask_p;
374 : using UVContSubKernel<T>::debug_p;
375 :
376 : public:
377 :
378 : UVContEstimationDenoisingKernel( denoising::GslPolynomialModel<Double>* model,
379 : size_t nIter,
380 : Vector<bool> *lineFreeChannelMask=nullptr);
381 :
382 : void changeFitOrder(size_t order);
383 :
384 : void defaultKernel(Vector<T> &inputVector,
385 : Vector<T> &outputVector);
386 :
387 : Complex kernelCore(Vector<T> &inputVector,
388 : Vector<bool> &inputFlags,
389 : Vector<float> &inputWeights,
390 : Vector<T> &outputVector);
391 :
392 : private:
393 :
394 : denoising::IterativelyReweightedLeastSquares fitter_p;
395 : };
396 :
397 :
398 :
399 : } //# NAMESPACE VI - END
400 :
401 : } //# NAMESPACE CASA - END
402 :
403 : #endif /* UVContSubTVI_H_ */
404 :
|
