Line data Source code
1 : //# BLParameterParser.cc: this code perses baseline fitting parameters
2 : //# Copyright (C) 2015
3 : //# National Astronomical Observatory of Japan
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 : #include <fstream>
28 : #include <iostream>
29 : #include <regex>
30 :
31 : #include <casacore/casa/Utilities/Assert.h>
32 : #include <singledish/SingleDish/BLParameterParser.h>
33 :
34 : using namespace std;
35 :
36 : #define _ORIGIN LogOrigin("BLParameterParser", __func__, WHERE)
37 :
38 : using namespace casacore;
39 : using namespace casacore;
40 : using namespace casacore;
41 : using namespace casacore;
42 : using namespace casacore;
43 : using namespace casacore;
44 : using namespace casacore;
45 : namespace casa {
46 :
47 102 : BLParameterParser::BLParameterParser(string const file_name)
48 : {
49 102 : initialize();
50 102 : parse(file_name);
51 102 : blparam_file_ = file_name;
52 102 : }
53 :
54 102 : BLParameterParser::~BLParameterParser()
55 : {
56 102 : Clearup();
57 102 : }
58 :
59 102 : void BLParameterParser::Clearup()
60 : {
61 102 : if (!bl_parameters_.empty()) {
62 : map<const pair<size_t, size_t> , BLParameterSet*>::iterator
63 89 : iter = bl_parameters_.begin();
64 430 : while (iter != bl_parameters_.end()) {
65 341 : delete (*iter).second;
66 341 : ++iter;
67 : }
68 89 : bl_parameters_.clear();
69 : }
70 102 : }
71 :
72 102 : void BLParameterParser::initialize()
73 : {
74 102 : if (!bl_parameters_.empty()) Clearup();
75 102 : baseline_types_.resize(0);
76 : // initialize max orders
77 102 : size_t num_type = BaselineType_kNumElements;
78 510 : for (size_t i=0; i<num_type; ++i) {
79 408 : max_orders_[i] = 0;
80 : }
81 102 : }
82 :
83 257 : uint16_t BLParameterParser::get_max_order(size_t const bltype)
84 : {
85 479 : for (size_t i = 0; i < baseline_types_.size(); ++i) {
86 479 : if (bltype == baseline_types_[i]) {
87 : //bltype is in file
88 257 : return max_orders_[bltype];
89 : }
90 : }
91 : // bltype is not in file
92 0 : throw(AipsError("The baseline type is not in file."));
93 : }
94 :
95 102 : void BLParameterParser::parse(string const file_name)
96 : {
97 204 : LogIO os(_ORIGIN);
98 102 : ifstream ifs(file_name.c_str(),ifstream::in);
99 102 : AlwaysAssert(!ifs.fail(), AipsError);
100 102 : string linestr;
101 587 : while (getline(ifs, linestr)) {
102 : // parse line here to construct BLParameterSet
103 : // The order should be
104 : // ROW,POL,MASK,NITERATION,CLIP_THRES,LF,LF_THRES,LEDGE,REDGE,CHANAVG,BL_TYPE,ORDER,N_PIECE,NWAVE
105 : // size_t,1,string,uint16_t,float,bool,float,size_t,size_t,size_t,sinusoidal,uint16_t,size_t,vector<size_t>
106 : //skip line starting with '#'
107 485 : if (linestr.empty() || linestr[0]=='#') continue;
108 350 : BLParameterSet *bl_param = new BLParameterSet();
109 : size_t row_idx, pol_idx;
110 350 : ConvertLineToParam(linestr, row_idx, pol_idx, *bl_param);
111 350 : bl_parameters_[make_pair(row_idx, pol_idx)] = bl_param;
112 : //Parameter summary output (Debugging purpose only)
113 : if (false) {
114 : os << "Summary of parsed Parameter" << LogIO::POST;
115 : os << "[ROW" << row_idx << ", POL" << pol_idx << "]" << LogIO::POST;
116 : bl_param->PrintSummary();
117 : }
118 : // update bealine_types_ list
119 350 : size_t curr_type_idx = static_cast<size_t>(bl_param->baseline_type);
120 350 : bool new_type = true;
121 637 : for (size_t i=0; i<baseline_types_.size(); ++i){
122 406 : if (bl_param->baseline_type==baseline_types_[i]){
123 119 : new_type = false;
124 119 : break;
125 : }
126 : }
127 350 : if (new_type) baseline_types_.push_back(bl_param->baseline_type);
128 : // update max_orders_
129 350 : size_t curr_order = GetTypeOrder(*bl_param);
130 350 : if (curr_order > max_orders_[curr_type_idx])
131 231 : max_orders_[curr_type_idx] = curr_order;
132 : }
133 102 : }
134 :
135 360 : void BLParameterParser::SplitLine(string const &linestr,
136 : char const separator,
137 : vector<string> &strvec)
138 : {
139 360 : istringstream iss(linestr);
140 360 : string selem;
141 360 : size_t ielem(0);
142 5442 : while(getline(iss, selem, separator)) {
143 5082 : if (ielem >= strvec.size()) {
144 210 : strvec.resize(ielem+1);
145 : }
146 5082 : strvec[ielem] = selem;
147 5082 : ++ielem;
148 : }
149 360 : }
150 :
151 350 : void BLParameterParser::ConvertLineToParam(string const &linestr,
152 : size_t &rowid, size_t &polid,
153 : BLParameterSet ¶mset){
154 : //split a line by ',' and make a vector of strings
155 700 : std::vector<string> svec(BLParameters_kNumElements,"");
156 350 : SplitLine(linestr,',', svec);
157 : // parse mandatory data
158 350 : rowid = ConvertString<size_t>(svec[BLParameters_kRow]);
159 350 : polid = ConvertString<size_t>(svec[BLParameters_kPol]);
160 350 : paramset.baseline_mask = svec[BLParameters_kMask];
161 350 : size_t num_piece=USHRT_MAX;//SIZE_MAX;
162 350 : string const bltype_str = svec[BLParameters_kBaselineType];
163 350 : if (bltype_str == "cspline")
164 : {
165 78 : if (svec[BLParameters_kNPiece].size()==0)
166 0 : throw(AipsError("Baseline type 'cspline' requires num_piece value."));
167 78 : num_piece = ConvertString<size_t>(svec[BLParameters_kNPiece]);
168 : // CreateBaselineContext does not supoort elements number > max(uint16_t)
169 78 : if (num_piece > USHRT_MAX)
170 0 : throw(AipsError("num_piece is larger than max of uint16_t"));
171 78 : paramset.npiece = num_piece;
172 78 : paramset.baseline_type = static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kCubicSpline);
173 : }
174 272 : else if (bltype_str == "sinusoid")
175 : {
176 : //Find the occurrence of "[n1,n2...]"
177 10 : std::regex pattern("\\[(\\d+(?:,\\d+)*)\\]");
178 10 : std::smatch matches;
179 10 : std::regex_search(linestr, matches, pattern);
180 10 : if (matches.size()< 1)
181 0 : throw(AipsError("Incorrect format or empty nwave list given: " +
182 0 : std::string(matches[0]) +
183 0 : ". Please refer the sdbaseline documentation. Ex. [1,2]"));
184 : // Split, convert and fill in the paramset_nwave
185 10 : std::vector<string> tmp_nwave;
186 : // regex result: matches[0] is the entire match
187 : // matches[1] is the content inside the square brackets
188 10 : SplitLine(matches[1], ',', tmp_nwave);
189 42 : for(const auto& i : tmp_nwave)
190 32 : paramset.nwave.emplace_back(ConvertString<size_t>(i));
191 : //Sort just in case and erase duplicates
192 10 : std::sort(paramset.nwave.begin(), paramset.nwave.end());
193 10 : paramset.nwave.erase(std::unique(paramset.nwave.begin(),paramset.nwave.end()), paramset.nwave.end());
194 10 : paramset.baseline_type = static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kSinusoid);
195 10 : }
196 : else
197 : { // poly or chebyshev
198 262 : if (svec[BLParameters_kOrder].size()==0)
199 0 : throw(AipsError("Baseline type 'poly' and 'chebyshev' require order value."));
200 262 : paramset.baseline_type = bltype_str == "chebyshev" ?
201 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kChebyshev) :
202 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kPolynomial);
203 262 : paramset.order = ConvertString<uint16_t>(svec[BLParameters_kOrder]);
204 : }
205 : // parse clipping parameters
206 350 : if (svec[BLParameters_kNumIteration].size() == 0)
207 0 : throw(AipsError("Number of maximum clip iteration is mandatory"));
208 350 : paramset.num_fitting_max = ConvertString<uint16_t>(svec[BLParameters_kNumIteration]) + 1;
209 350 : if (svec[BLParameters_kClipThreshold].size()>0)
210 : paramset.clip_threshold_sigma
211 340 : = ConvertString<float>(svec[BLParameters_kClipThreshold]);
212 : // parse line finder parameter
213 350 : LineFinderParameter &lf_param = paramset.line_finder;
214 350 : lf_param.use_line_finder = svec[BLParameters_kLineFinder]=="true" ? true : false;;
215 350 : if (lf_param.use_line_finder)
216 : { // use line finder
217 14 : if (svec[BLParameters_kLFThreshold].size()>0)
218 : {
219 14 : lf_param.threshold = ConvertString<float>(svec[BLParameters_kLFThreshold]);
220 : }
221 14 : vector<size_t> edge(2,0);
222 14 : if (svec[BLParameters_kLeftEdge].size() > 0)
223 14 : lf_param.edge[0] = ConvertString<size_t>(svec[BLParameters_kLeftEdge]);
224 14 : if (svec[BLParameters_kRightEdge].size() > 0)
225 14 : lf_param.edge[1] = ConvertString<size_t>(svec[BLParameters_kRightEdge]);
226 14 : if (svec[BLParameters_kChanAverageLim].size()>0)
227 : {
228 14 : lf_param.chan_avg_limit = ConvertString<size_t>(svec[BLParameters_kChanAverageLim]);
229 : }
230 14 : }
231 350 : }
232 :
233 350 : uint16_t BLParameterParser::GetTypeOrder(BLParameterSet const &bl_param)
234 : {
235 350 : size_t const type = static_cast<size_t>(bl_param.baseline_type);
236 350 : switch (type)
237 : {
238 262 : case BaselineType_kPolynomial:
239 : case BaselineType_kChebyshev:
240 262 : return bl_param.order;
241 : break;
242 78 : case BaselineType_kCubicSpline:
243 78 : AlwaysAssert(bl_param.npiece<=USHRT_MAX, AipsError);//UINT16_MAX);
244 78 : return static_cast<uint16_t>(bl_param.npiece);
245 : break;
246 10 : case BaselineType_kSinusoid:
247 10 : return static_cast<uint16_t>(bl_param.nwave[bl_param.nwave.size() - 1]); //<== must be max of nwave elements
248 : break;
249 0 : default:
250 0 : throw(AipsError("Unsupported baseline type."));
251 : }
252 : }
253 :
254 1960 : bool BLParameterParser::GetFitParameter(size_t const rowid,size_t const polid, BLParameterSet &bl_param)
255 : {
256 : map<const pair<size_t, size_t> ,BLParameterSet*>::iterator
257 1960 : iter = bl_parameters_.begin();
258 1960 : iter = bl_parameters_.find(make_pair(rowid, polid));
259 1960 : if (iter==bl_parameters_.end()) {
260 : // no matching element
261 1679 : return false;
262 : }
263 281 : bl_param = *(*iter).second;
264 281 : return true;
265 : }
266 :
267 13 : BLTableParser::BLTableParser(string const file_name) : BLParameterParser(file_name)
268 : {
269 13 : initialize();
270 13 : blparam_file_ = file_name;
271 13 : bt_ = new BaselineTable(file_name);
272 13 : parse();
273 13 : }
274 :
275 13 : BLTableParser::~BLTableParser()
276 : {
277 13 : delete bt_;
278 13 : bt_ = 0;
279 13 : }
280 :
281 13 : void BLTableParser::initialize()
282 : {
283 13 : baseline_types_.resize(0);
284 : // initialize max orders
285 13 : size_t num_type = BaselineType_kNumElements;
286 65 : for (size_t i = 0; i < num_type; ++i) {
287 52 : max_orders_[i] = 0;
288 : }
289 13 : }
290 :
291 73 : uint16_t BLTableParser::GetTypeOrder(size_t const &baseline_type,
292 : uInt const irow, uInt const ipol)
293 : {
294 73 : switch (baseline_type)
295 : {
296 51 : case BaselineType_kPolynomial:
297 : case BaselineType_kChebyshev:
298 : {
299 51 : return static_cast<uint16_t>(bt_->getFPar(irow, ipol, baseline_type));
300 : break;
301 : }
302 16 : case BaselineType_kCubicSpline:
303 : {
304 16 : uInt npiece = bt_->getFPar(irow, ipol, baseline_type);
305 16 : AlwaysAssert(npiece <= USHRT_MAX, AipsError);//UINT16_MAX);
306 16 : return static_cast<uint16_t>(npiece);
307 : break;
308 : }
309 6 : case BaselineType_kSinusoid:
310 : {
311 : // Return the last element of nwave, hence the highest order
312 6 : return static_cast<uint16_t>(bt_->getFPar(irow, ipol, baseline_type));
313 : break;
314 : }
315 0 : default:
316 0 : throw(AipsError("Unsupported baseline type."));
317 : }
318 : }
319 :
320 13 : void BLTableParser::parse()
321 : {
322 13 : uInt const npol = 2;
323 13 : size_t nrow = bt_->nrow();
324 13 : std::map<string, std::map<double, uInt> > sorted;
325 :
326 50 : for (uInt irow = 0; irow < nrow; ++irow) {
327 37 : stringstream ss;
328 37 : ss << bt_->getSpw(irow) << ","
329 37 : << bt_->getAntenna(irow) << ","
330 37 : << bt_->getBeam(irow);
331 37 : sorted[ss.str()][bt_->getTime(irow)] = irow;
332 111 : for (uInt ipol = 0; ipol < npol; ++ipol) {
333 74 : if (!bt_->getApply(irow, ipol)) continue;
334 :
335 : auto curr_type = static_cast<LIBSAKURA_SYMBOL(LSQFitType)>
336 73 : (bt_->getBaselineType(irow, ipol));
337 73 : bool new_type = true;
338 113 : for (size_t i = 0; i < baseline_types_.size(); ++i){
339 87 : if (curr_type == baseline_types_[i]){
340 47 : new_type = false;
341 47 : break;
342 : }
343 : }
344 :
345 73 : if (new_type) baseline_types_.push_back(curr_type);
346 : // update max_orders_
347 73 : size_t curr_order = GetTypeOrder(curr_type, irow, ipol);
348 73 : if (curr_order > max_orders_[curr_type]) {
349 33 : max_orders_[curr_type] = curr_order;
350 : }
351 : }
352 37 : }
353 :
354 50 : for (std::map<string, std::map<double, uInt> >::iterator it = sorted.begin(); it != sorted.end(); ++it) {
355 37 : sortedTimes_[it->first].clear();
356 37 : timeSortedIdx_[it->first].clear();
357 37 : numRows_[it->first] = 0;
358 37 : size_t num_rows = 0;
359 74 : for (std::map<double, uInt>::iterator it0 = sorted[it->first].begin(); it0 != sorted[it->first].end(); ++it0) {
360 37 : sortedTimes_[it->first].push_back(it0->first);
361 37 : timeSortedIdx_[it->first].push_back(it0->second);
362 37 : num_rows++;
363 : }
364 37 : numRows_[it->first] = num_rows;
365 : }
366 13 : }
367 :
368 41 : bool BLTableParser::GetFitParameterIdx(double const time, double const interval,
369 : size_t const scanid, size_t const beamid,
370 : size_t const antid, size_t const spwid,
371 : size_t &idx)
372 : {
373 41 : bool found = false;
374 41 : stringstream ss;
375 41 : ss << spwid << "," << antid << "," << beamid;
376 41 : string key = ss.str();
377 41 : if (numRows_.count(key) == 0) {
378 1 : return false;
379 : }
380 40 : uInt idx_top = 0;
381 40 : uInt idx_end = numRows_[key] - 1;
382 40 : uInt idx_mid = (idx_top + idx_end)/2;
383 :
384 40 : const auto & sorted_times_key = sortedTimes_[key];
385 :
386 40 : double time_top = sorted_times_key[idx_top];
387 40 : double time_end = sorted_times_key[idx_end];
388 40 : double time_mid = sorted_times_key[idx_mid];
389 40 : if ((time < time_top)||(time_end < time)) { // out of range.
390 3 : return false;
391 : }
392 37 : size_t idx0 = 0;
393 : //binary search in time-sorted table.
394 37 : while (idx_end - idx_top > 1) {
395 0 : if (time_top == time) {
396 0 : idx0 = idx_top;
397 0 : found = true;
398 0 : break;
399 0 : } else if (time_mid == time) {
400 0 : idx0 = idx_mid;
401 0 : found = true;
402 0 : break;
403 0 : } else if (time_end == time) {
404 0 : idx0 = idx_end;
405 0 : found = true;
406 0 : break;
407 0 : } else if (time_mid < time) {
408 0 : idx_top = idx_mid;
409 0 : time_top = time_mid;
410 : } else {
411 0 : idx_end = idx_mid;
412 0 : time_end = time_mid;
413 : }
414 0 : idx_mid = (idx_top + idx_end)/2;
415 0 : time_mid = sorted_times_key[idx_mid];
416 : }
417 37 : if (!found) {
418 : //only one or two candidates should exist now.
419 : //check if there is one with time identical to the specified time
420 37 : for (size_t i = idx_top; i <= idx_end; ++i) {
421 74 : if ((time - interval / 2 < sorted_times_key[i]) &&
422 37 : (sorted_times_key[i] < time + interval / 2)) {
423 37 : idx0 = i;
424 37 : found = true;
425 37 : break;
426 : }
427 : }
428 37 : if (!found) {
429 0 : return false;
430 : }
431 : }
432 : //finally, validate the candidate using its ids.
433 37 : idx = timeSortedIdx_[key][idx0];
434 37 : found = (scanid == bt_->getScan(idx))&&
435 74 : (beamid == bt_->getBeam(idx))&&
436 37 : (spwid == bt_->getSpw(idx));
437 :
438 37 : return found;
439 41 : }
440 :
441 : /*
442 : //simply search from the first line. not using binary search.
443 : bool BLTableParser::GetFitParameterIdx(double const time, double const interval,
444 : size_t const scanid, size_t const beamid,
445 : size_t const spwid, size_t &idx)
446 : {
447 : bool found = false;
448 : uInt idx_end = bt_->nrow() - 1;
449 : double bt_time;
450 : uInt bt_scanid;
451 : uInt bt_beamid;
452 : uInt bt_spwid;
453 : for (uInt i = 0; i <= idx_end; ++i) {
454 : bt_time = bt_->getTime(i);
455 : bt_scanid = bt_->getScan(i);
456 : bt_beamid = bt_->getBeam(i);
457 : bt_spwid = bt_->getSpw(i);
458 : if ((bt_time-interval < time) && (time < bt_time+interval) &&
459 : (scanid == bt_scanid) && (beamid == bt_beamid) &&
460 : (spwid == bt_spwid)) {
461 : idx = i;
462 : found = true;
463 : break;
464 : }
465 : }
466 : return found;
467 : }
468 : */
469 :
470 74 : void BLTableParser::GetFitParameterByIdx(size_t const idx, size_t const ipol,
471 : bool &apply,
472 : Vector<double> &coeff,
473 : Vector<size_t> &boundary,
474 : std::vector<bool> &masklist,
475 : BLParameterSet &bl_param)
476 : {
477 74 : apply = bt_->getApply(idx, ipol);
478 74 : if (!apply) return;
479 :
480 73 : size_t const type = bt_->getBaselineType(idx, ipol);
481 73 : bl_param.baseline_type = type;
482 :
483 146 : Vector<Int> fpar(bt_->getFuncParam(idx)[0]);
484 : /*
485 : The apply mode may have polarizations with different baseline
486 : To avoid the confusion and unify code, the coeff will be resized in all cases
487 : */
488 73 : Vector<float> tmp_result(bt_->getResult(idx).row(ipol));
489 73 : size_t coeff_size = 0;
490 :
491 73 : switch (type) {
492 51 : case BaselineType_kPolynomial:
493 : case BaselineType_kChebyshev:
494 51 : bl_param.order = fpar[ipol];
495 51 : coeff_size = bl_param.order + 1;
496 51 : break;
497 :
498 16 : case BaselineType_kCubicSpline:
499 16 : bl_param.npiece = fpar[ipol];
500 16 : boundary.resize(bl_param.npiece + 1);
501 16 : convertArray(boundary, Vector<float>(bt_->getFuncFParam(idx).row(ipol)));
502 16 : coeff_size = bl_param.npiece * 4;
503 16 : break;
504 :
505 6 : case BaselineType_kSinusoid:
506 6 : bl_param.nwave = bt_->getFParSinusoid(idx, ipol);
507 6 : if (bl_param.nwave[0] == 0)
508 0 : coeff_size = bl_param.nwave.size() * 2 - 1;
509 : else
510 6 : coeff_size = bl_param.nwave.size() * 2;
511 6 : break;
512 0 : default:
513 0 : throw(AipsError("Unsupported baseline type."));
514 : }
515 : // Resize coeff to prepare array
516 73 : coeff.resize(coeff_size);
517 : // If there are 0s any place other than first
518 : // Cut off zeros in the coeff array
519 73 : if (tmp_result.size() > coeff_size)
520 2 : tmp_result.resize(coeff_size, true);
521 73 : convertArray(coeff, tmp_result);
522 :
523 73 : masklist = bt_->getMask(idx, ipol);
524 73 : }
525 :
526 : } //# NAMESPACE CASA - END
|
