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

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