source: trunk/FACT++/src/fitsdump.cc@ 17165

Last change on this file since 17165 was 17147, checked in by tbretz, 11 years ago
Added more options to format the output stream. Added filter support to the stats option and fixed some issues with the calculation.
File size: 33.1 KB
Line 
1//****************************************************************
2/** @class FitsDumper
3
4 @brief Dumps contents of fits tables to stdout or a file
5
6 */
7 //****************************************************************
8#include "Configuration.h"
9
10#include <float.h>
11
12#include <map>
13#include <fstream>
14
15#include <boost/regex.hpp>
16
17#include "tools.h"
18#include "Time.h"
19#include "externals/factfits.h"
20
21#ifdef HAVE_ROOT
22#include "TFormula.h"
23#endif
24
25using namespace std;
26
27struct MyColumn
28{
29 string name;
30
31 fits::Table::Column col;
32
33 uint32_t first;
34 uint32_t last;
35
36 void *ptr;
37};
38
39struct minMaxStruct
40{
41 double min;
42 double max;
43 long double average;
44 long double squared;
45 long numValues;
46 minMaxStruct() : min(FLT_MAX), max(-FLT_MAX), average(0), squared(0), numValues(0) { }
47
48 void add(long double val)
49 {
50 average += val;
51 squared += val*val;
52
53 if (val<min)
54 min = val;
55
56 if (val>max)
57 max = val;
58
59 numValues++;
60 }
61};
62
63
64class FitsDumper : public factfits
65{
66private:
67 string fFilename;
68
69 // Convert CCfits::ValueType into a human readable string
70 string ValueTypeToStr(char type) const;
71
72 /// Lists all columns of an open file
73 void List();
74 void ListHeader(const string& filename);
75 void ListKeywords(ostream &);
76
77 vector<MyColumn> InitColumns(vector<string> list);
78 vector<MyColumn> InitColumnsRoot(vector<string> &list);
79
80 double GetDouble(const MyColumn &, size_t) const;
81 int64_t GetInteger(const MyColumn &, size_t) const;
82 string Format(const string &fmt, const double &val) const;
83 string Format(const string &fmt, const MyColumn &, size_t) const;
84
85 ///Display the selected columns values VS time
86 void Dump(ofstream &, const vector<string> &, const vector<MyColumn> &, const string &, size_t, size_t, const string &);
87 void DumpRoot(ofstream &, const vector<string> &, const string &, size_t, size_t, const string &);
88 void DumpMinMax(ofstream &, const vector<MyColumn> &, size_t, size_t, bool);
89 void DumpStats(ofstream &, const vector<MyColumn> &, const string &, size_t, size_t);
90
91public:
92 FitsDumper(const string &fname, const string &tablename);
93
94 ///Configures the fitsLoader from the config file and/or command arguments.
95 int Exec(Configuration& conf);
96};
97
98// --------------------------------------------------------------------------
99//
100//! Constructor
101//! @param out
102//! the ostream where to redirect the outputs
103//
104FitsDumper::FitsDumper(const string &fname, const string &tablename) : factfits(fname, tablename), fFilename(fname)
105{
106}
107
108string FitsDumper::ValueTypeToStr(char type) const
109{
110 switch (type)
111 {
112 case 'L': return "bool(8)";
113 case 'A': return "char(8)";
114 case 'B': return "byte(8)";
115 case 'I': return "short(16)";
116 case 'J': return "int(32)";
117 case 'K': return "int(64)";
118 case 'E': return "float(32)";
119 case 'D': return "double(64)";
120 default:
121 return "unknown";
122 }
123}
124
125void FitsDumper::List()
126{
127 const fits::Table::Keys &fKeyMap = GetKeys();
128 const fits::Table::Columns &fColMap = GetColumns();
129
130 cout << "\nFile: " << fFilename << "\n";
131
132 cout << " " << fKeyMap.find("EXTNAME")->second.value << " [";
133 cout << fKeyMap.find("NAXIS2")->second.value << "]\n";
134
135 for (auto it = fColMap.begin(); it != fColMap.end(); it++)
136 {
137 cout << " " << it->first << "[" << it->second.num << "] (" << it->second.unit << ":" << ValueTypeToStr(it->second.type) << ") ";
138 for (auto jt = fKeyMap.begin(); jt != fKeyMap.end(); jt++)
139 if (jt->second.value == it->first)
140 cout << "/ " << jt->second.comment << endl;
141 }
142
143 cout << endl;
144}
145
146void FitsDumper::ListKeywords(ostream &fout)
147{
148 const fits::Table::Keys &fKeyMap = GetKeys();
149
150 for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++)
151 {
152 fout << "## " << ::left << setw(8) << it->first << "= ";
153
154 if (it->second.type=='T')
155 fout << ::left << setw(20) << ("'"+it->second.value+"'");
156 else
157 fout << ::right << setw(20) << it->second.value;
158
159 if (!it->second.comment.empty())
160 fout << " / " << it->second.comment;
161 fout << '\n';
162 }
163
164 fout << flush;
165}
166
167void FitsDumper::ListHeader(const string& filename)
168{
169 ofstream fout(filename=="-"?"/dev/stdout":filename);
170 if (!fout)
171 {
172 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
173 return;
174 }
175
176 const fits::Table::Keys &fKeyMap = GetKeys();
177
178 fout << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n";
179 if (fKeyMap.find("COMMENT") != fKeyMap.end())
180 fout << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n";
181
182 ListKeywords(fout);
183 fout << endl;
184}
185
186vector<MyColumn> FitsDumper::InitColumns(vector<string> names)
187{
188 static const boost::regex expr("([[:word:].]+)(\\[([[:digit:]]+)?(:)?([[:digit:]]+)?\\])?");
189
190 const fits::Table::Columns &fColMap = GetColumns();
191
192 if (names.empty())
193 for (auto it=fColMap.begin(); it!=fColMap.end(); it++)
194 if (it->second.num>0)
195 names.push_back(it->first);
196
197 vector<MyColumn> vec;
198
199 for (auto it=names.begin(); it!=names.end(); it++)
200 {
201 boost::smatch what;
202 if (!boost::regex_match(*it, what, expr, boost::match_extra))
203 {
204 cerr << "Couldn't parse expression '" << *it << "' " << endl;
205 return vector<MyColumn>();
206 }
207
208 const string name = what[1];
209
210 const auto iter = fColMap.find(name);
211 if (iter==fColMap.end())
212 {
213 cerr << "ERROR - Column '" << name << "' not found in table." << endl;
214 return vector<MyColumn>();
215 }
216
217 const fits::Table::Column &col = iter->second;
218
219 const string val0 = what[3];
220 const string delim = what[4];
221 const string val1 = what[5];
222
223 const uint32_t first = atol(val0.c_str());
224 const uint32_t last = (val0.empty() && delim.empty()) ? col.num-1 : (val1.empty() ? first : atoi(val1.c_str()));
225
226 if (first>=col.num)
227 {
228 cerr << "ERROR - First index " << first << " for column " << name << " exceeds number of elements " << col.num << endl;
229 return vector<MyColumn>();
230 }
231
232 if (last>=col.num)
233 {
234 cerr << "ERROR - Last index " << last << " for column " << name << " exceeds number of elements " << col.num << endl;
235 return vector<MyColumn>();
236 }
237
238 if (first>last)
239 {
240 cerr << "ERROR - Last index " << last << " for column " << name << " exceeds first index " << first << endl;
241 return vector<MyColumn>();
242 }
243
244 MyColumn mycol;
245
246 mycol.name = name;
247 mycol.col = col;
248 mycol.first = first;
249 mycol.last = last;
250 mycol.ptr = SetPtrAddress(name);
251
252 vec.push_back(mycol);
253 }
254
255 return vec;
256}
257
258double FitsDumper::GetDouble(const MyColumn &it, size_t i) const
259{
260 switch (it.col.type)
261 {
262 case 'A':
263 return reinterpret_cast<const char*>(it.ptr)[i];
264
265 case 'L':
266 return reinterpret_cast<const bool*>(it.ptr)[i];
267
268 case 'B':
269 return (unsigned int)reinterpret_cast<const uint8_t*>(it.ptr)[i];
270
271 case 'I':
272 return reinterpret_cast<const int16_t*>(it.ptr)[i];
273
274 case 'J':
275 return reinterpret_cast<const int32_t*>(it.ptr)[i];
276
277 case 'K':
278 return reinterpret_cast<const int64_t*>(it.ptr)[i];
279
280 case 'E':
281 return reinterpret_cast<const float*>(it.ptr)[i];
282
283 case 'D':
284 return reinterpret_cast<const double*>(it.ptr)[i];
285 }
286
287 return 0;
288}
289
290int64_t FitsDumper::GetInteger(const MyColumn &it, size_t i) const
291{
292 switch (it.col.type)
293 {
294 case 'A':
295 return reinterpret_cast<const char*>(it.ptr)[i];
296
297 case 'L':
298 return reinterpret_cast<const bool*>(it.ptr)[i];
299
300 case 'B':
301 return (unsigned int)reinterpret_cast<const uint8_t*>(it.ptr)[i];
302
303 case 'I':
304 return reinterpret_cast<const int16_t*>(it.ptr)[i];
305
306 case 'J':
307 return reinterpret_cast<const int32_t*>(it.ptr)[i];
308
309 case 'K':
310 return reinterpret_cast<const int64_t*>(it.ptr)[i];
311
312 case 'E':
313 return reinterpret_cast<const float*>(it.ptr)[i];
314
315 case 'D':
316 return reinterpret_cast<const double*>(it.ptr)[i];
317 }
318
319 return 0;
320}
321
322string FitsDumper::Format(const string &format, const MyColumn &col, size_t i) const
323{
324 switch (*format.rbegin())
325 {
326 case 'd':
327 case 'i':
328 case 'o':
329 case 'u':
330 case 'x':
331 case 'X':
332 return Tools::Form(format.c_str(), GetDouble(col, i));
333
334 case 'e':
335 case 'E':
336 case 'f':
337 case 'F':
338 case 'g':
339 case 'G':
340 case 'a':
341 case 'A':
342 return Tools::Form(format.c_str(), GetInteger(col, i));
343
344 case 'h':
345 {
346 string rc = Tools::Scientific(GetDouble(col, i));
347 *remove_if(rc.begin(), rc.end(), ::isspace)=0;
348 return rc;
349 }
350 }
351
352 return "";
353}
354
355string FitsDumper::Format(const string &format, const double &val) const
356{
357 switch (*format.rbegin())
358 {
359 case 'd':
360 case 'i':
361 case 'o':
362 case 'u':
363 case 'x':
364 case 'X':
365 return Tools::Form(format.c_str(), int64_t(val));
366
367 case 'e':
368 case 'E':
369 case 'f':
370 case 'F':
371 case 'g':
372 case 'G':
373 case 'a':
374 case 'A':
375 return Tools::Form(format.c_str(), val);
376
377 case 'h':
378 {
379 string rc = Tools::Scientific(val);
380 *remove_if(rc.begin(), rc.end(), ::isspace)=0;
381 return rc;
382 }
383 }
384
385 return "";
386}
387
388
389// --------------------------------------------------------------------------
390//
391//! Perform the actual dump, based on the current parameters
392//
393void FitsDumper::Dump(ofstream &fout, const vector<string> &format, const vector<MyColumn> &cols, const string &filter, size_t first, size_t limit, const string &filename)
394{
395 const fits::Table::Keys &fKeyMap = GetKeys();
396
397#ifdef HAVE_ROOT
398 TFormula select;
399 if (!filter.empty() && select.Compile(filter.c_str()))
400 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+filter+"'");
401#endif
402
403 fout << "## --------------------------------------------------------------------------\n";
404 fout << "## Fits file: \t" << fFilename << '\n';
405 if (filename!="-")
406 fout << "## File: \t" << filename << '\n';
407 fout << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
408 fout << "## NumRows: \t" << GetInt("NAXIS2") << '\n';
409 fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
410#ifdef HAVE_ROOT
411 if (!filter.empty())
412 fout << "## Selection: \t" << select.GetExpFormula() << '\n';
413#endif
414 fout << "## --------------------------------------------------------------------------\n";
415 ListKeywords(fout);
416 fout << "## --------------------------------------------------------------------------\n";
417 fout << "#\n";
418
419 size_t num = 0;
420 for (auto it=cols.begin(); it!=cols.end(); it++)
421 {
422 fout << "# " << it->name;
423
424 if (it->first==it->last)
425 {
426 if (it->first!=0)
427 fout << "[" << it->first << "]";
428 }
429 else
430 fout << "[" << it->first << ":" << it->last << "]";
431
432 if (!it->col.unit.empty())
433 fout << ": " << it->col.unit;
434 fout << '\n';
435
436 num += it->last-it->first+1;
437 }
438 fout << "#" << endl;
439
440 // -----------------------------------------------------------------
441
442#ifdef HAVE_ROOT
443 vector<Double_t> data(num+1);
444#endif
445
446 const size_t last = limit ? first + limit : size_t(-1);
447
448 while (GetRow(first++))
449 {
450 const size_t row = GetRow();
451 if (row==GetNumRows() || row==last)
452 break;
453
454 size_t p = 0;
455
456#ifdef HAVE_ROOT
457 data[p++] = first-1;
458#endif
459
460 ostringstream sout;
461 sout.precision(fout.precision());
462 sout.flags(fout.flags());
463
464 uint32_t col = 0;
465 for (auto it=cols.begin(); it!=cols.end(); it++, col++)
466 {
467 string msg;
468 for (uint32_t i=it->first; i<=it->last; i++, p++)
469 {
470 if (col<format.size())
471 sout << Format("%"+format[col], *it, i) << " ";
472 else
473 {
474 switch (it->col.type)
475 {
476 case 'A':
477 msg += reinterpret_cast<const char*>(it->ptr)[i];
478 break;
479 case 'B':
480 sout << (unsigned int)reinterpret_cast<const unsigned char*>(it->ptr)[i] << " ";
481 break;
482 case 'L':
483 sout << reinterpret_cast<const bool*>(it->ptr)[i] << " ";
484 break;
485 case 'I':
486 sout << reinterpret_cast<const int16_t*>(it->ptr)[i] << " ";
487 break;
488 case 'J':
489 sout << reinterpret_cast<const int32_t*>(it->ptr)[i] << " ";
490 break;
491 case 'K':
492 sout << reinterpret_cast<const int64_t*>(it->ptr)[i] << " ";
493 break;
494 case 'E':
495 sout << reinterpret_cast<const float*>(it->ptr)[i] << " ";
496 break;
497 case 'D':
498 sout << reinterpret_cast<const double*>(it->ptr)[i] << " ";
499 break;
500 default:
501 ;
502 }
503 }
504#ifdef HAVE_ROOT
505 if (!filter.empty())
506 data[p] = GetDouble(*it, i);
507#endif
508 }
509
510 if (it->col.type=='A')
511 sout << "'" << msg.c_str() << "' ";
512 }
513#ifdef HAVE_ROOT
514 if (!filter.empty() && select.EvalPar(0, data.data())<0.5)
515 continue;
516#endif
517 fout << sout.str() << endl;
518 }
519}
520
521vector<MyColumn> FitsDumper::InitColumnsRoot(vector<string> &names)
522{
523 static const boost::regex expr("[^\\[]([[:word:].]+)(\\[([[:digit:]]+)\\])?");
524
525 const fits::Table::Columns &cols = GetColumns();
526
527 vector<MyColumn> vec;
528
529 for (auto it=names.begin(); it!=names.end(); it++)
530 {
531 if (it->empty())
532 continue;
533
534 *it = ' '+*it;
535
536 string::const_iterator ibeg = it->begin();
537 string::const_iterator iend = it->end();
538
539 boost::smatch what;
540 while (boost::regex_search(ibeg, iend, what, expr, boost::match_extra))
541 {
542 const string all = what[0];
543 const string name = what[1];
544 const size_t idx = atol(string(what[3]).c_str());
545
546 // Check if found colum is valid
547 const auto ic = cols.find(name);
548 if (ic==cols.end())
549 {
550 ibeg++;
551 //cout << "Column '" << name << "' does not exist." << endl;
552 //return vector<MyColumn>();
553 continue;
554 }
555 if (idx>=ic->second.num)
556 {
557 cout << "Column '" << name << "' has no index " << idx << "." << endl;
558 return vector<MyColumn>();
559 }
560
561 // find index if column already exists
562 size_t p = 0;
563 for (; p<vec.size(); p++)
564 if (vec[p].name==name)
565 break;
566
567 ostringstream id;
568 id << '[' << p << ']';
569
570 it->replace(ibeg-it->begin()+what.position(1), what.length()-1, id.str());
571
572 ibeg = what[0].first+3;
573 iend = it->end();
574
575 if (p<vec.size())
576 continue;
577
578 // Column not found, add new column
579 MyColumn mycol;
580
581 mycol.name = name;
582 mycol.col = ic->second;
583 mycol.first = idx;
584 mycol.last = idx;
585 mycol.ptr = SetPtrAddress(name);
586
587 vec.push_back(mycol);
588 }
589 }
590
591 ostringstream id;
592 id << '[' << vec.size() << ']';
593
594 for (auto it=names.begin(); it!=names.end(); it++)
595 {
596 while (1)
597 {
598 auto p = it->find_first_of('#');
599 if (p==string::npos)
600 break;
601
602 it->replace(p, 1, id.str());
603 }
604 }
605
606 //cout << endl;
607 //for (size_t i=0; i<vec.size(); i++)
608 // cout << "val[" << i << "] = " << vec[i].name << '[' << vec[i].first << ']' << endl;
609 //cout << endl;
610
611 return vec;
612}
613
614void FitsDumper::DumpRoot(ofstream &fout, const vector<string> &cols, const string &filter, size_t first, size_t limit, const string &filename)
615{
616#ifdef HAVE_ROOT
617 vector<string> names(cols);
618 names.insert(names.begin(), filter);
619
620 const vector<MyColumn> vec = InitColumnsRoot(names);
621 if (vec.empty())
622 return;
623
624 vector<TFormula> form(names.size());
625
626 auto ifo = form.begin();
627 for (auto it=names.begin(); it!=names.end(); it++, ifo++)
628 {
629 if (!it->empty() && ifo->Compile(it->c_str()))
630 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+*it+"'");
631 }
632
633 const fits::Table::Keys &fKeyMap = GetKeys();
634
635 fout << "## --------------------------------------------------------------------------\n";
636 fout << "## Fits file: \t" << fFilename << '\n';
637 if (filename!="-")
638 fout << "## File: \t" << filename << '\n';
639 fout << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
640 fout << "## NumRows: \t" << GetInt("NAXIS2") << '\n';
641 fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
642 fout << "## --------------------------------------------------------------------------\n";
643 ListKeywords(fout);
644 fout << "## --------------------------------------------------------------------------\n";
645 fout << "##\n";
646 if (!filter.empty())
647 fout << "## Selection: " << form[0].GetExpFormula() << "\n##\n";
648
649 size_t num = 0;
650 for (auto it=vec.begin(); it!=vec.end(); it++, num++)
651 {
652 fout << "## [" << num << "] = " << it->name;
653
654 if (it->first==it->last)
655 {
656 if (it->first!=0)
657 fout << "[" << it->first << "]";
658 }
659 else
660 fout << "[" << it->first << ":" << it->last << "]";
661
662 if (!it->col.unit.empty())
663 fout << ": " << it->col.unit;
664 fout << '\n';
665 }
666 fout << "##\n";
667 fout << "## --------------------------------------------------------------------------\n";
668 fout << "#\n";
669
670 fout << "# ";
671 for (auto it=form.begin()+1; it!=form.end(); it++)
672 fout << " \"" << it->GetExpFormula() << "\"";
673 fout << "\n#" << endl;
674
675 // -----------------------------------------------------------------
676
677 vector<Double_t> data(vec.size()+1);
678
679 const size_t last = limit ? first + limit : size_t(-1);
680
681 while (GetRow(first++))
682 {
683 const size_t row = GetRow();
684 if (row==GetNumRows() || row==last)
685 break;
686
687 size_t p = 0;
688 for (auto it=vec.begin(); it!=vec.end(); it++, p++)
689 data[p] = GetDouble(*it, it->first);
690
691 data[p] = first;
692
693 if (!filter.empty() && form[0].EvalPar(0, data.data())<0.5)
694 continue;
695
696 for (auto iform=form.begin()+1; iform!=form.end(); iform++)
697 fout << iform->EvalPar(0, data.data()) << " ";
698
699 fout << endl;
700 }
701#endif
702}
703
704void FitsDumper::DumpMinMax(ofstream &fout, const vector<MyColumn> &cols, size_t first, size_t limit, bool fNoZeroPlease)
705{
706 vector<minMaxStruct> statData(cols.size());
707
708 // Loop over all columns in our list of requested columns
709 const size_t last = limit ? first + limit : size_t(-1);
710
711 while (GetRow(first++))
712 {
713 const size_t row = GetRow();
714 if (row==GetNumRows() || row==last)
715 break;
716
717 auto statsIt = statData.begin();
718
719 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
720 {
721 if ((it->name=="UnixTimeUTC" || it->name=="PCTime") && it->first==0 && it->last==1)
722 {
723 const uint32_t *val = reinterpret_cast<const uint32_t*>(it->ptr);
724 if (fNoZeroPlease && val[0]==0 && val[1]==0)
725 continue;
726
727 statsIt->add(Time(val[0], val[1]).Mjd());
728 continue;
729 }
730
731 for (uint32_t i=it->first; i<=it->last; i++)
732 {
733 const double cValue = GetDouble(*it, i);
734
735 if (fNoZeroPlease && cValue == 0)
736 continue;
737
738 statsIt->add(cValue);
739 }
740 }
741 }
742
743 // okay. So now I've got ALL the data, loaded.
744 // let's do the summing and averaging in a safe way (i.e. avoid overflow
745 // of variables as much as possible)
746 auto statsIt = statData.begin();
747 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
748 {
749 fout << "\n[" << it->name << ':' << it->first;
750 if (it->first!=it->last)
751 fout << ':' << it->last;
752 fout << "]\n";
753
754 if (statsIt->numValues == 0)
755 {
756 fout << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
757 continue;
758 }
759
760 const long &num = statsIt->numValues;
761
762 long double &avg = statsIt->average;
763 long double &rms = statsIt->squared;
764
765 avg /= num;
766 rms = sqrt(rms/num - avg*avg);
767
768 fout << "Min: " << statsIt->min << '\n';
769 fout << "Max: " << statsIt->max << '\n';
770 fout << "Avg: " << avg << '\n';
771 fout << "Rms: " << rms << endl;
772 }
773}
774
775template<typename T>
776void displayStats(vector<char> &array, ofstream& out)
777{
778 const size_t numElems = array.size()/sizeof(T);
779 if (numElems == 0)
780 {
781 out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
782 return;
783 }
784
785 T *val = reinterpret_cast<T*>(array.data());
786
787 sort(val, val+numElems);
788
789 out << "Min: " << double(val[0]) << '\n';
790 out << "Max: " << double(val[numElems-1]) << '\n';
791
792 if (numElems%2 == 0)
793 out << "Med: " << (double(val[numElems/2-1]) + double(val[numElems/2]))/2 << '\n';
794 else
795 out << "Med: " << double(val[numElems/2]) << '\n';
796
797 long double avg = 0;
798 long double rms = 0;
799 for (uint32_t i=0;i<numElems;i++)
800 {
801 const long double v = val[i];
802 avg += v;
803 rms += v*v;
804 }
805
806 avg /= numElems;
807 rms = sqrt(rms/numElems - avg*avg);
808
809 out << "Avg: " << avg << '\n';
810 out << "Rms: " << rms << endl;
811
812}
813
814void FitsDumper::DumpStats(ofstream &fout, const vector<MyColumn> &cols, const string &filter, size_t first, size_t limit)
815{
816#ifdef HAVE_ROOT
817 TFormula select;
818 if (!filter.empty() && select.Compile(filter.c_str()))
819 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+filter+"'");
820#endif
821
822 // Loop over all columns in our list of requested columns
823 vector<vector<char>> statData;
824
825 const size_t num = limit==0 || GetNumRows()<limit ? GetNumRows() : limit;
826
827 for (auto it=cols.begin(); it!=cols.end(); it++)
828 statData.emplace_back(vector<char>(it->col.size*num*(it->last-it->first+1)));
829
830#ifdef HAVE_ROOT
831 vector<Double_t> data(num+1);
832#endif
833
834 // Loop over all columns in our list of requested columns
835 const size_t last = limit ? first + limit : size_t(-1);
836
837 uint64_t counter = 0;
838
839 while (GetRow(first++))
840 {
841 const size_t row = GetRow();
842 if (row==GetNumRows() || row==last)
843 break;
844
845#ifdef HAVE_ROOT
846 if (!filter.empty())
847 {
848 size_t p = 0;
849
850 data[p++] = first-1;
851
852 for (auto it=cols.begin(); it!=cols.end(); it++)
853 for (uint32_t i=it->first; i<=it->last; i++, p++)
854 data[p] = GetDouble(*it, i);
855
856 if (select.EvalPar(0, data.data())<0.5)
857 continue;
858 }
859#endif
860
861 auto statsIt = statData.begin();
862 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
863 {
864 const char *src = reinterpret_cast<const char*>(it->ptr);
865 const size_t sz = (it->last-it->first+1)*it->col.size;
866 memcpy(statsIt->data()+counter*sz, src+it->first*it->col.size, sz);
867 }
868
869 counter++;
870 }
871
872 auto statsIt = statData.begin();
873 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
874 {
875 fout << "\n[" << it->name << ':' << it->first;
876 if (it->last!=it->first)
877 fout << ':' << it->last;
878 fout << "]\n";
879
880 const size_t sz = (it->last-it->first+1)*it->col.size;
881 statsIt->resize(counter*sz);
882
883 switch (it->col.type)
884 {
885 case 'L':
886 displayStats<bool>(*statsIt, fout);
887 break;
888 case 'B':
889 displayStats<char>(*statsIt, fout);
890 break;
891 case 'I':
892 displayStats<int16_t>(*statsIt, fout);
893 break;
894 case 'J':
895 displayStats<int32_t>(*statsIt, fout);
896 break;
897 case 'K':
898 displayStats<int64_t>(*statsIt, fout);
899 break;
900 case 'E':
901 displayStats<float>(*statsIt, fout);
902 break;
903 case 'D':
904 displayStats<double>(*statsIt, fout);
905 break;
906 default:
907 ;
908 }
909 }
910}
911
912// --------------------------------------------------------------------------
913//
914//! Retrieves the configuration parameters
915//! @param conf
916//! the configuration object
917//
918int FitsDumper::Exec(Configuration& conf)
919{
920 if (conf.Get<bool>("list"))
921 List();
922
923 if (conf.Get<bool>("header"))
924 ListHeader(conf.Get<string>("outfile"));
925
926
927 if (conf.Get<bool>("header") || conf.Get<bool>("list"))
928 return 1;
929
930 // ------------------------------------------------------------
931
932 if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
933 {
934 cerr << "Invalid combination of options: cannot do stats and minmax." << endl;
935 return -1;
936 }
937 if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
938 {
939 cerr << "Invalid combination of options: nozero only works with minmax." << endl;
940 return -1;
941 }
942
943 if (conf.Get<bool>("scientific") && conf.Get<bool>("fixed"))
944 {
945 cerr << "Switched --scientific and --fixed are mutually exclusive." << endl;
946 return -1;
947 }
948
949 if (conf.Has("%") && conf.Has("%%"))
950 {
951 cerr << "Switched --% and --%% are mutually exclusive." << endl;
952 return -1;
953 }
954
955 // ------------------------------------------------------------
956
957 const string filename = conf.Get<string>("outfile");
958
959 ofstream fout(filename=="-"?"/dev/stdout":filename);
960 if (!fout)
961 {
962 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
963 return false;
964 }
965 fout.precision(conf.Get<int>("precision"));
966 if (conf.Get<bool>("fixed"))
967 fout << fixed;
968 if (conf.Get<bool>("scientific"))
969 fout << scientific;
970
971 const string filter = conf.Has("filter") ? conf.Get<string>("filter") : "";
972 const size_t first = conf.Get<size_t>("first");
973 const size_t limit = conf.Get<size_t>("limit");
974
975#ifdef HAVE_ROOT
976 if (conf.Get<bool>("root"))
977 {
978 DumpRoot(fout, conf.Vec<string>("col"), filter, first, limit, filename);
979 return 0;
980 }
981#endif
982
983 const vector<string> format = conf.Vec<string>("%");
984 for (auto it=format.begin(); it<format.end(); it++)
985 {
986 static const boost::regex expr("-?[0-9]*[.]?[0-9]*[diouxXeEfFgGaAh]");
987
988 boost::smatch what;
989 if (!boost::regex_match(*it, what, expr, boost::match_extra))
990 {
991 cerr << "Format '" << *it << "' not supported." << endl;
992 return -1;
993 }
994 }
995
996 const vector<MyColumn> cols = InitColumns(conf.Vec<string>("col"));
997 if (cols.empty())
998 return false;
999
1000 if (conf.Get<bool>("minmax"))
1001 {
1002 DumpMinMax(fout, cols, first, limit, conf.Get<bool>("nozero"));
1003 return 0;
1004 }
1005
1006 if (conf.Get<bool>("stat"))
1007 {
1008 DumpStats(fout, cols, filter, first, limit);
1009 return 0;
1010 }
1011
1012 Dump(fout, format, cols, filter, first, limit, filename);
1013
1014 return 0;
1015}
1016
1017void PrintUsage()
1018{
1019 cout <<
1020 "fitsdump is a tool to dump data from a FITS table as ascii.\n"
1021 "\n"
1022 "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
1023 " or: fitsdump [OPTIONS]\n"
1024 "\n"
1025 "Addressing a column:\n"
1026 " ColumnName: Will address all fields of a column\n"
1027 " ColumnName[n]: Will address the n-th field of a column (starts with 0)\n"
1028 " ColumnName[n1:n2]: Will address all fields between n1 and including n2\n"
1029#ifdef HAVE_ROOT
1030 "\n"
1031 "Selecting a column:\n"
1032 " Commandline option: --filter\n"
1033 " Explanation: Such a selection is evaluated using TFormula, hence, every "
1034 "mathematical operation allowed in TFormula is allowed there, too. "
1035 "The reference is the column index as printed in the output stream, "
1036 "starting with 1. The index 0 is reserved for the row number.\n"
1037#endif
1038 ;
1039 cout << endl;
1040}
1041
1042void PrintHelp()
1043{
1044#ifdef HAVE_ROOT
1045 cout <<
1046 "\n\n"
1047 "Examples:\n"
1048 "In --root mode, fitsdump support TFormula's syntax for all columns and the filter "
1049 "You can then refer to a column or a (single) index of the column just by its name "
1050 "If the index is omitted, 0 is assumed. Note that the [x:y] syntax in this mode is "
1051 "not supported\n"
1052 "\n"
1053 " fitsdump Zd --filter=\"[0]>20 && cos([1])*TMath::RadToDeg()<45\"\n"
1054 "\n"
1055 "The columns can also be addressed with their names\n"
1056 "\n"
1057 " fitsdump -r \"(Zd+Err)*TMath::DegToRad()\" --filter=\"[0]<25 && [1]<0.05\"\n"
1058 "\n"
1059 "is identical to\n"
1060 "\n"
1061 " fitsdump -r \"(Zd[0]+Err[0])*TMath::DegToRad()\" --filter=\"[0]<25 && [1]<0.05\"\n"
1062 "\n"
1063 "A special placeholder exists for the row number\n"
1064 "\n"
1065 " fitsdump -r \"#\" --filter=\"#>10 && #<100\"\n"
1066 "\n"
1067 "To format a single column you can do\n"
1068 "\n"
1069 " fitsdump col1 -%.1f col2 -%d\n"
1070 "\n"
1071 "A special format is provided converting to 'human readable format'\n"
1072 "\n"
1073 " fitsdump col1 -%h\n"
1074 "\n";
1075 cout << endl;
1076#endif
1077}
1078
1079
1080void SetupConfiguration(Configuration& conf)
1081{
1082 po::options_description configs("Fitsdump options");
1083 configs.add_options()
1084 ("fitsfile,f", var<string>()
1085#if BOOST_VERSION >= 104200
1086 ->required()
1087#endif
1088 , "Name of FITS file")
1089 ("col,c", vars<string>(), "List of columns to dump\narg is a list of columns, separated by a space.\nAdditionnally, a list of sub-columns can be added\ne.g. Data[3] will dump sub-column 3 of column Data\nData[3:4] will dump sub-columns 3 and 4\nOmitting this argument dump the entire column\nnota: all indices start at zero")
1090 ("outfile,o", var<string>("-"), "Name of output file (-:/dev/stdout)")
1091 ("precision,p", var<int>(20), "Precision of ofstream")
1092 ("list,l", po_switch(), "List all tables and columns in file")
1093 ("header,h", po_switch(), "Dump header of given table")
1094 ("stat,s", po_switch(), "Perform statistics instead of dump")
1095 ("minmax,m", po_switch(), "Calculates min and max of data")
1096 ("nozero,z", po_switch(), "skip 0 values for stats")
1097 ("fixed", po_switch(), "Switch output stream to floating point values in fixed-point notation")
1098 ("scientific", po_switch(), "Switch output stream to floating point values in scientific notation")
1099 ("%,%", vars<string>(), "Format for the output (currently not available in root-mode)")
1100 ("force", po_switch(), "Force reading the fits file even if END key is missing")
1101 ("first", var<size_t>(0), "First number of row to read")
1102 ("limit", var<size_t>(0), "Limit for the maximum number of rows to read (0=unlimited)")
1103 ("tablename,t", var<string>(""), "Name of the table to open. If not specified, first binary table is opened")
1104#ifdef HAVE_ROOT
1105 ("root,r", po_switch(), "Enable root mode")
1106 ("filter,f", var<string>(""), "Filter to restrict the selection of events (e.g. '[0]>10 && [0]<20'; does not work with stat and minmax yet)")
1107#endif
1108 ;
1109
1110 po::positional_options_description p;
1111 p.add("fitsfile", 1); // The first positional options
1112 p.add("col", -1); // All others
1113
1114 conf.AddOptions(configs);
1115 conf.SetArgumentPositions(p);
1116}
1117
1118int main(int argc, const char** argv)
1119{
1120 Configuration conf(argv[0]);
1121 conf.SetPrintUsage(PrintUsage);
1122 SetupConfiguration(conf);
1123
1124 if (!conf.DoParse(argc, argv, PrintHelp))
1125 return -1;
1126
1127 if (!conf.Has("fitsfile"))
1128 {
1129 cerr << "Filename required." << endl;
1130 return -1;
1131 }
1132
1133 FitsDumper loader(conf.Get<string>("fitsfile"), conf.Get<string>("tablename"));
1134 if (!loader)
1135 {
1136 cerr << "ERROR - Opening " << conf.Get<string>("fitsfile");
1137 cerr << " failed: " << strerror(errno) << endl;
1138 return -1;
1139 }
1140
1141 return loader.Exec(conf);
1142}
Note: See TracBrowser for help on using the repository browser.