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

Last change on this file since 18826 was 18602, checked in by tbretz, 8 years ago
The latest change broke the compiling with root.
File size: 34.6 KB
Line 
1//****************************************************************
2/** @class FitsDumper
3
4 @brief Dumps contents of fits tables to stdout or a file
5
6 */
7 //****************************************************************
8#include "Configuration.h"
9
10#include <float.h>
11
12#include <map>
13#include <fstream>
14
15#include <boost/regex.hpp>
16
17#include "tools.h"
18#include "Time.h"
19#include "externals/factfits.h"
20
21#ifdef HAVE_ROOT
22#include "TFormula.h"
23#endif
24
25using namespace std;
26
27struct MyColumn
28{
29 string name;
30
31 fits::Table::Column col;
32
33 uint32_t first;
34 uint32_t last;
35
36 void *ptr;
37};
38
39struct minMaxStruct
40{
41 double min;
42 double max;
43 long double average;
44 long double squared;
45 long numValues;
46 minMaxStruct() : min(FLT_MAX), max(-FLT_MAX), average(0), squared(0), numValues(0) { }
47
48 void add(long double val)
49 {
50 average += val;
51 squared += val*val;
52
53 if (val<min)
54 min = val;
55
56 if (val>max)
57 max = val;
58
59 numValues++;
60 }
61};
62
63
64class FitsDumper : public factfits
65{
66private:
67 string fFilename;
68
69 // Convert CCfits::ValueType into a human readable string
70 string ValueTypeToStr(char type) const;
71
72 /// Lists all columns of an open file
73 void List();
74 void ListFileContent() const;
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//
410#ifdef HAVE_ROOT
411void FitsDumper::Dump(ostream &fout, const vector<string> &format, const vector<MyColumn> &cols, const string &filter, size_t first, size_t limit, const string &filename)
412#else
413void FitsDumper::Dump(ostream &fout, const vector<string> &format, const vector<MyColumn> &cols, const string &, size_t first, size_t limit, const string &filename)
414#endif
415{
416 const fits::Table::Keys &fKeyMap = GetKeys();
417
418#ifdef HAVE_ROOT
419 TFormula select;
420 if (!filter.empty() && select.Compile(filter.c_str()))
421 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+filter+"'");
422#endif
423
424 fout << "## --------------------------------------------------------------------------\n";
425 fout << "## Fits file: \t" << fFilename << '\n';
426 if (filename!="-")
427 fout << "## File: \t" << filename << '\n';
428 fout << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
429 fout << "## NumRows: \t" << GetNumRows() << '\n';
430 fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
431#ifdef HAVE_ROOT
432 if (!filter.empty())
433 fout << "## Selection: \t" << select.GetExpFormula() << '\n';
434#endif
435 fout << "## --------------------------------------------------------------------------\n";
436 ListKeywords(fout);
437 fout << "## --------------------------------------------------------------------------\n";
438 fout << "#\n";
439
440 size_t num = 0;
441 for (auto it=cols.begin(); it!=cols.end(); it++)
442 {
443 fout << "# " << it->name;
444
445 if (it->first==it->last)
446 {
447 if (it->first!=0)
448 fout << "[" << it->first << "]";
449 }
450 else
451 fout << "[" << it->first << ":" << it->last << "]";
452
453 if (!it->col.unit.empty())
454 fout << ": " << it->col.unit;
455 fout << '\n';
456
457 num += it->last-it->first+1;
458 }
459 fout << "#" << endl;
460
461 // -----------------------------------------------------------------
462
463#ifdef HAVE_ROOT
464 vector<Double_t> data(num+1);
465#endif
466
467 const size_t last = limit ? first + limit : size_t(-1);
468
469 while (GetRow(first++))
470 {
471 const size_t row = GetRow();
472 if (row==GetNumRows() || row==last)
473 break;
474
475 size_t p = 0;
476
477#ifdef HAVE_ROOT
478 data[p++] = first-1;
479#endif
480
481 ostringstream sout;
482 sout.precision(fout.precision());
483 sout.flags(fout.flags());
484
485 uint32_t col = 0;
486 for (auto it=cols.begin(); it!=cols.end(); it++, col++)
487 {
488 string msg;
489 for (uint32_t i=it->first; i<=it->last; i++, p++)
490 {
491 if (col<format.size())
492 sout << Format("%"+format[col], *it, i) << " ";
493 else
494 {
495 switch (it->col.type)
496 {
497 case 'A':
498 msg += reinterpret_cast<const char*>(it->ptr)[i];
499 break;
500 case 'B':
501 sout << (unsigned int)reinterpret_cast<const unsigned char*>(it->ptr)[i] << " ";
502 break;
503 case 'L':
504 sout << reinterpret_cast<const bool*>(it->ptr)[i] << " ";
505 break;
506 case 'I':
507 sout << reinterpret_cast<const int16_t*>(it->ptr)[i] << " ";
508 break;
509 case 'J':
510 sout << reinterpret_cast<const int32_t*>(it->ptr)[i] << " ";
511 break;
512 case 'K':
513 sout << reinterpret_cast<const int64_t*>(it->ptr)[i] << " ";
514 break;
515 case 'E':
516 sout << reinterpret_cast<const float*>(it->ptr)[i] << " ";
517 break;
518 case 'D':
519 sout << reinterpret_cast<const double*>(it->ptr)[i] << " ";
520 break;
521 default:
522 ;
523 }
524 }
525#ifdef HAVE_ROOT
526 if (!filter.empty())
527 data[p] = GetDouble(*it, i);
528#endif
529 }
530
531 if (it->col.type=='A')
532 sout << "'" << msg.c_str() << "' ";
533 }
534#ifdef HAVE_ROOT
535 if (!filter.empty() && select.EvalPar(0, data.data())<0.5)
536 continue;
537#endif
538 fout << sout.str() << endl;
539 }
540}
541
542vector<MyColumn> FitsDumper::InitColumnsRoot(vector<string> &names)
543{
544 static const boost::regex expr("[^\\[]([[:word:].]+)(\\[([[:digit:]]+)\\])?");
545
546 const fits::Table::Columns &cols = GetColumns();
547
548 vector<MyColumn> vec;
549
550 for (auto it=names.begin(); it!=names.end(); it++)
551 {
552 if (it->empty())
553 continue;
554
555 *it = ' '+*it;
556
557 string::const_iterator ibeg = it->begin();
558 string::const_iterator iend = it->end();
559
560 boost::smatch what;
561 while (boost::regex_search(ibeg, iend, what, expr, boost::match_extra))
562 {
563 const string all = what[0];
564 const string name = what[1];
565 const size_t idx = atol(string(what[3]).c_str());
566
567 // Check if found colum is valid
568 const auto ic = cols.find(name);
569 if (ic==cols.end())
570 {
571 ibeg++;
572 //cout << "Column '" << name << "' does not exist." << endl;
573 //return vector<MyColumn>();
574 continue;
575 }
576 if (idx>=ic->second.num)
577 {
578 cout << "Column '" << name << "' has no index " << idx << "." << endl;
579 return vector<MyColumn>();
580 }
581
582 // find index if column already exists
583 size_t p = 0;
584 for (; p<vec.size(); p++)
585 if (vec[p].name==name)
586 break;
587
588 const string id = '['+to_string(p)+']';
589
590 // Replace might reallocate the memory. Therefore, we cannot use what[0].first
591 // directly but have to store the offset
592 const size_t offset = what[0].first - it->begin();
593
594 it->replace(ibeg-it->begin()+what.position(1), what.length()-1, id);
595
596 ibeg = it->begin() + offset + id.size();
597 iend = it->end();
598
599 if (p<vec.size())
600 continue;
601
602 // Column not found, add new column
603 MyColumn mycol;
604
605 mycol.name = name;
606 mycol.col = ic->second;
607 mycol.first = idx;
608 mycol.last = idx;
609 mycol.ptr = SetPtrAddress(name);
610
611 vec.push_back(mycol);
612 }
613 }
614
615 ostringstream id;
616 id << '[' << vec.size() << ']';
617
618 for (auto it=names.begin(); it!=names.end(); it++)
619 {
620 while (1)
621 {
622 auto p = it->find_first_of('#');
623 if (p==string::npos)
624 break;
625
626 it->replace(p, 1, id.str());
627 }
628 }
629
630 //cout << endl;
631 //for (size_t i=0; i<vec.size(); i++)
632 // cout << "val[" << i << "] = " << vec[i].name << '[' << vec[i].first << ']' << endl;
633 //cout << endl;
634
635 return vec;
636}
637
638#ifdef HAVE_ROOT
639void FitsDumper::DumpRoot(ostream &fout, const vector<string> &cols, const string &filter, size_t first, size_t limit, const string &filename)
640#else
641void FitsDumper::DumpRoot(ostream &, const vector<string> &, const string &, size_t, size_t, const string &)
642#endif
643{
644#ifdef HAVE_ROOT
645 vector<string> names(cols);
646 names.insert(names.begin(), filter);
647
648 const vector<MyColumn> vec = InitColumnsRoot(names);
649 if (vec.empty())
650 return;
651
652 vector<TFormula> form(names.size());
653
654 auto ifo = form.begin();
655 for (auto it=names.begin(); it!=names.end(); it++, ifo++)
656 {
657 if (!it->empty() && ifo->Compile(it->c_str()))
658 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+*it+"'");
659 }
660
661 const fits::Table::Keys &fKeyMap = GetKeys();
662
663 fout << "## --------------------------------------------------------------------------\n";
664 fout << "## Fits file: \t" << fFilename << '\n';
665 if (filename!="-")
666 fout << "## File: \t" << filename << '\n';
667 fout << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
668 fout << "## NumRows: \t" << GetNumRows() << '\n';
669 fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
670 fout << "## --------------------------------------------------------------------------\n";
671 ListKeywords(fout);
672 fout << "## --------------------------------------------------------------------------\n";
673 fout << "##\n";
674 if (!filter.empty())
675 fout << "## Selection: " << form[0].GetExpFormula() << "\n##\n";
676
677 size_t num = 0;
678 for (auto it=vec.begin(); it!=vec.end(); it++, num++)
679 {
680 fout << "## [" << num << "] = " << it->name;
681
682 if (it->first==it->last)
683 {
684 if (it->first!=0)
685 fout << "[" << it->first << "]";
686 }
687 else
688 fout << "[" << it->first << ":" << it->last << "]";
689
690 if (!it->col.unit.empty())
691 fout << ": " << it->col.unit;
692 fout << '\n';
693 }
694 fout << "##\n";
695 fout << "## --------------------------------------------------------------------------\n";
696 fout << "#\n";
697
698 fout << "# ";
699 for (auto it=form.begin()+1; it!=form.end(); it++)
700 fout << " \"" << it->GetExpFormula() << "\"";
701 fout << "\n#" << endl;
702
703 // -----------------------------------------------------------------
704
705 vector<Double_t> data(vec.size()+1);
706
707 const size_t last = limit ? first + limit : size_t(-1);
708
709 while (GetRow(first++))
710 {
711 const size_t row = GetRow();
712 if (row==GetNumRows() || row==last)
713 break;
714
715 size_t p = 0;
716 for (auto it=vec.begin(); it!=vec.end(); it++, p++)
717 data[p] = GetDouble(*it, it->first);
718
719 data[p] = first;
720
721 if (!filter.empty() && form[0].EvalPar(0, data.data())<0.5)
722 continue;
723
724 for (auto iform=form.begin()+1; iform!=form.end(); iform++)
725 fout << iform->EvalPar(0, data.data()) << " ";
726
727 fout << endl;
728 }
729#endif
730}
731
732void FitsDumper::DumpMinMax(ostream &fout, const vector<MyColumn> &cols, size_t first, size_t limit, bool fNoZeroPlease)
733{
734 vector<minMaxStruct> statData(cols.size());
735
736 // Loop over all columns in our list of requested columns
737 const size_t last = limit ? first + limit : size_t(-1);
738
739 while (GetRow(first++))
740 {
741 const size_t row = GetRow();
742 if (row==GetNumRows() || row==last)
743 break;
744
745 auto statsIt = statData.begin();
746
747 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
748 {
749 if ((it->name=="UnixTimeUTC" || it->name=="PCTime") && it->first==0 && it->last==1)
750 {
751 const uint32_t *val = reinterpret_cast<const uint32_t*>(it->ptr);
752 if (fNoZeroPlease && val[0]==0 && val[1]==0)
753 continue;
754
755 statsIt->add(Time(val[0], val[1]).Mjd());
756 continue;
757 }
758
759 for (uint32_t i=it->first; i<=it->last; i++)
760 {
761 const double cValue = GetDouble(*it, i);
762
763 if (fNoZeroPlease && cValue == 0)
764 continue;
765
766 statsIt->add(cValue);
767 }
768 }
769 }
770
771 // okay. So now I've got ALL the data, loaded.
772 // let's do the summing and averaging in a safe way (i.e. avoid overflow
773 // of variables as much as possible)
774 auto statsIt = statData.begin();
775 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
776 {
777 fout << "\n[" << it->name << ':' << it->first;
778 if (it->first!=it->last)
779 fout << ':' << it->last;
780 fout << "]\n";
781
782 if (statsIt->numValues == 0)
783 {
784 fout << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
785 continue;
786 }
787
788 const long &num = statsIt->numValues;
789
790 long double &avg = statsIt->average;
791 long double &rms = statsIt->squared;
792
793 avg /= num;
794 rms /= num;
795 rms += avg*avg;
796 rms = rms<0 ? 0 : sqrt(rms);
797
798 fout << "Min: " << statsIt->min << '\n';
799 fout << "Max: " << statsIt->max << '\n';
800 fout << "Avg: " << avg << '\n';
801 fout << "Rms: " << rms << endl;
802 }
803}
804
805template<typename T>
806void displayStats(vector<char> &array, ostream& out)
807{
808 const size_t numElems = array.size()/sizeof(T);
809 if (numElems == 0)
810 {
811 out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
812 return;
813 }
814
815 T *val = reinterpret_cast<T*>(array.data());
816
817 sort(val, val+numElems);
818
819 out << "Min: " << double(val[0]) << '\n';
820 out << "Max: " << double(val[numElems-1]) << '\n';
821
822 if (numElems%2 == 0)
823 out << "Med: " << (double(val[numElems/2-1]) + double(val[numElems/2]))/2 << '\n';
824 else
825 out << "Med: " << double(val[numElems/2]) << '\n';
826
827 long double avg = 0;
828 long double rms = 0;
829 for (uint32_t i=0;i<numElems;i++)
830 {
831 const long double v = val[i];
832 avg += v;
833 rms += v*v;
834 }
835
836 avg /= numElems;
837 rms /= numElems;
838 rms -= avg*avg;
839 rms = rms<0 ? 0 : sqrt(rms);
840
841
842 out << "Avg: " << avg << '\n';
843 out << "Rms: " << rms << endl;
844}
845
846#ifdef HAVE_ROOT
847void FitsDumper::DumpStats(ostream &fout, const vector<MyColumn> &cols, const string &filter, size_t first, size_t limit)
848#else
849void FitsDumper::DumpStats(ostream &fout, const vector<MyColumn> &cols, const string &, size_t first, size_t limit)
850#endif
851{
852#ifdef HAVE_ROOT
853 TFormula select;
854 if (!filter.empty() && select.Compile(filter.c_str()))
855 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+filter+"'");
856#endif
857
858 // Loop over all columns in our list of requested columns
859 vector<vector<char>> statData;
860
861 const size_t rows = limit==0 || GetNumRows()<limit ? GetNumRows() : limit;
862
863 for (auto it=cols.begin(); it!=cols.end(); it++)
864 statData.emplace_back(vector<char>(it->col.size*rows*(it->last-it->first+1)));
865
866#ifdef HAVE_ROOT
867 size_t num = 0;
868 for (auto it=cols.begin(); it!=cols.end(); it++)
869 num += it->last-it->first+1;
870
871 vector<Double_t> data(num+1);
872#endif
873
874 // Loop over all columns in our list of requested columns
875 const size_t last = limit ? first + limit : size_t(-1);
876
877 uint64_t counter = 0;
878
879 while (GetRow(first++))
880 {
881 const size_t row = GetRow();
882 if (row==GetNumRows() || row==last)
883 break;
884
885#ifdef HAVE_ROOT
886 if (!filter.empty())
887 {
888 size_t p = 0;
889
890 data[p++] = first-1;
891
892 for (auto it=cols.begin(); it!=cols.end(); it++)
893 for (uint32_t i=it->first; i<=it->last; i++, p++)
894 data[p] = GetDouble(*it, i);
895
896 if (select.EvalPar(0, data.data())<0.5)
897 continue;
898 }
899#endif
900
901 auto statsIt = statData.begin();
902 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
903 {
904 const char *src = reinterpret_cast<const char*>(it->ptr);
905 const size_t sz = (it->last-it->first+1)*it->col.size;
906 memcpy(statsIt->data()+counter*sz, src+it->first*it->col.size, sz);
907 }
908
909 counter++;
910 }
911
912 auto statsIt = statData.begin();
913 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
914 {
915 fout << "\n[" << it->name << ':' << it->first;
916 if (it->last!=it->first)
917 fout << ':' << it->last;
918 fout << "]\n";
919
920 const size_t sz = (it->last-it->first+1)*it->col.size;
921 statsIt->resize(counter*sz);
922
923 switch (it->col.type)
924 {
925 case 'L':
926 displayStats<bool>(*statsIt, fout);
927 break;
928 case 'B':
929 displayStats<char>(*statsIt, fout);
930 break;
931 case 'I':
932 displayStats<int16_t>(*statsIt, fout);
933 break;
934 case 'J':
935 displayStats<int32_t>(*statsIt, fout);
936 break;
937 case 'K':
938 displayStats<int64_t>(*statsIt, fout);
939 break;
940 case 'E':
941 displayStats<float>(*statsIt, fout);
942 break;
943 case 'D':
944 displayStats<double>(*statsIt, fout);
945 break;
946 default:
947 ;
948 }
949 }
950}
951
952// --------------------------------------------------------------------------
953//
954//! Retrieves the configuration parameters
955//! @param conf
956//! the configuration object
957//
958int FitsDumper::Exec(Configuration& conf)
959{
960 if (conf.Get<bool>("list"))
961 List();
962
963 if (conf.Get<bool>("filecontent"))
964 ListFileContent();
965
966 if (conf.Get<bool>("header"))
967 ListHeader(conf.Get<string>("outfile"));
968
969
970 if (conf.Get<bool>("header") || conf.Get<bool>("list") || conf.Get<bool>("filecontent"))
971 return 1;
972
973 // ------------------------------------------------------------
974
975 if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
976 {
977 cerr << "Invalid combination of options: cannot do stats and minmax." << endl;
978 return -1;
979 }
980 if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
981 {
982 cerr << "Invalid combination of options: nozero only works with minmax." << endl;
983 return -1;
984 }
985
986 if (conf.Get<bool>("scientific") && conf.Get<bool>("fixed"))
987 {
988 cerr << "Switched --scientific and --fixed are mutually exclusive." << endl;
989 return -1;
990 }
991
992 if (conf.Has("%") && conf.Has("%%"))
993 {
994 cerr << "Switched --% and --%% are mutually exclusive." << endl;
995 return -1;
996 }
997
998 // ------------------------------------------------------------
999
1000 const string filename = conf.Get<string>("outfile");
1001
1002 ostream fout(cout.rdbuf());
1003
1004 ofstream sout;
1005 if (filename!="-")
1006 {
1007 sout.open(filename);
1008 if (!sout)
1009 {
1010 cerr << "Cannot open output stream " << filename << ": " << strerror(errno) << endl;
1011 return false;
1012 }
1013 fout.rdbuf(sout.rdbuf());
1014 }
1015
1016 fout.precision(conf.Get<int>("precision"));
1017 if (conf.Get<bool>("fixed"))
1018 fout << fixed;
1019 if (conf.Get<bool>("scientific"))
1020 fout << scientific;
1021
1022 const string filter = conf.Has("filter") ? conf.Get<string>("filter") : "";
1023 const size_t first = conf.Get<size_t>("first");
1024 const size_t limit = conf.Get<size_t>("limit");
1025
1026#ifdef HAVE_ROOT
1027 if (conf.Get<bool>("root"))
1028 {
1029 DumpRoot(fout, conf.Vec<string>("col"), filter, first, limit, filename);
1030 return 0;
1031 }
1032#endif
1033
1034 const vector<string> format = conf.Vec<string>("%");
1035 for (auto it=format.begin(); it<format.end(); it++)
1036 {
1037 static const boost::regex expr("-?[0-9]*[.]?[0-9]*[diouxXeEfFgGaAh]");
1038
1039 boost::smatch what;
1040 if (!boost::regex_match(*it, what, expr, boost::match_extra))
1041 {
1042 cerr << "Format '" << *it << "' not supported." << endl;
1043 return -1;
1044 }
1045 }
1046
1047 const vector<MyColumn> cols = InitColumns(conf.Vec<string>("col"));
1048 if (cols.empty())
1049 return false;
1050
1051 if (conf.Get<bool>("minmax"))
1052 {
1053 DumpMinMax(fout, cols, first, limit, conf.Get<bool>("nozero"));
1054 return 0;
1055 }
1056
1057 if (conf.Get<bool>("stat"))
1058 {
1059 DumpStats(fout, cols, filter, first, limit);
1060 return 0;
1061 }
1062
1063 Dump(fout, format, cols, filter, first, limit, filename);
1064
1065 return 0;
1066}
1067
1068void PrintUsage()
1069{
1070 cout <<
1071 "fitsdump is a tool to dump data from a FITS table as ascii.\n"
1072 "\n"
1073 "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
1074 " or: fitsdump [OPTIONS]\n"
1075 "\n"
1076 "Addressing a column:\n"
1077 " ColumnName: Will address all fields of a column\n"
1078 " ColumnName[n]: Will address the n-th field of a column (starts with 0)\n"
1079 " ColumnName[n1:n2]: Will address all fields between n1 and including n2\n"
1080#ifdef HAVE_ROOT
1081 "\n"
1082 "Selecting a column:\n"
1083 " Commandline option: --filter\n"
1084 " Explanation: Such a selection is evaluated using TFormula, hence, every "
1085 "mathematical operation allowed in TFormula is allowed there, too. "
1086 "The reference is the column index as printed in the output stream, "
1087 "starting with 1. The index 0 is reserved for the row number.\n"
1088#endif
1089 ;
1090 cout << endl;
1091}
1092
1093void PrintHelp()
1094{
1095#ifdef HAVE_ROOT
1096 cout <<
1097 "\n\n"
1098 "Examples:\n"
1099 "In --root mode, fitsdump support TFormula's syntax for all columns and the filter "
1100 "You can then refer to a column or a (single) index of the column just by its name "
1101 "If the index is omitted, 0 is assumed. Note that the [x:y] syntax in this mode is "
1102 "not supported\n"
1103 "\n"
1104 " fitsdump Zd --filter=\"[0]>20 && cos([1])*TMath::RadToDeg()<45\"\n"
1105 "\n"
1106 "The columns can also be addressed with their names\n"
1107 "\n"
1108 " fitsdump -r \"(Zd+Err)*TMath::DegToRad()\" --filter=\"[0]<25 && [1]<0.05\"\n"
1109 "\n"
1110 "is identical to\n"
1111 "\n"
1112 " fitsdump -r \"(Zd[0]+Err[0])*TMath::DegToRad()\" --filter=\"[0]<25 && [1]<0.05\"\n"
1113 "\n"
1114 "A special placeholder exists for the row number\n"
1115 "\n"
1116 " fitsdump -r \"#\" --filter=\"#>10 && #<100\"\n"
1117 "\n"
1118 "To format a single column you can do\n"
1119 "\n"
1120 " fitsdump col1 -%.1f col2 -%d\n"
1121 "\n"
1122 "A special format is provided converting to 'human readable format'\n"
1123 "\n"
1124 " fitsdump col1 -%h\n"
1125 "\n";
1126 cout << endl;
1127#endif
1128}
1129
1130
1131void SetupConfiguration(Configuration& conf)
1132{
1133 po::options_description configs("Fitsdump options");
1134 configs.add_options()
1135 ("filecontent", po_switch(), "List the number of tables in the file, along with their name")
1136 ("header,h", po_switch(), "Dump header of given table")
1137 ("list,l", po_switch(), "List all tables and columns in file")
1138 ("fitsfile", var<string>()
1139#if BOOST_VERSION >= 104200
1140 ->required()
1141#endif
1142 , "Name of FITS file")
1143 ("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")
1144 ("outfile,o", var<string>("-"), "Name of output file (-:/dev/stdout)")
1145 ("precision,p", var<int>(20), "Precision of ofstream")
1146 ("stat,s", po_switch(), "Perform statistics instead of dump")
1147 ("minmax,m", po_switch(), "Calculates min and max of data")
1148 ("nozero,z", po_switch(), "skip 0 values for stats")
1149 ("fixed", po_switch(), "Switch output stream to floating point values in fixed-point notation")
1150 ("scientific", po_switch(), "Switch output stream to floating point values in scientific notation")
1151 ("%,%", vars<string>(), "Format for the output (currently not available in root-mode)")
1152 ("force", po_switch(), "Force reading the fits file even if END key is missing")
1153 ("first", var<size_t>(size_t(0)), "First number of row to read")
1154 ("limit", var<size_t>(size_t(0)), "Limit for the maximum number of rows to read (0=unlimited)")
1155 ("tablename,t", var<string>(""), "Name of the table to open. If not specified, first binary table is opened")
1156#ifdef HAVE_ROOT
1157 ("root,r", po_switch(), "Enable root mode")
1158 ("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)")
1159#endif
1160 ;
1161
1162 po::positional_options_description p;
1163 p.add("fitsfile", 1); // The first positional options
1164 p.add("col", -1); // All others
1165
1166 conf.AddOptions(configs);
1167 conf.SetArgumentPositions(p);
1168}
1169
1170int main(int argc, const char** argv)
1171{
1172 Configuration conf(argv[0]);
1173 conf.SetPrintUsage(PrintUsage);
1174 SetupConfiguration(conf);
1175
1176 if (!conf.DoParse(argc, argv, PrintHelp))
1177 return -1;
1178
1179 if (!conf.Has("fitsfile"))
1180 {
1181 cerr << "Filename required." << endl;
1182 return -1;
1183 }
1184
1185 FitsDumper loader(conf.Get<string>("fitsfile"), conf.Get<string>("tablename"));
1186 if (!loader)
1187 {
1188 cerr << "ERROR - Opening " << conf.Get<string>("fitsfile");
1189 cerr << " failed: " << strerror(errno) << endl;
1190 return -1;
1191 }
1192
1193 return loader.Exec(conf);
1194}
Note: See TracBrowser for help on using the repository browser.