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

Last change on this file since 12912 was 12912, checked in by tbretz, 9 years ago
Transfer precision to stream which is used for formattting
File size: 20.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 "Time.h"
18#include "externals/fits.h"
19
20#ifdef HAVE_ROOT
21#include "TFormula.h"
22#endif
23
24using namespace std;
25
26struct MyColumn
27{
28 string name;
29
30 fits::Table::Column col;
31
32 uint32_t first;
33 uint32_t last;
34
35 void *ptr;
36};
37
38struct minMaxStruct
39{
40 double min;
41 double max;
42 long double average;
43 long double squared;
44 long numValues;
45 minMaxStruct() : min(FLT_MAX), max(-FLT_MAX), average(0), squared(0), numValues(0) { }
46
47 void add(double val)
48 {
49 average += val;
50 squared += val*val;
51
52 if (val<min)
53 min = val;
54
55 if (val>max)
56 max = val;
57
58 numValues++;
59 }
60};
61
62
63class FitsDumper : public fits
64{
65private:
66 string fFilename;
67
68 // Convert CCfits::ValueType into a human readable string
69 string ValueTypeToStr(char type) const;
70
71 /// Lists all columns of an open file
72 void List();
73 void ListHeader(const string& filename);
74 void ListKeywords(ostream &);
75
76 vector<MyColumn> InitColumns(vector<string> list);
77
78 double GetDouble(const MyColumn &, size_t) const;
79
80 ///Display the selected columns values VS time
81 void Dump(ofstream &, const vector<MyColumn> &, const string &, size_t, size_t, const string &);
82 void DumpMinMax(ofstream &, const vector<MyColumn> &, size_t, size_t, bool);
83 void DumpStats(ofstream &, const vector<MyColumn> &, size_t, size_t);
84
85public:
86 FitsDumper(const string &fname);
87
88 ///Configures the fitsLoader from the config file and/or command arguments.
89 int Exec(Configuration& conf);
90};
91
92// --------------------------------------------------------------------------
93//
94//! Constructor
95//! @param out
96//! the ostream where to redirect the outputs
97//
98FitsDumper::FitsDumper(const string &fname) : fits(fname), fFilename(fname)
99{
100}
101
102string FitsDumper::ValueTypeToStr(char type) const
103{
104 switch (type)
105 {
106 case 'L': return "bool(8)";
107 case 'A': return "char(8)";
108 case 'B': return "byte(8)";
109 case 'I': return "short(16)";
110 case 'J': return "int(32)";
111 case 'K': return "int(64)";
112 case 'E': return "float(32)";
113 case 'D': return "double(64)";
114 default:
115 return "unknown";
116 }
117}
118
119void FitsDumper::List()
120{
121 const fits::Table::Keys &fKeyMap = GetKeys();
122 const fits::Table::Columns &fColMap = GetColumns();
123
124 cout << "\nFile: " << fFilename << "\n";
125
126 cout << " " << fKeyMap.find("EXTNAME")->second.value << " [";
127 cout << fKeyMap.find("NAXIS2")->second.value << "]\n";
128
129 for (auto it = fColMap.begin(); it != fColMap.end(); it++)
130 {
131 cout << " " << it->first << "[" << it->second.num << "] (" << it->second.unit << ":" << ValueTypeToStr(it->second.type) << ") ";
132 for (auto jt = fKeyMap.begin(); jt != fKeyMap.end(); jt++)
133 if (jt->second.value == it->first)
134 cout << "/ " << jt->second.comment << endl;
135 }
136
137 cout << endl;
138}
139
140void FitsDumper::ListKeywords(ostream &fout)
141{
142 const fits::Table::Keys &fKeyMap = GetKeys();
143
144 for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++)
145 {
146 fout << "## " << ::left << setw(8) << it->first << "= ";
147
148 if (it->second.type=='T')
149 fout << ::left << setw(20) << ("'"+it->second.value+"'");
150 else
151 fout << ::right << setw(20) << it->second.value;
152
153 if (it->second.comment.size()>0)
154 fout << " / " << it->second.comment;
155 fout << '\n';
156 }
157
158 fout << flush;
159}
160
161void FitsDumper::ListHeader(const string& filename)
162{
163 ofstream fout(filename=="-"?"/dev/stdout":filename);
164 if (!fout)
165 {
166 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
167 return;
168 }
169
170 const fits::Table::Keys &fKeyMap = GetKeys();
171
172 fout << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n";
173 if (fKeyMap.find("COMMENT") != fKeyMap.end())
174 fout << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n";
175
176 ListKeywords(fout);
177 fout << endl;
178}
179
180vector<MyColumn> FitsDumper::InitColumns(vector<string> names)
181{
182 static const boost::regex expr("([[:word:].]+)(\\[([[:digit:]]+)?(:)?([[:digit:]]+)?\\])?");
183
184 const fits::Table::Columns &fColMap = GetColumns();
185
186 if (names.size()==0)
187 for (auto it=fColMap.begin(); it!=fColMap.end(); it++)
188 if (it->second.num>0)
189 names.push_back(it->first);
190
191 vector<MyColumn> vec;
192
193 for (auto it=names.begin(); it!=names.end(); it++)
194 {
195 boost::smatch what;
196 if (!boost::regex_match(*it, what, expr, boost::match_extra))
197 {
198 cerr << "Couldn't parse expression '" << *it << "' " << endl;
199 return vector<MyColumn>();
200 }
201
202 const string name = what[1];
203
204 const auto iter = fColMap.find(name);
205 if (iter==fColMap.end())
206 {
207 cerr << "ERROR - Column '" << name << "' not found in table." << endl;
208 return vector<MyColumn>();
209 }
210
211 const fits::Table::Column &col = iter->second;
212
213 const string val0 = what[3];
214 const string delim = what[4];
215 const string val1 = what[5];
216
217 const uint32_t first = val0.empty() ? 0 : atoi(val0.c_str());
218 const uint32_t last = (val0.empty() && delim.empty()) ? col.num-1 : (val1.empty() ? first : atoi(val1.c_str()));
219
220 if (first>=col.num)
221 {
222 cerr << "ERROR - First index " << first << " for column " << name << " exceeds number of elements " << col.num << endl;
223 return vector<MyColumn>();
224 }
225
226 if (last>=col.num)
227 {
228 cerr << "ERROR - Last index " << last << " for column " << name << " exceeds number of elements " << col.num << endl;
229 return vector<MyColumn>();
230 }
231
232 if (first>last)
233 {
234 cerr << "ERROR - Last index " << last << " for column " << name << " exceeds first index " << first << endl;
235 return vector<MyColumn>();
236 }
237
238 MyColumn mycol;
239
240 mycol.name = name;
241 mycol.col = col;
242 mycol.first = first;
243 mycol.last = last;
244
245 vec.push_back(mycol);
246 }
247
248 for (auto it=vec.begin(); it!=vec.end(); it++)
249 it->ptr = SetPtrAddress(it->name);
250
251 return vec;
252}
253
254double FitsDumper::GetDouble(const MyColumn &it, size_t i) const
255{
256 switch (it.col.type)
257 {
258 case 'A':
259 return reinterpret_cast<const char*>(it.ptr)[i];
260
261 case 'L':
262 return reinterpret_cast<const bool*>(it.ptr)[i];
263
264 case 'B':
265 return (unsigned int)reinterpret_cast<const uint8_t*>(it.ptr)[i];
266
267 case 'I':
268 return reinterpret_cast<const int16_t*>(it.ptr)[i];
269
270 case 'J':
271 return reinterpret_cast<const int32_t*>(it.ptr)[i];
272
273 case 'K':
274 return reinterpret_cast<const int64_t*>(it.ptr)[i];
275
276 case 'E':
277 return reinterpret_cast<const float*>(it.ptr)[i];
278
279 case 'D':
280 return reinterpret_cast<const double*>(it.ptr)[i];
281 }
282
283 return 0;
284}
285
286// --------------------------------------------------------------------------
287//
288//! Perform the actual dump, based on the current parameters
289//
290void FitsDumper::Dump(ofstream &fout, const vector<MyColumn> &cols, const string &filter, size_t first, size_t limit, const string &filename)
291{
292 const fits::Table::Keys &fKeyMap = GetKeys();
293
294 TFormula select;
295 if (!filter.empty() && select.Compile(filter.c_str()))
296 throw runtime_error("Syntax Error: TFormula::Compile failed for '"+filter+"'");
297
298 fout << "## --------------------------------------------------------------------------\n";
299 fout << "## Fits file: \t" << fFilename << '\n';
300 if (filename!="-")
301 fout << "## File: \t" << filename << '\n';
302 fout << "## Table: \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
303 fout << "## NumRows: \t" << GetInt("NAXIS2") << '\n';
304 fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
305 if (!filter.empty())
306 fout << "## Selection: \t" << select.GetExpFormula() << '\n';
307 fout << "## --------------------------------------------------------------------------\n";
308 ListKeywords(fout);
309 fout << "## --------------------------------------------------------------------------\n";
310 fout << "#\n";
311
312 size_t num = 0;
313 for (auto it=cols.begin(); it!=cols.end(); it++)
314 {
315 fout << "# " << it->name;
316
317 if (it->first==it->last)
318 {
319 if (it->first!=0)
320 fout << "[" << it->first << "]";
321 }
322 else
323 fout << "[" << it->first << ":" << it->last << "]";
324
325 fout << ": " << it->col.unit << '\n';
326
327 num += it->last-it->first+1;
328 }
329 fout << "#" << endl;
330
331 // -----------------------------------------------------------------
332
333 vector<Double_t> data(num);
334
335 const size_t last = limit ? first + limit : size_t(-1);
336
337 while (GetRow(first++))
338 {
339 const size_t row = GetRow();
340 if (row==GetNumRows() || row==last)
341 break;
342
343 size_t p = 0;
344
345 ostringstream out;
346 out.precision(fout.precision());
347 for (auto it=cols.begin(); it!=cols.end(); it++)
348 {
349 string msg;
350 for (uint32_t i=it->first; i<=it->last; i++, p++)
351 {
352 switch (it->col.type)
353 {
354 case 'A':
355 msg += reinterpret_cast<const char*>(it->ptr)[i];
356 break;
357 case 'B':
358 out << (unsigned int)reinterpret_cast<const unsigned char*>(it->ptr)[i] << " ";
359 break;
360 case 'L':
361 out << reinterpret_cast<const bool*>(it->ptr)[i] << " ";
362 break;
363 case 'I':
364 out << reinterpret_cast<const int16_t*>(it->ptr)[i] << " ";
365 break;
366 case 'J':
367 out << reinterpret_cast<const int32_t*>(it->ptr)[i] << " ";
368 break;
369 case 'K':
370 out << reinterpret_cast<const int64_t*>(it->ptr)[i] << " ";
371 break;
372 case 'E':
373 out << reinterpret_cast<const float*>(it->ptr)[i] << " ";
374 break;
375 case 'D':
376 out << reinterpret_cast<const double*>(it->ptr)[i] << " ";
377 break;
378 default:
379 ;
380 }
381
382 if (!filter.empty())
383 data[p] = GetDouble(*it, i);
384 }
385
386 if (it->col.type=='A')
387 out << "'" << msg << "' ";
388 }
389
390 if (!filter.empty() && select.EvalPar(0, data.data())<0.5)
391 continue;
392
393 fout << out.str() << endl;
394 }
395}
396
397void FitsDumper::DumpMinMax(ofstream &fout, const vector<MyColumn> &cols, size_t first, size_t limit, bool fNoZeroPlease)
398{
399 vector<minMaxStruct> statData(cols.size());
400
401 // Loop over all columns in our list of requested columns
402 const size_t last = limit ? first + limit : size_t(-1);
403
404 while (GetRow(first++))
405 {
406 const size_t row = GetRow();
407 if (row==GetNumRows() || row==last)
408 break;
409
410 auto statsIt = statData.begin();
411
412 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
413 {
414 if ((it->name=="UnixTimeUTC" || it->name=="PCTime") && it->first==0 && it->last==1)
415 {
416 const uint32_t *val = reinterpret_cast<const uint32_t*>(it->ptr);
417 if (fNoZeroPlease && val[0]==0 && val[1]==0)
418 continue;
419
420 statsIt->add(Time(val[0], val[1]).Mjd());
421 continue;
422 }
423
424 for (uint32_t i=it->first; i<=it->last; i++)
425 {
426 const double cValue = GetDouble(*it, i);
427
428 if (fNoZeroPlease && cValue == 0)
429 continue;
430
431 statsIt->add(cValue);
432 }
433 }
434 }
435
436 // okay. So now I've got ALL the data, loaded.
437 // let's do the summing and averaging in a safe way (i.e. avoid overflow
438 // of variables as much as possible)
439 auto statsIt = statData.begin();
440 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
441 {
442 fout << "\n[" << it->name << ':' << it->first;
443 if (it->first!=it->last)
444 fout << ':' << it->last;
445 fout << "]\n";
446
447 if (statsIt->numValues == 0)
448 {
449 fout << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
450 continue;
451 }
452
453 const long &num = statsIt->numValues;
454
455 long double &avg = statsIt->average;
456 long double &rms = statsIt->squared;
457
458 avg /= num;
459 rms = sqrt(rms/num - avg*avg);
460
461 fout << "Min: " << statsIt->min << '\n';
462 fout << "Max: " << statsIt->max << '\n';
463 fout << "Avg: " << avg << '\n';
464 fout << "Rms: " << rms << endl;
465 }
466}
467
468template<typename T>
469void displayStats(vector<char> &array, ofstream& out)
470{
471 const size_t numElems = array.size()/sizeof(T);
472 if (numElems == 0)
473 {
474 out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
475 return;
476 }
477
478 T *val = reinterpret_cast<T*>(array.data());
479
480 sort(val, val+numElems);
481
482 out << "Min: " << double(val[0]) << '\n';
483 out << "Max: " << double(val[numElems-1]) << '\n';
484
485 if (numElems>2)
486 {
487 if (numElems%2 == 0)
488 out << "Med: " << (double(val[numElems/2]) + double(val[numElems/2+1]))/2 << '\n';
489 else
490 out << "Med: " << double(val[numElems/2+1]) << '\n';
491 }
492
493 long double avg = 0;
494 long double rms = 0;
495 for (uint32_t i=0;i<numElems;i++)
496 {
497 avg += double(val[i]);
498 rms += double(val[i])*double(val[i]);
499 }
500
501 avg /= numElems;
502 rms = sqrt(rms/numElems - avg*avg);
503
504 out << "Avg: " << avg << '\n';
505 out << "Rms: " << rms << endl;
506
507}
508
509void FitsDumper::DumpStats(ofstream &fout, const vector<MyColumn> &cols, size_t first, size_t limit)
510{
511 // Loop over all columns in our list of requested columns
512 vector<vector<char>> statData;
513
514 const size_t num = limit==0 || GetNumRows()<limit ? GetNumRows() : limit;
515
516 for (auto it=cols.begin(); it!=cols.end(); it++)
517 statData.push_back(vector<char>(it->col.size*num*(it->last-it->first+1)));
518
519 // Loop over all columns in our list of requested columns
520 const size_t last = limit ? first + limit : size_t(-1);
521
522 while (GetRow(first++))
523 {
524 const size_t row = GetRow();
525 if (row==GetNumRows() || row==last)
526 break;
527
528
529 auto statsIt = statData.begin();
530 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
531 {
532 const char *src = reinterpret_cast<const char*>(it->ptr);
533 const size_t sz = (it->last-it->first+1)*it->col.size;
534 memcpy(statsIt->data()+row*sz, src+it->first*it->col.size, sz);
535 }
536 }
537
538 auto statsIt = statData.begin();
539 for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
540 {
541 fout << "\n[" << it->name << ':' << it->first;
542 if (it->last!=it->first)
543 fout << ':' << it->last;
544 fout << "]\n";
545
546 switch (it->col.type)
547 {
548 case 'L':
549 displayStats<bool>(*statsIt, fout);
550 break;
551 case 'B':
552 displayStats<char>(*statsIt, fout);
553 break;
554 case 'I':
555 displayStats<int16_t>(*statsIt, fout);
556 break;
557 case 'J':
558 displayStats<int32_t>(*statsIt, fout);
559 break;
560 case 'K':
561 displayStats<int64_t>(*statsIt, fout);
562 break;
563 case 'E':
564 displayStats<float>(*statsIt, fout);
565 break;
566 case 'D':
567 displayStats<double>(*statsIt, fout);
568 break;
569 default:
570 ;
571 }
572 }
573}
574
575// --------------------------------------------------------------------------
576//
577//! Retrieves the configuration parameters
578//! @param conf
579//! the configuration object
580//
581int FitsDumper::Exec(Configuration& conf)
582{
583 if (conf.Get<bool>("list"))
584 List();
585
586 if (conf.Get<bool>("header"))
587 ListHeader(conf.Get<string>("outfile"));
588
589
590 if (conf.Get<bool>("header") || conf.Get<bool>("list"))
591 return 1;
592
593 // ------------------------------------------------------------
594
595 if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
596 {
597 cerr << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl;
598 return -1;
599 }
600 if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
601 {
602 cerr << "Invalid combination of options: nozero only works with minmax. Aborting" << endl;
603 return -1;
604 }
605
606 // ------------------------------------------------------------
607
608 const string filename = conf.Get<string>("outfile");
609
610 ofstream fout(filename=="-"?"/dev/stdout":filename);
611 if (!fout)
612 {
613 cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
614 return false;
615 }
616 fout.precision(conf.Get<int>("precision"));
617
618 const vector<MyColumn> cols = InitColumns(conf.Vec<string>("col"));
619 if (cols.size()==0)
620 return false;
621
622 const string filter = conf.Has("filter") ? conf.Get<string>("filter") : "";
623
624 const size_t first = conf.Get<size_t>("first");
625 const size_t limit = conf.Get<size_t>("limit");
626
627 if (conf.Get<bool>("minmax"))
628 {
629 DumpMinMax(fout, cols, first, limit, conf.Get<bool>("nozero"));
630 return 0;
631 }
632
633 if (conf.Get<bool>("stat"))
634 {
635 DumpStats(fout, cols, first, limit);
636 return 0;
637 }
638
639 Dump(fout, cols, filter, first, limit, filename);
640
641 return 0;
642}
643
644void PrintUsage()
645{
646 cout <<
647 "fitsdump is a tool to dump data from a FITS table as ascii.\n"
648 "\n"
649 "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
650 " or: fitsdump [OPTIONS]\n";
651 cout << endl;
652}
653
654void PrintHelp()
655{
656 //
657}
658
659void SetupConfiguration(Configuration& conf)
660{
661 po::options_description configs("Fitsdump options");
662 configs.add_options()
663 ("fitsfile,f", var<string>()
664#if BOOST_VERSION >= 104200
665 ->required()
666#endif
667 , "Name of FITS file")
668 ("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")
669 ("outfile,o", var<string>("/dev/stdout"), "Name of output file (-:/dev/stdout)")
670 ("precision,p", var<int>(20), "Precision of ofstream")
671 ("list,l", po_switch(), "List all tables and columns in file")
672 ("header,h", po_switch(), "Dump header of given table")
673 ("stat,s", po_switch(), "Perform statistics instead of dump")
674 ("minmax,m", po_switch(), "Calculates min and max of data")
675 ("nozero,z", po_switch(), "skip 0 values for stats")
676 ("force", po_switch(), "Force reading the fits file even if END key is missing")
677 ("first", var<size_t>(size_t(0)), "First number of row to read")
678 ("limit", var<size_t>(size_t(0)), "Limit for the maximum number of rows to read (0=unlimited)")
679#ifdef HAVE_ROOT
680 ("filter,r", var<string>(""), "Filter to restrict the selection of events (e.g. '[0]>10 && [0]<20'; does not work with stat and minmax yet)")
681#endif
682 ;
683
684 po::positional_options_description p;
685 p.add("fitsfile", 1); // The first positional options
686 p.add("col", -1); // All others
687
688 conf.AddOptions(configs);
689 conf.SetArgumentPositions(p);
690}
691
692int main(int argc, const char** argv)
693{
694 Configuration conf(argv[0]);
695 conf.SetPrintUsage(PrintUsage);
696 SetupConfiguration(conf);
697
698 if (!conf.DoParse(argc, argv, PrintHelp))
699 return -1;
700
701 if (!conf.Has("fitsfile"))
702 {
703 cerr << "Filename required." << endl;
704 return -1;
705 }
706
707 FitsDumper loader(conf.Get<string>("fitsfile"));
708 if (!loader)
709 {
710 cerr << "ERROR - Opening " << conf.Get<string>("fitsfile");
711 cerr << " failed: " << strerror(errno) << endl;
712 return -1;
713 }
714
715 return loader.Exec(conf);
716}
Note: See TracBrowser for help on using the repository browser.