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

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