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

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