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