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 0 : BLParameterParser::BLParameterParser(string const file_name)
47 : {
48 0 : initialize();
49 0 : parse(file_name);
50 0 : blparam_file_ = file_name;
51 0 : }
52 :
53 0 : BLParameterParser::~BLParameterParser()
54 : {
55 0 : Clearup();
56 0 : }
57 :
58 0 : void BLParameterParser::Clearup()
59 : {
60 0 : if (!bl_parameters_.empty()) {
61 : map<const pair<size_t, size_t> , BLParameterSet*>::iterator
62 0 : iter = bl_parameters_.begin();
63 0 : while (iter != bl_parameters_.end()) {
64 0 : delete (*iter).second;
65 0 : ++iter;
66 : }
67 0 : bl_parameters_.clear();
68 : }
69 0 : }
70 :
71 0 : void BLParameterParser::initialize()
72 : {
73 0 : if (!bl_parameters_.empty()) Clearup();
74 0 : baseline_types_.resize(0);
75 : // initialize max orders
76 0 : size_t num_type = BaselineType_kNumElements;
77 0 : for (size_t i=0; i<num_type; ++i) {
78 0 : max_orders_[i] = 0;
79 : }
80 0 : }
81 :
82 0 : uint16_t BLParameterParser::get_max_order(size_t const bltype)
83 : {
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 << "]"
116 : << LogIO::POST;
117 : bl_param->PrintSummary();
118 : }
119 : // update bealine_types_ list
120 0 : size_t curr_type_idx = static_cast<size_t>(bl_param->baseline_type);
121 0 : bool new_type = true;
122 0 : for (size_t i=0; i<baseline_types_.size(); ++i){
123 0 : if (bl_param->baseline_type==baseline_types_[i]){
124 0 : new_type = false;
125 0 : break;
126 : }
127 : }
128 0 : if (new_type) baseline_types_.push_back(bl_param->baseline_type);
129 : // update max_orders_
130 0 : size_t curr_order = GetTypeOrder(*bl_param);
131 0 : if (curr_order > max_orders_[curr_type_idx])
132 0 : max_orders_[curr_type_idx] = curr_order;
133 : }
134 0 : }
135 :
136 0 : void BLParameterParser::SplitLine(string const &linestr,
137 : char const separator,
138 : vector<string> &strvec)
139 : {
140 0 : istringstream iss(linestr);
141 0 : string selem;
142 0 : size_t ielem(0);
143 0 : while(getline(iss, selem, separator)) {
144 0 : if (ielem >= strvec.size()) {
145 0 : strvec.resize(ielem+1);
146 : }
147 0 : strvec[ielem] = selem;
148 0 : ++ielem;
149 : }
150 0 : }
151 :
152 0 : 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 0 : std::vector<string> svec(BLParameters_kNumElements,"");
157 0 : SplitLine(linestr,',', svec);
158 : // parse mandatory data
159 0 : rowid = ConvertString<size_t>(svec[BLParameters_kRow]);
160 0 : polid = ConvertString<size_t>(svec[BLParameters_kPol]);
161 0 : paramset.baseline_mask = svec[BLParameters_kMask];
162 0 : size_t num_piece=USHRT_MAX;//SIZE_MAX;
163 0 : string const bltype_str = svec[BLParameters_kBaselineType];
164 0 : if (bltype_str == "cspline")
165 : {
166 0 : if (svec[BLParameters_kNPiece].size()==0)
167 0 : throw(AipsError("Baseline type 'cspline' requires num_piece value."));
168 0 : num_piece = ConvertString<size_t>(svec[BLParameters_kNPiece]);
169 : // CreateBaselineContext does not supoort elements number > max(uint16_t)
170 0 : if (num_piece > USHRT_MAX)
171 0 : throw(AipsError("num_piece is larger than max of uint16_t"));
172 0 : paramset.npiece = num_piece;
173 0 : paramset.baseline_type = static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kCubicSpline);
174 : }
175 0 : 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 0 : if (svec[BLParameters_kOrder].size()==0)
183 0 : throw(AipsError("Baseline type 'poly' and 'chebyshev' require order value."));
184 0 : paramset.baseline_type = bltype_str == "chebyshev" ?
185 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kChebyshev) :
186 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(BaselineType_kPolynomial);
187 0 : paramset.order = ConvertString<uint16_t>(svec[BLParameters_kOrder]);
188 : }
189 : // parse clipping parameters
190 0 : if (svec[BLParameters_kNumIteration].size() == 0)
191 0 : throw(AipsError("Number of maximum clip iteration is mandatory"));
192 : paramset.num_fitting_max
193 0 : = ConvertString<uint16_t>(svec[BLParameters_kNumIteration]) + 1;
194 0 : if (svec[BLParameters_kClipThreshold].size()>0)
195 : paramset.clip_threshold_sigma
196 0 : = ConvertString<float>(svec[BLParameters_kClipThreshold]);
197 : // parse line finder parameter
198 0 : LineFinderParameter &lf_param = paramset.line_finder;
199 0 : lf_param.use_line_finder = svec[BLParameters_kLineFinder]=="true" ? true : false;;
200 0 : if (lf_param.use_line_finder)
201 : { // use line finder
202 0 : if (svec[BLParameters_kLFThreshold].size()>0)
203 : {
204 : lf_param.threshold
205 0 : = ConvertString<float>(svec[BLParameters_kLFThreshold]);
206 : }
207 0 : vector<size_t> edge(2,0);
208 0 : if (svec[BLParameters_kLeftEdge].size() > 0)
209 0 : lf_param.edge[0] = ConvertString<size_t>(svec[BLParameters_kLeftEdge]);
210 0 : if (svec[BLParameters_kRightEdge].size() > 0)
211 0 : lf_param.edge[1] = ConvertString<size_t>(svec[BLParameters_kRightEdge]);
212 0 : if (svec[BLParameters_kChanAverageLim].size()>0)
213 : {
214 : lf_param.chan_avg_limit
215 0 : = ConvertString<size_t>(svec[BLParameters_kChanAverageLim]);
216 : }
217 : }
218 0 : }
219 :
220 0 : uint16_t BLParameterParser::GetTypeOrder(BLParameterSet const &bl_param)
221 : {
222 0 : size_t const type = static_cast<size_t>(bl_param.baseline_type);
223 0 : switch (type)
224 : {
225 0 : case BaselineType_kPolynomial:
226 : case BaselineType_kChebyshev:
227 0 : return bl_param.order;
228 : break;
229 0 : case BaselineType_kCubicSpline:
230 0 : AlwaysAssert(bl_param.npiece<=USHRT_MAX, AipsError);//UINT16_MAX);
231 0 : 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 0 : 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 0 : iter = bl_parameters_.begin();
245 0 : iter = bl_parameters_.find(make_pair(rowid, polid));
246 0 : if (iter==bl_parameters_.end()) {
247 : // no matching element
248 0 : return false;
249 : }
250 0 : bl_param = *(*iter).second;
251 0 : return true;
252 : }
253 :
254 0 : BLTableParser::BLTableParser(string const file_name) : BLParameterParser(file_name)
255 : {
256 0 : initialize();
257 0 : blparam_file_ = file_name;
258 0 : bt_ = new BaselineTable(file_name);
259 0 : parse();
260 0 : }
261 :
262 0 : BLTableParser::~BLTableParser()
263 : {
264 0 : delete bt_;
265 0 : bt_ = 0;
266 0 : }
267 :
268 0 : void BLTableParser::initialize()
269 : {
270 0 : baseline_types_.resize(0);
271 : // initialize max orders
272 0 : size_t num_type = BaselineType_kNumElements;
273 0 : for (size_t i = 0; i < num_type; ++i) {
274 0 : max_orders_[i] = 0;
275 : }
276 0 : }
277 :
278 0 : uint16_t BLTableParser::GetTypeOrder(size_t const &baseline_type,
279 : uInt const irow, uInt const ipol)
280 : {
281 0 : switch (baseline_type)
282 : {
283 0 : case BaselineType_kPolynomial:
284 : case BaselineType_kChebyshev:
285 : {
286 0 : return static_cast<uint16_t>(bt_->getFPar(irow, ipol));
287 : break;
288 : }
289 0 : case BaselineType_kCubicSpline:
290 : {
291 0 : uInt npiece = bt_->getFPar(irow, ipol);
292 0 : AlwaysAssert(npiece <= USHRT_MAX, AipsError);//UINT16_MAX);
293 0 : 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 0 : void BLTableParser::parse()
305 : {
306 0 : uInt const npol = 2;
307 0 : size_t nrow = bt_->nrow();
308 0 : std::map<string, std::map<double, uInt> > sorted;
309 0 : for (uInt irow = 0; irow < nrow; ++irow) {
310 0 : stringstream ss;
311 0 : ss << bt_->getSpw(irow) << ","
312 0 : << bt_->getAntenna(irow) << ","
313 0 : << bt_->getBeam(irow);
314 0 : sorted[ss.str()][bt_->getTime(irow)] = irow;
315 0 : for (uInt ipol = 0; ipol < npol; ++ipol) {
316 0 : if (!bt_->getApply(irow, ipol)) continue;
317 : LIBSAKURA_SYMBOL(LSQFitType) curr_type_idx =
318 0 : static_cast<LIBSAKURA_SYMBOL(LSQFitType)>(bt_->getBaselineType(irow, ipol));
319 0 : bool new_type = true;
320 0 : for (size_t i = 0; i < baseline_types_.size(); ++i){
321 0 : if (curr_type_idx == baseline_types_[i]){
322 0 : new_type = false;
323 0 : break;
324 : }
325 : }
326 0 : if (new_type) baseline_types_.push_back(curr_type_idx);
327 : // update max_orders_
328 0 : size_t curr_order = GetTypeOrder(curr_type_idx, irow, ipol);
329 0 : if (curr_order > max_orders_[curr_type_idx]) {
330 0 : max_orders_[curr_type_idx] = curr_order;
331 : }
332 : }
333 : }
334 :
335 0 : for (std::map<string, std::map<double, uInt> >::iterator it = sorted.begin(); it != sorted.end(); ++it) {
336 0 : sortedTimes_[it->first].clear();
337 0 : timeSortedIdx_[it->first].clear();
338 0 : numRows_[it->first] = 0;
339 0 : size_t num_rows = 0;
340 0 : for (std::map<double, uInt>::iterator it0 = sorted[it->first].begin(); it0 != sorted[it->first].end(); ++it0) {
341 0 : sortedTimes_[it->first].push_back(it0->first);
342 0 : timeSortedIdx_[it->first].push_back(it0->second);
343 0 : num_rows++;
344 : }
345 0 : numRows_[it->first] = num_rows;
346 : }
347 0 : }
348 :
349 0 : 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 0 : bool found = false;
355 0 : stringstream ss;
356 0 : ss << spwid << "," << antid << "," << beamid;
357 0 : string key = ss.str();
358 0 : if (numRows_.count(key) == 0) {
359 0 : return false;
360 : }
361 0 : uInt idx_top = 0;
362 0 : uInt idx_end = numRows_[key] - 1;
363 0 : uInt idx_mid = (idx_top + idx_end)/2;
364 0 : double time_top = sortedTimes_[key][idx_top];
365 0 : double time_end = sortedTimes_[key][idx_end];
366 0 : double time_mid = sortedTimes_[key][idx_mid];
367 0 : if ((time < time_top)||(time_end < time)) { // out of range.
368 0 : return false;
369 : }
370 0 : size_t idx0 = 0;
371 : //binary search in time-sorted table.
372 0 : 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 0 : 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 0 : for (size_t i = idx_top; i <= idx_end; ++i) {
399 0 : if ((time - interval < sortedTimes_[key][i])&&
400 0 : (sortedTimes_[key][i] < time + interval)) {
401 0 : idx0 = i;
402 0 : found = true;
403 0 : break;
404 : }
405 : }
406 0 : if (!found) {
407 0 : return false;
408 : }
409 : }
410 : //finally, validate the candidate using its ids.
411 0 : idx = timeSortedIdx_[key][idx0];
412 0 : found = (scanid == bt_->getScan(idx))&&
413 0 : (beamid == bt_->getBeam(idx))&&
414 0 : (spwid == bt_->getSpw(idx));
415 :
416 0 : 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 0 : 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 0 : apply = bt_->getApply(idx, ipol);
456 0 : if (!apply) return;
457 0 : size_t const type = bt_->getBaselineType(idx, ipol);
458 0 : bl_param.baseline_type = type;
459 0 : Vector<Int> fpar(bt_->getFuncParam(idx)[0]);
460 0 : switch (type) {
461 0 : case BaselineType_kPolynomial:
462 : case BaselineType_kChebyshev:
463 0 : bl_param.order = fpar[ipol];
464 0 : coeff.resize(bl_param.order + 1);
465 0 : for (size_t i = 0; i < coeff.size(); ++i) {
466 0 : coeff[i] = bt_->getResult(idx).row(ipol)[i];
467 : }
468 0 : break;
469 0 : case BaselineType_kCubicSpline:
470 0 : bl_param.npiece = fpar[ipol];
471 0 : boundary.resize(bl_param.npiece + 1);
472 0 : for (size_t i = 0; i < boundary.size(); ++i) {
473 0 : boundary[i] = bt_->getFuncFParam(idx).row(ipol)[i];
474 : }
475 0 : coeff.resize(bl_param.npiece * 4);
476 0 : for (size_t i = 0; i < coeff.size(); ++i) {
477 0 : coeff[i] = bt_->getResult(idx).row(ipol)[i];
478 : }
479 0 : break;
480 0 : default:
481 0 : throw(AipsError("Unsupported baseline type."));
482 : }
483 0 : masklist = bt_->getMask(idx, ipol);
484 : }
485 :
486 :
487 : } //# NAMESPACE CASA - END
|
