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

Last change on this file since 12912 was 12912, checked in by tbretz, 8 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.