source: branches/testFACT++branch/src/fitsdump.cc@ 18302

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