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

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