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 0 : BLParameterParser::BLParameterParser(string const file_name)
48 : {
49 0 : initialize();
50 0 : parse(file_name);
51 0 : blparam_file_ = file_name;
52 0 : }
53 :
54 0 : BLParameterParser::~BLParameterParser()
55 : {
56 0 : Clearup();
57 0 : }
58 :
59 0 : void BLParameterParser::Clearup()
60 : {
61 0 : if (!bl_parameters_.empty()) {
62 : map<const pair<size_t, size_t> , BLParameterSet*>::iterator
63 0 : iter = bl_parameters_.begin();
64 0 : while (iter != bl_parameters_.end()) {
65 0 : delete (*iter).second;
66 0 : ++iter;
67 : }
68 0 : bl_parameters_.clear();
69 : }
70 0 : }
71 :
72 0 : void BLParameterParser::initialize()
73 : {
74 0 : if (!bl_parameters_.empty()) Clearup();
75 0 : baseline_types_.resize(0);
76 : // initialize max orders
77 0 : size_t num_type = BaselineType_kNumElements;
78 0 : for (size_t i=0; i<num_type; ++i) {
79 0 : max_orders_[i] = 0;
80 : }
81 0 : }
82 :
83 0 : uint16_t BLParameterParser::get_max_order(size_t const bltype)
84 : {
85 0 : for (size_t i = 0; i < baseline_types_.size(); ++i) {
86 0 : if (bltype == baseline_types_[i]) {
87 : //bltype is in file
88 0 : 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 0 : void BLParameterParser::parse(string const file_name)
96 : {
97 0 : LogIO os(_ORIGIN);
98 0 : ifstream ifs(file_name.c_str(),ifstream::in);
99 0 : AlwaysAssert(!ifs.fail(), AipsError);
100 0 : string linestr;
101 0 : 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 0 : if (linestr.empty() || linestr[0]=='#') continue;
108 0 : BLParameterSet *bl_param = new BLParameterSet();
109 : size_t row_idx, pol_idx;
110 0 : ConvertLineToParam(linestr, row_idx, pol_idx, *bl_param);
111 0 : 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 0 : size_t curr_type_idx = static_cast<size_t>(bl_param->baseline_type);
120 0 : bool new_type = true;
121 0 : for (size_t i=0; i<baseline_types_.size(); ++i){
122 0 : if (bl_param->baseline_type==baseline_types_[i]){
123 0 : new_type = false;
124 0 : break;
125 : }
126 : }
127 0 : if (new_type) baseline_types_.push_back(bl_param->baseline_type);
128 : // update max_orders_
129 0 : size_t curr_order = GetTypeOrder(*bl_param);
130 0 : if (curr_order > max_orders_[curr_type_idx])
131 0 : max_orders_[curr_type_idx] = curr_order;
132 : }
133 0 : }
134 :
135 0 : void BLParameterParser::SplitLine(string const &linestr,
136 : char const separator,
137 : vector<string> &strvec)
138 : {
139 0 : istringstream iss(linestr);
140 0 : string selem;
141 0 : size_t ielem(0);
142 0 : while(getline(iss, selem, separator)) {
143 0 : if (ielem >= strvec.size()) {
144 0 : strvec.resize(ielem+1);
145 : }
146 0 : strvec[ielem] = selem;
147 0 : ++ielem;
148 : }
149 0 : }
150 :
151 0 : 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 0 : std::vector<string> svec(BLParameters_kNumElements,"");
156 0 : SplitLine(linestr,',', svec);
157 : // parse mandatory data
158 0 : rowid = ConvertString<size_t>(svec[BLParameters_kRow]);
159 0 : polid = ConvertString<size_t>(svec[BLParameters_kPol]);
160 0 : paramset.baseline_mask = svec[BLParameters_kMask];
161 0 : size_t num_piece=USHRT_MAX;//SIZE_MAX;
162 0 : string const bltype_str = svec[BLParameters_kBaselineType];
163 0 : if (bltype_str == "cspline")
164 : {
165 0 : if (svec[BLParameters_kNPiece].size()==0)
166 0 : throw(AipsError("Baseline type 'cspline' requires num_piece value."));
167 0 : num_piece = ConvertString<size_t>(svec[BLParameters_kNPiece]);
168 : // CreateBaselineContext does not supoort elements number > max(uint16_t)
169 0 : if (num_piece > USHRT_MAX)
170 0 : throw(AipsError("num_piece is larger than max of uint16_t"));
171 0 : paramset.npiece = num_piece;
172 0 : paramset.baseline_type = static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kCubicSpline);
173 : }
174 0 : else if (bltype_str == "sinusoid")
175 : {
176 : //Find the occurrence of "[n1,n2...]"
177 0 : std::regex pattern("\\[(\\d+(?:,\\d+)*)\\]");
178 0 : std::smatch matches;
179 0 : std::regex_search(linestr, matches, pattern);
180 0 : 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 0 : 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 0 : SplitLine(matches[1], ',', tmp_nwave);
189 0 : for(const auto& i : tmp_nwave)
190 0 : paramset.nwave.emplace_back(ConvertString<size_t>(i));
191 : //Sort just in case and erase duplicates
192 0 : std::sort(paramset.nwave.begin(), paramset.nwave.end());
193 0 : paramset.nwave.erase(std::unique(paramset.nwave.begin(),paramset.nwave.end()), paramset.nwave.end());
194 0 : paramset.baseline_type = static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kSinusoid);
195 0 : }
196 : else
197 : { // poly or chebyshev
198 0 : if (svec[BLParameters_kOrder].size()==0)
199 0 : throw(AipsError("Baseline type 'poly' and 'chebyshev' require order value."));
200 0 : paramset.baseline_type = bltype_str == "chebyshev" ?
201 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kChebyshev) :
202 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kPolynomial);
203 0 : paramset.order = ConvertString<uint16_t>(svec[BLParameters_kOrder]);
204 : }
205 : // parse clipping parameters
206 0 : if (svec[BLParameters_kNumIteration].size() == 0)
207 0 : throw(AipsError("Number of maximum clip iteration is mandatory"));
208 0 : paramset.num_fitting_max = ConvertString<uint16_t>(svec[BLParameters_kNumIteration]) + 1;
209 0 : if (svec[BLParameters_kClipThreshold].size()>0)
210 : paramset.clip_threshold_sigma
211 0 : = ConvertString<float>(svec[BLParameters_kClipThreshold]);
212 : // parse line finder parameter
213 0 : LineFinderParameter &lf_param = paramset.line_finder;
214 0 : lf_param.use_line_finder = svec[BLParameters_kLineFinder]=="true" ? true : false;;
215 0 : if (lf_param.use_line_finder)
216 : { // use line finder
217 0 : if (svec[BLParameters_kLFThreshold].size()>0)
218 : {
219 0 : lf_param.threshold = ConvertString<float>(svec[BLParameters_kLFThreshold]);
220 : }
221 0 : vector<size_t> edge(2,0);
222 0 : if (svec[BLParameters_kLeftEdge].size() > 0)
223 0 : lf_param.edge[0] = ConvertString<size_t>(svec[BLParameters_kLeftEdge]);
224 0 : if (svec[BLParameters_kRightEdge].size() > 0)
225 0 : lf_param.edge[1] = ConvertString<size_t>(svec[BLParameters_kRightEdge]);
226 0 : if (svec[BLParameters_kChanAverageLim].size()>0)
227 : {
228 0 : lf_param.chan_avg_limit = ConvertString<size_t>(svec[BLParameters_kChanAverageLim]);
229 : }
230 0 : }
231 0 : }
232 :
233 0 : uint16_t BLParameterParser::GetTypeOrder(BLParameterSet const &bl_param)
234 : {
235 0 : size_t const type = static_cast<size_t>(bl_param.baseline_type);
236 0 : switch (type)
237 : {
238 0 : case BaselineType_kPolynomial:
239 : case BaselineType_kChebyshev:
240 0 : return bl_param.order;
241 : break;
242 0 : case BaselineType_kCubicSpline:
243 0 : AlwaysAssert(bl_param.npiece<=USHRT_MAX, AipsError);//UINT16_MAX);
244 0 : return static_cast<uint16_t>(bl_param.npiece);
245 : break;
246 0 : case BaselineType_kSinusoid:
247 0 : 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 0 : 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 0 : iter = bl_parameters_.begin();
258 0 : iter = bl_parameters_.find(make_pair(rowid, polid));
259 0 : if (iter==bl_parameters_.end()) {
260 : // no matching element
261 0 : return false;
262 : }
263 0 : bl_param = *(*iter).second;
264 0 : return true;
265 : }
266 :
267 0 : BLTableParser::BLTableParser(string const file_name) : BLParameterParser(file_name)
268 : {
269 0 : initialize();
270 0 : blparam_file_ = file_name;
271 0 : bt_ = new BaselineTable(file_name);
272 0 : parse();
273 0 : }
274 :
275 0 : BLTableParser::~BLTableParser()
276 : {
277 0 : delete bt_;
278 0 : bt_ = 0;
279 0 : }
280 :
281 0 : void BLTableParser::initialize()
282 : {
283 0 : baseline_types_.resize(0);
284 : // initialize max orders
285 0 : size_t num_type = BaselineType_kNumElements;
286 0 : for (size_t i = 0; i < num_type; ++i) {
287 0 : max_orders_[i] = 0;
288 : }
289 0 : }
290 :
291 0 : uint16_t BLTableParser::GetTypeOrder(size_t const &baseline_type,
292 : uInt const irow, uInt const ipol)
293 : {
294 0 : switch (baseline_type)
295 : {
296 0 : case BaselineType_kPolynomial:
297 : case BaselineType_kChebyshev:
298 : {
299 0 : return static_cast<uint16_t>(bt_->getFPar(irow, ipol, baseline_type));
300 : break;
301 : }
302 0 : case BaselineType_kCubicSpline:
303 : {
304 0 : uInt npiece = bt_->getFPar(irow, ipol, baseline_type);
305 0 : AlwaysAssert(npiece <= USHRT_MAX, AipsError);//UINT16_MAX);
306 0 : return static_cast<uint16_t>(npiece);
307 : break;
308 : }
309 0 : case BaselineType_kSinusoid:
310 : {
311 : // Return the last element of nwave, hence the highest order
312 0 : 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 0 : void BLTableParser::parse()
321 : {
322 0 : uInt const npol = 2;
323 0 : size_t nrow = bt_->nrow();
324 0 : std::map<string, std::map<double, uInt> > sorted;
325 :
326 0 : for (uInt irow = 0; irow < nrow; ++irow) {
327 0 : stringstream ss;
328 0 : ss << bt_->getSpw(irow) << ","
329 0 : << bt_->getAntenna(irow) << ","
330 0 : << bt_->getBeam(irow);
331 0 : sorted[ss.str()][bt_->getTime(irow)] = irow;
332 0 : for (uInt ipol = 0; ipol < npol; ++ipol) {
333 0 : if (!bt_->getApply(irow, ipol)) continue;
334 :
335 : auto curr_type = static_cast<LIBSAKURA_SYMBOL(LSQFitType)>
336 0 : (bt_->getBaselineType(irow, ipol));
337 0 : bool new_type = true;
338 0 : for (size_t i = 0; i < baseline_types_.size(); ++i){
339 0 : if (curr_type == baseline_types_[i]){
340 0 : new_type = false;
341 0 : break;
342 : }
343 : }
344 :
345 0 : if (new_type) baseline_types_.push_back(curr_type);
346 : // update max_orders_
347 0 : size_t curr_order = GetTypeOrder(curr_type, irow, ipol);
348 0 : if (curr_order > max_orders_[curr_type]) {
349 0 : max_orders_[curr_type] = curr_order;
350 : }
351 : }
352 0 : }
353 :
354 0 : for (std::map<string, std::map<double, uInt> >::iterator it = sorted.begin(); it != sorted.end(); ++it) {
355 0 : sortedTimes_[it->first].clear();
356 0 : timeSortedIdx_[it->first].clear();
357 0 : numRows_[it->first] = 0;
358 0 : size_t num_rows = 0;
359 0 : for (std::map<double, uInt>::iterator it0 = sorted[it->first].begin(); it0 != sorted[it->first].end(); ++it0) {
360 0 : sortedTimes_[it->first].push_back(it0->first);
361 0 : timeSortedIdx_[it->first].push_back(it0->second);
362 0 : num_rows++;
363 : }
364 0 : numRows_[it->first] = num_rows;
365 : }
366 0 : }
367 :
368 0 : 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 0 : bool found = false;
374 0 : stringstream ss;
375 0 : ss << spwid << "," << antid << "," << beamid;
376 0 : string key = ss.str();
377 0 : if (numRows_.count(key) == 0) {
378 0 : return false;
379 : }
380 0 : uInt idx_top = 0;
381 0 : uInt idx_end = numRows_[key] - 1;
382 0 : uInt idx_mid = (idx_top + idx_end)/2;
383 :
384 0 : const auto & sorted_times_key = sortedTimes_[key];
385 :
386 0 : double time_top = sorted_times_key[idx_top];
387 0 : double time_end = sorted_times_key[idx_end];
388 0 : double time_mid = sorted_times_key[idx_mid];
389 0 : if ((time < time_top)||(time_end < time)) { // out of range.
390 0 : return false;
391 : }
392 0 : size_t idx0 = 0;
393 : //binary search in time-sorted table.
394 0 : 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 0 : 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 0 : for (size_t i = idx_top; i <= idx_end; ++i) {
421 0 : if ((time - interval / 2 < sorted_times_key[i]) &&
422 0 : (sorted_times_key[i] < time + interval / 2)) {
423 0 : idx0 = i;
424 0 : found = true;
425 0 : break;
426 : }
427 : }
428 0 : if (!found) {
429 0 : return false;
430 : }
431 : }
432 : //finally, validate the candidate using its ids.
433 0 : idx = timeSortedIdx_[key][idx0];
434 0 : found = (scanid == bt_->getScan(idx))&&
435 0 : (beamid == bt_->getBeam(idx))&&
436 0 : (spwid == bt_->getSpw(idx));
437 :
438 0 : return found;
439 0 : }
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 0 : 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 0 : apply = bt_->getApply(idx, ipol);
478 0 : if (!apply) return;
479 :
480 0 : size_t const type = bt_->getBaselineType(idx, ipol);
481 0 : bl_param.baseline_type = type;
482 :
483 0 : 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 0 : Vector<float> tmp_result(bt_->getResult(idx).row(ipol));
489 0 : size_t coeff_size = 0;
490 :
491 0 : switch (type) {
492 0 : case BaselineType_kPolynomial:
493 : case BaselineType_kChebyshev:
494 0 : bl_param.order = fpar[ipol];
495 0 : coeff_size = bl_param.order + 1;
496 0 : break;
497 :
498 0 : case BaselineType_kCubicSpline:
499 0 : bl_param.npiece = fpar[ipol];
500 0 : boundary.resize(bl_param.npiece + 1);
501 0 : convertArray(boundary, Vector<float>(bt_->getFuncFParam(idx).row(ipol)));
502 0 : coeff_size = bl_param.npiece * 4;
503 0 : break;
504 :
505 0 : case BaselineType_kSinusoid:
506 0 : bl_param.nwave = bt_->getFParSinusoid(idx, ipol);
507 0 : if (bl_param.nwave[0] == 0)
508 0 : coeff_size = bl_param.nwave.size() * 2 - 1;
509 : else
510 0 : coeff_size = bl_param.nwave.size() * 2;
511 0 : break;
512 0 : default:
513 0 : throw(AipsError("Unsupported baseline type."));
514 : }
515 : // Resize coeff to prepare array
516 0 : coeff.resize(coeff_size);
517 : // If there are 0s any place other than first
518 : // Cut off zeros in the coeff array
519 0 : if (tmp_result.size() > coeff_size)
520 0 : tmp_result.resize(coeff_size, true);
521 0 : convertArray(coeff, tmp_result);
522 :
523 0 : masklist = bt_->getMask(idx, ipol);
524 0 : }
525 :
526 : } //# NAMESPACE CASA - END
|
