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: aips2-request@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 :
30 : #include <casacore/casa/Utilities/Assert.h>
31 : #include <singledish/SingleDish/BLParameterParser.h>
32 :
33 : using namespace std;
34 :
35 : #define _ORIGIN LogOrigin("BLParameterParser", __func__, WHERE)
36 :
37 : using namespace casacore;
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 : namespace casa {
45 :
46 91 : BLParameterParser::BLParameterParser(string const file_name)
47 : {
48 91 : initialize();
49 91 : parse(file_name);
50 91 : blparam_file_ = file_name;
51 91 : }
52 :
53 91 : BLParameterParser::~BLParameterParser()
54 : {
55 91 : Clearup();
56 91 : }
57 :
58 91 : void BLParameterParser::Clearup()
59 : {
60 91 : if (!bl_parameters_.empty()) {
61 : map<const pair<size_t, size_t> , BLParameterSet*>::iterator
62 81 : iter = bl_parameters_.begin();
63 382 : while (iter != bl_parameters_.end()) {
64 301 : delete (*iter).second;
65 301 : ++iter;
66 : }
67 81 : bl_parameters_.clear();
68 : }
69 91 : }
70 :
71 91 : void BLParameterParser::initialize()
72 : {
73 91 : if (!bl_parameters_.empty()) Clearup();
74 91 : baseline_types_.resize(0);
75 : // initialize max orders
76 91 : size_t num_type = BaselineType_kNumElements;
77 455 : for (size_t i=0; i<num_type; ++i) {
78 364 : max_orders_[i] = 0;
79 : }
80 91 : }
81 :
82 228 : uint16_t BLParameterParser::get_max_order(size_t const bltype)
83 : {
84 :
85 420 : for (size_t i = 0; i < baseline_types_.size(); ++i) {
86 420 : if (bltype == baseline_types_[i]) {
87 : //bltype is in file
88 228 : 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 91 : void BLParameterParser::parse(string const file_name)
96 : {
97 273 : LogIO os(_ORIGIN);
98 182 : ifstream ifs(file_name.c_str(),ifstream::in);
99 91 : AlwaysAssert(!ifs.fail(), AipsError);
100 182 : string linestr;
101 524 : 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 433 : if (linestr.empty() || linestr[0]=='#') continue;
108 304 : BLParameterSet *bl_param = new BLParameterSet();
109 : size_t row_idx, pol_idx;
110 304 : ConvertLineToParam(linestr, row_idx, pol_idx, *bl_param);
111 304 : 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 << "]"
116 : << LogIO::POST;
117 : bl_param->PrintSummary();
118 : }
119 : // update bealine_types_ list
120 304 : size_t curr_type_idx = static_cast<size_t>(bl_param->baseline_type);
121 304 : bool new_type = true;
122 537 : for (size_t i=0; i<baseline_types_.size(); ++i){
123 329 : if (bl_param->baseline_type==baseline_types_[i]){
124 96 : new_type = false;
125 96 : break;
126 : }
127 : }
128 304 : if (new_type) baseline_types_.push_back(bl_param->baseline_type);
129 : // update max_orders_
130 304 : size_t curr_order = GetTypeOrder(*bl_param);
131 304 : if (curr_order > max_orders_[curr_type_idx])
132 202 : max_orders_[curr_type_idx] = curr_order;
133 : }
134 91 : }
135 :
136 304 : void BLParameterParser::SplitLine(string const &linestr,
137 : char const separator,
138 : vector<string> &strvec)
139 : {
140 608 : istringstream iss(linestr);
141 608 : string selem;
142 304 : size_t ielem(0);
143 4694 : while(getline(iss, selem, separator)) {
144 4390 : if (ielem >= strvec.size()) {
145 156 : strvec.resize(ielem+1);
146 : }
147 4390 : strvec[ielem] = selem;
148 4390 : ++ielem;
149 : }
150 304 : }
151 :
152 304 : 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 912 : std::vector<string> svec(BLParameters_kNumElements,"");
157 304 : SplitLine(linestr,',', svec);
158 : // parse mandatory data
159 304 : rowid = ConvertString<size_t>(svec[BLParameters_kRow]);
160 304 : polid = ConvertString<size_t>(svec[BLParameters_kPol]);
161 304 : paramset.baseline_mask = svec[BLParameters_kMask];
162 304 : size_t num_piece=USHRT_MAX;//SIZE_MAX;
163 608 : string const bltype_str = svec[BLParameters_kBaselineType];
164 304 : 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 238 : else if (bltype_str == "sinusoid")
176 : {
177 : // sinusoid is not supported yet
178 0 : throw(AipsError("Unsupported baseline type, sinusoid"));
179 : }
180 : else
181 : { // poly or chebyshev
182 238 : if (svec[BLParameters_kOrder].size()==0)
183 0 : throw(AipsError("Baseline type 'poly' and 'chebyshev' require order value."));
184 238 : paramset.baseline_type = bltype_str == "chebyshev" ?
185 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kChebyshev) :
186 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kPolynomial);
187 238 : paramset.order = ConvertString<uint16_t>(svec[BLParameters_kOrder]);
188 : }
189 : // parse clipping parameters
190 304 : if (svec[BLParameters_kNumIteration].size() == 0)
191 0 : throw(AipsError("Number of maximum clip iteration is mandatory"));
192 : paramset.num_fitting_max
193 304 : = ConvertString<uint16_t>(svec[BLParameters_kNumIteration]) + 1;
194 304 : if (svec[BLParameters_kClipThreshold].size()>0)
195 : paramset.clip_threshold_sigma
196 304 : = ConvertString<float>(svec[BLParameters_kClipThreshold]);
197 : // parse line finder parameter
198 304 : LineFinderParameter &lf_param = paramset.line_finder;
199 304 : lf_param.use_line_finder = svec[BLParameters_kLineFinder]=="true" ? true : false;;
200 304 : if (lf_param.use_line_finder)
201 : { // use line finder
202 11 : if (svec[BLParameters_kLFThreshold].size()>0)
203 : {
204 : lf_param.threshold
205 11 : = ConvertString<float>(svec[BLParameters_kLFThreshold]);
206 : }
207 22 : vector<size_t> edge(2,0);
208 11 : if (svec[BLParameters_kLeftEdge].size() > 0)
209 11 : lf_param.edge[0] = ConvertString<size_t>(svec[BLParameters_kLeftEdge]);
210 11 : if (svec[BLParameters_kRightEdge].size() > 0)
211 11 : lf_param.edge[1] = ConvertString<size_t>(svec[BLParameters_kRightEdge]);
212 11 : if (svec[BLParameters_kChanAverageLim].size()>0)
213 : {
214 : lf_param.chan_avg_limit
215 11 : = ConvertString<size_t>(svec[BLParameters_kChanAverageLim]);
216 : }
217 : }
218 304 : }
219 :
220 304 : uint16_t BLParameterParser::GetTypeOrder(BLParameterSet const &bl_param)
221 : {
222 304 : size_t const type = static_cast<size_t>(bl_param.baseline_type);
223 304 : switch (type)
224 : {
225 238 : case BaselineType_kPolynomial:
226 : case BaselineType_kChebyshev:
227 238 : return bl_param.order;
228 : break;
229 66 : case BaselineType_kCubicSpline:
230 66 : AlwaysAssert(bl_param.npiece<=USHRT_MAX, AipsError);//UINT16_MAX);
231 66 : return static_cast<uint16_t>(bl_param.npiece);
232 : break;
233 : // case BaselineType_kSinusoidal:
234 : // return static_cast<size_t>(bl_param.nwave.size()); <== must be max of nwave elements
235 : // break;
236 0 : default:
237 0 : throw(AipsError("Unsupported baseline type."));
238 : }
239 : }
240 :
241 1914 : bool BLParameterParser::GetFitParameter(size_t const rowid,size_t const polid, BLParameterSet &bl_param)
242 : {
243 : map<const pair<size_t, size_t> ,BLParameterSet*>::iterator
244 1914 : iter = bl_parameters_.begin();
245 1914 : iter = bl_parameters_.find(make_pair(rowid, polid));
246 1914 : if (iter==bl_parameters_.end()) {
247 : // no matching element
248 1673 : return false;
249 : }
250 241 : bl_param = *(*iter).second;
251 241 : return true;
252 : }
253 :
254 10 : BLTableParser::BLTableParser(string const file_name) : BLParameterParser(file_name)
255 : {
256 10 : initialize();
257 10 : blparam_file_ = file_name;
258 10 : bt_ = new BaselineTable(file_name);
259 10 : parse();
260 10 : }
261 :
262 10 : BLTableParser::~BLTableParser()
263 : {
264 10 : delete bt_;
265 10 : bt_ = 0;
266 10 : }
267 :
268 10 : void BLTableParser::initialize()
269 : {
270 10 : baseline_types_.resize(0);
271 : // initialize max orders
272 10 : size_t num_type = BaselineType_kNumElements;
273 50 : for (size_t i = 0; i < num_type; ++i) {
274 40 : max_orders_[i] = 0;
275 : }
276 10 : }
277 :
278 61 : uint16_t BLTableParser::GetTypeOrder(size_t const &baseline_type,
279 : uInt const irow, uInt const ipol)
280 : {
281 61 : switch (baseline_type)
282 : {
283 51 : case BaselineType_kPolynomial:
284 : case BaselineType_kChebyshev:
285 : {
286 51 : return static_cast<uint16_t>(bt_->getFPar(irow, ipol));
287 : break;
288 : }
289 10 : case BaselineType_kCubicSpline:
290 : {
291 10 : uInt npiece = bt_->getFPar(irow, ipol);
292 10 : AlwaysAssert(npiece <= USHRT_MAX, AipsError);//UINT16_MAX);
293 10 : return static_cast<uint16_t>(npiece);
294 : break;
295 : }
296 : // case BaselineType_kSinusoidal:
297 : // return static_cast<size_t>(nwave.size());
298 : // break;
299 0 : default:
300 0 : throw(AipsError("Unsupported baseline type."));
301 : }
302 : }
303 :
304 10 : void BLTableParser::parse()
305 : {
306 10 : uInt const npol = 2;
307 10 : size_t nrow = bt_->nrow();
308 20 : std::map<string, std::map<double, uInt> > sorted;
309 41 : for (uInt irow = 0; irow < nrow; ++irow) {
310 62 : stringstream ss;
311 31 : ss << bt_->getSpw(irow) << ","
312 31 : << bt_->getAntenna(irow) << ","
313 31 : << bt_->getBeam(irow);
314 31 : sorted[ss.str()][bt_->getTime(irow)] = irow;
315 93 : for (uInt ipol = 0; ipol < npol; ++ipol) {
316 62 : if (!bt_->getApply(irow, ipol)) continue;
317 : LIBSAKURA_SYMBOL(LSQFitType) curr_type_idx =
318 61 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(bt_->getBaselineType(irow, ipol));
319 61 : bool new_type = true;
320 95 : for (size_t i = 0; i < baseline_types_.size(); ++i){
321 75 : if (curr_type_idx == baseline_types_[i]){
322 41 : new_type = false;
323 41 : break;
324 : }
325 : }
326 61 : if (new_type) baseline_types_.push_back(curr_type_idx);
327 : // update max_orders_
328 61 : size_t curr_order = GetTypeOrder(curr_type_idx, irow, ipol);
329 61 : if (curr_order > max_orders_[curr_type_idx]) {
330 26 : max_orders_[curr_type_idx] = curr_order;
331 : }
332 : }
333 : }
334 :
335 41 : for (std::map<string, std::map<double, uInt> >::iterator it = sorted.begin(); it != sorted.end(); ++it) {
336 31 : sortedTimes_[it->first].clear();
337 31 : timeSortedIdx_[it->first].clear();
338 31 : numRows_[it->first] = 0;
339 31 : size_t num_rows = 0;
340 62 : for (std::map<double, uInt>::iterator it0 = sorted[it->first].begin(); it0 != sorted[it->first].end(); ++it0) {
341 31 : sortedTimes_[it->first].push_back(it0->first);
342 31 : timeSortedIdx_[it->first].push_back(it0->second);
343 31 : num_rows++;
344 : }
345 31 : numRows_[it->first] = num_rows;
346 : }
347 10 : }
348 :
349 32 : bool BLTableParser::GetFitParameterIdx(double const time, double const interval,
350 : size_t const scanid, size_t const beamid,
351 : size_t const antid, size_t const spwid,
352 : size_t &idx)
353 : {
354 32 : bool found = false;
355 64 : stringstream ss;
356 32 : ss << spwid << "," << antid << "," << beamid;
357 64 : string key = ss.str();
358 32 : if (numRows_.count(key) == 0) {
359 1 : return false;
360 : }
361 31 : uInt idx_top = 0;
362 31 : uInt idx_end = numRows_[key] - 1;
363 31 : uInt idx_mid = (idx_top + idx_end)/2;
364 31 : double time_top = sortedTimes_[key][idx_top];
365 31 : double time_end = sortedTimes_[key][idx_end];
366 31 : double time_mid = sortedTimes_[key][idx_mid];
367 31 : if ((time < time_top)||(time_end < time)) { // out of range.
368 0 : return false;
369 : }
370 31 : size_t idx0 = 0;
371 : //binary search in time-sorted table.
372 31 : while (idx_end - idx_top > 1) {
373 0 : if (time_top == time) {
374 0 : idx0 = idx_top;
375 0 : found = true;
376 0 : break;
377 0 : } else if (time_mid == time) {
378 0 : idx0 = idx_mid;
379 0 : found = true;
380 0 : break;
381 0 : } else if (time_end == time) {
382 0 : idx0 = idx_end;
383 0 : found = true;
384 0 : break;
385 0 : } else if (time_mid < time) {
386 0 : idx_top = idx_mid;
387 0 : time_top = time_mid;
388 : } else {
389 0 : idx_end = idx_mid;
390 0 : time_end = time_mid;
391 : }
392 0 : idx_mid = (idx_top + idx_end)/2;
393 0 : time_mid = sortedTimes_[key][idx_mid];
394 : }
395 31 : if (!found) {
396 : //only one or two candidates should exist now.
397 : //check if there is one with time identical to the specified time
398 31 : for (size_t i = idx_top; i <= idx_end; ++i) {
399 62 : if ((time - interval < sortedTimes_[key][i])&&
400 31 : (sortedTimes_[key][i] < time + interval)) {
401 31 : idx0 = i;
402 31 : found = true;
403 31 : break;
404 : }
405 : }
406 31 : if (!found) {
407 0 : return false;
408 : }
409 : }
410 : //finally, validate the candidate using its ids.
411 31 : idx = timeSortedIdx_[key][idx0];
412 31 : found = (scanid == bt_->getScan(idx))&&
413 62 : (beamid == bt_->getBeam(idx))&&
414 31 : (spwid == bt_->getSpw(idx));
415 :
416 31 : return found;
417 : }
418 :
419 : /*
420 : //simply search from the first line. not using binary search.
421 : bool BLTableParser::GetFitParameterIdx(double const time, double const interval,
422 : size_t const scanid, size_t const beamid,
423 : size_t const spwid, size_t &idx)
424 : {
425 : bool found = false;
426 : uInt idx_end = bt_->nrow() - 1;
427 : double bt_time;
428 : uInt bt_scanid;
429 : uInt bt_beamid;
430 : uInt bt_spwid;
431 : for (uInt i = 0; i <= idx_end; ++i) {
432 : bt_time = bt_->getTime(i);
433 : bt_scanid = bt_->getScan(i);
434 : bt_beamid = bt_->getBeam(i);
435 : bt_spwid = bt_->getSpw(i);
436 : if ((bt_time-interval < time) && (time < bt_time+interval) &&
437 : (scanid == bt_scanid) && (beamid == bt_beamid) &&
438 : (spwid == bt_spwid)) {
439 : idx = i;
440 : found = true;
441 : break;
442 : }
443 : }
444 : return found;
445 : }
446 : */
447 :
448 62 : void BLTableParser::GetFitParameterByIdx(size_t const idx, size_t const ipol,
449 : bool &apply,
450 : Vector<double> &coeff,
451 : Vector<size_t> &boundary,
452 : std::vector<bool> &masklist,
453 : BLParameterSet &bl_param)
454 : {
455 62 : apply = bt_->getApply(idx, ipol);
456 62 : if (!apply) return;
457 61 : size_t const type = bt_->getBaselineType(idx, ipol);
458 61 : bl_param.baseline_type = type;
459 122 : Vector<Int> fpar(bt_->getFuncParam(idx)[0]);
460 61 : switch (type) {
461 51 : case BaselineType_kPolynomial:
462 : case BaselineType_kChebyshev:
463 51 : bl_param.order = fpar[ipol];
464 51 : coeff.resize(bl_param.order + 1);
465 226 : for (size_t i = 0; i < coeff.size(); ++i) {
466 175 : coeff[i] = bt_->getResult(idx).row(ipol)[i];
467 : }
468 51 : break;
469 10 : case BaselineType_kCubicSpline:
470 10 : bl_param.npiece = fpar[ipol];
471 10 : boundary.resize(bl_param.npiece + 1);
472 30 : for (size_t i = 0; i < boundary.size(); ++i) {
473 20 : boundary[i] = bt_->getFuncFParam(idx).row(ipol)[i];
474 : }
475 10 : coeff.resize(bl_param.npiece * 4);
476 50 : for (size_t i = 0; i < coeff.size(); ++i) {
477 40 : coeff[i] = bt_->getResult(idx).row(ipol)[i];
478 : }
479 10 : break;
480 0 : default:
481 0 : throw(AipsError("Unsupported baseline type."));
482 : }
483 61 : masklist = bt_->getMask(idx, ipol);
484 : }
485 :
486 :
487 : } //# NAMESPACE CASA - END
|
