Line data Source code
1 : //# SingleDishMS.h: this defines a class that handles single dish MSes
2 : //#
3 : //# Copyright (C) 2014,2015
4 : //# National Astronomical Observatory of Japan
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: casa-feedback@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 : #ifndef _CASA_SINGLEDISH_MS_H_
29 : #define _CASA_SINGLEDISH_MS_H_
30 :
31 : #include <iostream>
32 : #include <string>
33 :
34 : #include <casacore/casa/aipstype.h>
35 : #include <casacore/casa/Containers/Record.h>
36 : #include <libsakura/sakura.h>
37 : #include <casacore/ms/MeasurementSets/MeasurementSet.h>
38 : #include <msvis/MSVis/VisBuffer2.h>
39 : #include <casacore/scimath/Mathematics/FFTServer.h>
40 : #include <singledish/SingleDish/SDMSManager.h>
41 :
42 : namespace casa { //# NAMESPACE CASA - BEGIN
43 : class SingleDishMS {
44 : public:
45 : // Default constructor
46 : SingleDishMS();
47 : // Construct from casacore::MS name string
48 : SingleDishMS(string const& ms_name);
49 :
50 : // Destructor
51 : ~SingleDishMS();
52 :
53 : /*
54 : * Return the name of the MeasurementSet
55 : */
56 : string name() const {
57 : return msname_;
58 : }
59 :
60 : bool close();
61 :
62 : /*
63 : * Formats selection parameters for single dish processing.
64 : * @param [in] selection A casacore::Record consists of selection key and values
65 : * @param [in] verbose If true, print summary of selection logger
66 : */
67 : void setSelection(casacore::Record const& selection, bool const verbose = true);
68 :
69 : /*
70 : * Set time averaging parameters.
71 : */
72 : void setAverage(casacore::Record const& average, bool const verbose = true);
73 :
74 : /*
75 : * Set polarization averaging parameters.
76 : */
77 : void setPolAverage(casacore::Record const& average, bool const verbose = true);
78 :
79 : // Invoke baseline subtraction
80 : // (polynomial, write results in new casacore::MS)
81 : void subtractBaseline(string const& in_column_name,
82 : string const& out_ms_name,
83 : string const& out_bloutput_name,
84 : bool const& do_subtract,
85 : string const& in_spw,
86 : bool const& update_weight,
87 : string const& sigma_value,
88 : string const& blfunc,
89 : int const order,
90 : float const clip_threshold_sigma,
91 : int const num_fitting_max,
92 : bool const linefinding,
93 : float const threshold,
94 : int const avg_limit,
95 : int const minwidth,
96 : std::vector<int> const& edge);
97 :
98 : //Cubicspline
99 : void subtractBaselineCspline(string const& in_column_name,
100 : string const& out_ms_name,
101 : string const& out_bloutput_name,
102 : bool const& do_subtract,
103 : string const& in_spw,
104 : bool const& update_weight,
105 : string const& sigma_value,
106 : int const npiece,
107 : float const clip_threshold_sigma,
108 : int const num_fitting_max,
109 : bool const linefinding,
110 : float const threshold,
111 : int const avg_limit,
112 : int const minwidth,
113 : std::vector<int> const& edge);
114 :
115 : //Sinusoid
116 : void subtractBaselineSinusoid(string const& in_column_name,
117 : string const& out_ms_name,
118 : string const& out_bloutput_name,
119 : bool const& do_subtract,
120 : string const& in_spw,
121 : bool const& update_weight,
122 : string const& sigma_value,
123 : string const& addwn0,
124 : string const& rejwn0,
125 : bool const applyfft,
126 : string const fftmethod,
127 : string const fftthresh,
128 : float const clip_threshold_sigma,
129 : int const num_fitting_max,
130 : bool const linefinding,
131 : float const threshold,
132 : int const avg_limit,
133 : int const minwidth,
134 : std::vector<int> const& edge);
135 :
136 : // variable fitting parameters stored in a text file
137 : void subtractBaselineVariable(string const& in_column_name,
138 : string const& out_ms_name,
139 : string const& out_bloutput_name,
140 : bool const& do_subtract,
141 : string const& in_spw,
142 : bool const& update_weight,
143 : string const& sigma_value,
144 : string const& param_file,
145 : bool const& verbose = true);
146 :
147 : // apply baseline table
148 : void applyBaselineTable(string const& in_column_name,
149 : string const& in_bltable_name,
150 : string const& in_spw,
151 : bool const& update_weight,
152 : string const& sigma_value,
153 : string const& out_ms_name);
154 :
155 : // fit line profile
156 : void fitLine(string const& in_column_name, string const& in_spw,
157 : string const& in_pol, string const& fitfunc, string const& in_nfit,
158 : bool const linefinding, float const threshold, int const avg_limit,
159 : int const minwidth, std::vector<int> const& edge,
160 : string const& tempfile_name, string const& temp_out_ms_name);
161 :
162 : // smooth data with arbitrary smoothing kernel
163 : // smoothing kernels currently supported include gaussian and boxcar
164 : void smooth(string const& kernelType, float const kernelWidth,
165 : string const& columnName, string const& outMsName);
166 :
167 : // C++ implementation of sdatmcor
168 : void atmcor(casacore::Record const &config, string const &columnName, string const &outMsName);
169 :
170 : private:
171 : /////////////////////////
172 : /// Utility functions ///
173 : /////////////////////////
174 : /*
175 : * Initializes member variables: in_column_
176 : */
177 : void initialize();
178 : /*
179 : * Formats selection parameters for single dish processing.
180 : * @param [in] selection A casacore::Record consists of selection key and values
181 : */
182 : void format_selection(casacore::Record& selection);
183 :
184 : // retrieve a field by name from casacore::Record as casa::String.
185 : casacore::String get_field_as_casa_string(casacore::Record const& in_data,
186 : string const& field_name);
187 :
188 : bool prepare_for_process(string const& in_column_name,
189 : string const& out_ms_name);
190 :
191 : bool prepare_for_process(string const& in_column_name,
192 : string const& out_ms_name,
193 : casacore::Block<casacore::Int> const& sortColumns,
194 : bool const addDefaultSortCols = false);
195 : void finalize_process();
196 :
197 : // check column 'in' is in input casacore::MS and set to 'out' if it exists.
198 : // if not, out is set to casacore::MS::UNDEFINED_COLUMN
199 : bool set_column(casacore::MSMainEnums::PredefinedColumns const& in,
200 : casacore::MSMainEnums::PredefinedColumns& out);
201 : // Convert a casacore::Complex casacore::Array to casacore::Float Array
202 : void convertArrayC2F(casacore::Array<casacore::Float>& from,
203 : casacore::Array<casacore::Complex> const& to);
204 : // Split a string with given delimiter
205 : std::vector<string> split_string(string const& s,
206 : char delim);
207 : // examine if a file with specified name exists
208 : bool file_exists(string const& filename);
209 : /* Convert spw selection string to vectors of spwid and channel
210 : ranges by parsing msseltoindex output. Also creates some
211 : placeholder vectors to store mask and the muber of channels
212 : in each selected SPW.
213 : [in] in_spw: input SPW selection string
214 : [out] spw : a vector of selected SPW IDs. the number of
215 : elements is the number of selected SPWs
216 : [out] chan : a vector of selected SPW IDs and channel ranges
217 : returned by msseltoindex. [[SPWID, start, end, stride], ...]
218 : [out] nchan : a vector of length spw.size() initialized by zero
219 : [out] mask : an uninitialized vector of length spw.size()
220 : [out] nchan_set: a vector of length spw.size() initilazed by false
221 : */
222 : void parse_spw(string const& in_spw,
223 : casacore::Vector<casacore::Int>& spw,
224 : casacore::Matrix<casacore::Int>& chan,
225 : casacore::Vector<size_t>& nchan,
226 : casacore::Vector<casacore::Vector<casacore::Bool> >& mask,
227 : casacore::Vector<bool>& nchan_set);
228 : /* Go through chunk and set valudes to nchan and mask selection
229 : vectors of the SPW if not already done.
230 : [in] rec_spw: a vector of selected SPW IDs. the number of
231 : elements is the number of selected SPWs
232 : [in] data_spw: a vector of SPW IDs in current chunk. the
233 : number of elements is equals to the number
234 : of rows in the chunk.
235 : [in] rec_chan: a vector of selected SPW IDs and channel ranges
236 : returned by msseltoindex. [[SPWID, start, end, stride], ...]
237 : [in] num_chan: a vector of the number of channels in current chunk.
238 : the number of elements is equals to the number
239 : of rows in the chunk.
240 : [out] nchan: a vector of length spw.size(). the number of channels
241 : in the corresponding SPW is set when the loop traverses
242 : the SPW for the first time.
243 : [out] mask: a vector of length spw.size().
244 : [in,out] nchan_set: a boolean vector of length spw.size().
245 : the value indicates if nchan, and mask of
246 : the corresponding SPW is already set.
247 : [out] new_nchan: returns true if this is the first time finding
248 : SPWs with the same number of channels.
249 : */
250 : void get_nchan_and_mask(casacore::Vector<casacore::Int> const &rec_spw,
251 : casacore::Vector<casacore::Int> const &data_spw,
252 : casacore::Matrix<casacore::Int> const &rec_chan,
253 : size_t const num_chan,
254 : casacore::Vector<size_t>& nchan,
255 : casacore::Vector<casacore::Vector<casacore::Bool> >& mask,
256 : casacore::Vector<bool>& nchan_set,
257 : bool& new_nchan);
258 : void get_mask_from_rec(casacore::Int spwid,
259 : casacore::Matrix<casacore::Int> const &rec_chan,
260 : casacore::Vector<casacore::Bool> &mask,
261 : bool initialize = true);
262 : void get_masklist_from_mask(size_t const num_chan,
263 : bool const* mask,
264 : casacore::Vector<casacore::uInt>& masklist);
265 : // Create a set of baseline contexts (if necessary)
266 : void get_baseline_context(size_t const bltype,
267 : uint16_t order,
268 : size_t num_chan,
269 : casacore::Vector<size_t> const& nchan,
270 : casacore::Vector<bool> const& nchan_set,
271 : casacore::Vector<size_t>& ctx_indices,
272 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *>& bl_contexts);
273 : void get_baseline_context(size_t const bltype,
274 : uint16_t order,
275 : size_t num_chan,
276 : size_t ispw,
277 : casacore::Vector<size_t> &ctx_indices,
278 : std::vector<size_t> &ctx_nchans,
279 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts);
280 : // Destroy a set of baseline contexts
281 : void destroy_baseline_contexts(std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts);
282 : void check_sakura_status(string const &name, LIBSAKURA_SYMBOL(Status) const status);
283 : void check_baseline_status(LIBSAKURA_SYMBOL(LSQFitStatus) const bl_status);
284 : template<typename T, typename U>
285 3183504 : void set_matrix_for_bltable(size_t const num_pol,
286 : size_t const num_data_max,
287 : std::vector<std::vector<T> > const& in_data,
288 : casacore::Array<U> &out_data) {
289 6377918 : for (size_t ipol = 0; ipol < num_pol; ++ipol) {
290 10622003 : for (size_t i = 0; i < num_data_max; ++i) {
291 7427589 : out_data[i][ipol] = static_cast<U>(0);
292 : }
293 3194414 : size_t num_data = in_data[ipol].size();
294 10621697 : for (size_t i = 0; i < num_data; ++i) {
295 7427283 : out_data[i][ipol] = static_cast<U>(in_data[ipol][i]);
296 : }
297 : }
298 3183504 : }
299 : template<typename T, typename U>
300 : void set_array_for_bltable(size_t const ipol,
301 : size_t const num_data,
302 : T const *in_data,
303 : casacore::Array<U> &out_data) {
304 : for (size_t i = 0; i < num_data; ++i) {
305 : out_data[i][ipol] = static_cast<U>(in_data[i]);
306 : }
307 : }
308 : size_t get_num_coeff_bloutput(size_t const bltype,
309 : size_t order,
310 : size_t& num_coeff_max);
311 : std::vector<int> string_to_list(string const &wn_str, char const delim);
312 : void get_effective_nwave(std::vector<int> const& addwn,
313 : std::vector<int> const& rejwn,
314 : int const wn_ulimit,
315 : std::vector<int>& effwn);
316 : void finalise_effective_nwave(std::vector<int> const& blparam_eff_base,
317 : std::vector<int> const& blparam_exclude,
318 : int const& blparam_upperlimit,
319 : size_t const& num_chan,
320 : float const* spec,
321 : bool const* mask,
322 : bool const& applyfft,
323 : string const& fftmethod,
324 : string const& fftthresh,
325 : std::vector<size_t>& blparam_eff);
326 : void parse_fftthresh(string const& fftthresh_str,
327 : string& fftthresh_attr,
328 : float& fftthresh_sigma,
329 : int& fftthresh_top);
330 : void select_wavenumbers_via_fft(size_t const num_chan,
331 : float const *spec,
332 : bool const *mask,
333 : string const &fftmethod,
334 : string const &fftthresh_attr,
335 : float const fftthresh_sigma,
336 : int const fftthresh_top,
337 : int const blparam_upperlimit,
338 : std::vector<int> &blparam_fft);
339 : void exec_fft(size_t const num_chan,
340 : float const *in_spec,
341 : bool const *in_mask,
342 : bool const get_real_imag,
343 : bool const get_ampl_only,
344 : std::vector<float> &fourier_spec);
345 : void interpolate_constant(int const num_chan,
346 : float const *in_spec,
347 : bool const *in_mask,
348 : casacore::Vector<casacore::Float> &spec);
349 : void merge_wavenumbers(std::vector<int> const& blparam_eff_base,
350 : std::vector<int> const& blparam_fft,
351 : std::vector<int> const& blparam_exclude,
352 : std::vector<size_t>& blparam_eff);
353 : list<pair<size_t, size_t>> findLineAndGetRanges(size_t const num_data,
354 : float const data[/*num_data*/],
355 : bool mask[/*num_data*/],
356 : float const threshold,
357 : int const avg_limit,
358 : int const minwidth,
359 : std::vector<int> const& edge,
360 : bool const invert);
361 : void findLineAndGetMask(size_t const num_data,
362 : float const data[/*num_data*/],
363 : bool const in_mask[/*num_data*/],
364 : float const threshold,
365 : int const avg_limit,
366 : int const minwidth,
367 : std::vector<int> const& edge,
368 : bool const invert,
369 : bool out_mask[/*num_data*/]);
370 : template<typename Func0, typename Func1, typename Func2, typename Func3>
371 : void doSubtractBaseline(string const& in_column_name,
372 : string const& out_ms_name,
373 : string const& out_bloutput_name,
374 : bool const& do_subtract,
375 : string const& in_spw,
376 : bool const& update_weight,
377 : string const& sigma_value,
378 : LIBSAKURA_SYMBOL(Status)& status,
379 : std::vector<LIBSAKURA_SYMBOL(LSQFitContextFloat) *> &bl_contexts,
380 : size_t const bltype,
381 : std::vector<int> const& blparam,
382 : std::vector<int> const& blparam_exclude,
383 : bool const& applyfft,
384 : string const& fftmethod,
385 : string const& fftthresh,
386 : float const clip_threshold_sigma,
387 : int const num_fitting_max,
388 : bool const linefinding,
389 : float const threshold,
390 : int const avg_limit,
391 : int const minwidth,
392 : std::vector<int> const& edge,
393 : Func0 func0,
394 : Func1 func1,
395 : Func2 func2,
396 : Func3 func3,
397 : casacore::LogIO os);
398 :
399 : /////////////////////////////
400 : /// casacore::MS handling functions ///
401 : /////////////////////////////
402 : // retrieve a spectrum at the row and plane (polarization) from data cube
403 : void get_spectrum_from_cube(casacore::Cube<casacore::Float>& data_cube,
404 : size_t const row,
405 : size_t const plane,
406 : size_t const num_data,
407 : float out_data[/*num_data*/]);
408 : // set a spectrum at the row and plane (polarization) to data cube
409 : void set_spectrum_to_cube(casacore::Cube<casacore::Float>& data_cube,
410 : size_t const row,
411 : size_t const plane,
412 : size_t const num_data,
413 : float in_data[/*num_data*/]);
414 : // get data cube (npol*nchan*nvirow) in in_column_ from visbuffer
415 : // and convert it to float cube
416 : void get_data_cube_float(vi::VisBuffer2 const& vb,
417 : casacore::Cube<casacore::Float>& data_cube);
418 : // get weight matrix (npol*nvirow) from visbuffer
419 : void get_weight_matrix(vi::VisBuffer2 const& vb,
420 : casacore::Matrix<casacore::Float>& weight_matrix);
421 : // set a weight at the row and plane (polarization) to weight matrix
422 : void set_weight_to_matrix(casacore::Matrix<casacore::Float>& weight_matrix,
423 : size_t const row,
424 : size_t const plane,
425 : float in_weight);
426 : // get flag cube (npol*nchan*nvirow) from visbuffer
427 : void get_flag_cube(vi::VisBuffer2 const& vb,
428 : casacore::Cube<casacore::Bool>& flag_cube);
429 : // retrieve a flag at the row and plane (polarization) from flag cube
430 : void get_flag_from_cube(casacore::Cube<casacore::Bool>& flag_cube,
431 : size_t const row,
432 : size_t const plane,
433 : size_t const num_flag,
434 : bool out_flag[/*num_flag*/]);
435 : // set a flag at the row and plane (polarization) to flag cube
436 : void set_flag_to_cube(casacore::Cube<casacore::Bool> &flag_cube,
437 : size_t const row,
438 : size_t const plane,
439 : size_t const num_flag,
440 : bool in_flag[/*num_data*/]);
441 : // flag all channels in a supectrum in cube at the row and plane (polarization)
442 : void flag_spectrum_in_cube(casacore::Cube<casacore::Bool> &flag_cube,
443 : size_t const row,
444 : size_t const plane);
445 : // return true if all channels are flagged
446 : bool allchannels_flagged(size_t const num_flag,
447 : bool const* flag);
448 : // returns the number of channels with true in input mask
449 : size_t NValidMask(size_t const num_mask,
450 : bool const* mask);
451 :
452 : /////////////////////////////////
453 : /// casacore::Array execution functions ///
454 : /////////////////////////////////
455 :
456 : ////////////////////////
457 : /// Member variables ///
458 : ////////////////////////
459 : // the name of input MS
460 : string msname_;
461 : // columns to read and save data
462 : casacore::MSMainEnums::PredefinedColumns in_column_; //, out_column_;
463 : // casacore::Record of selection
464 : casacore::Record selection_;
465 : // casacore::Record of average
466 : casacore::Record average_;
467 : // Record of polarization average
468 : casacore::Record pol_average_;
469 : // SDMSManager
470 : SDMSManager *sdh_;
471 : // pointer to accessor function
472 : void (*visCubeAccessor_)(vi::VisBuffer2 const &vb, casacore::Cube<casacore::Float> &cube);
473 : // smoothing flag
474 : casacore::Bool doSmoothing_;
475 : // offline ATM correction flag
476 : casacore::Bool doAtmCor_;
477 : casacore::Record atmCorConfig_;
478 :
479 : string bloutputname_csv;
480 : string bloutputname_text;
481 : string bloutputname_table;
482 :
483 : //split the name
484 : void split_bloutputname(string str);
485 :
486 : public:
487 : static bool importAsap(string const &infile, string const &outfile, bool const parallel=false);
488 : static bool importNRO(string const &infile, string const &outfile, bool const parallel=false);
489 : };
490 : // class SingleDishMS -END
491 :
492 : } //# NAMESPACE CASA - END
493 :
494 : #endif /* _CASA_SINGLEDISH_MS_H_ */
|
