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

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