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

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