Changeset 12884


Ignore:
Timestamp:
02/09/12 18:36:02 (13 years ago)
Author:
tbretz
Message:
Fixed a bug which caused the first line incorrectly displayed; at this stage did a lot of code cleanup, simplifications, unifications of identical code and removed the plotting stuff to make it better maintainable.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/FACT++/src/fitsdump.cc

    r12883 r12884  
    88#include "Configuration.h"
    99
     10#include <float.h>
     11
    1012#include <map>
    1113#include <fstream>
    1214
     15#include <boost/regex.hpp>
     16
     17#include "Time.h"
    1318#include "externals/fits.h"
    1419
    15 //#define PLOTTING_PLEASE
    16 
    17 #ifdef PLOTTING_PLEASE
    18 #include <QPen>
    19 #include <QtGui>
    20 #include <QApplication>
    21 
    22 #include <qwt_plot.h>
    23 #include <qwt_plot_grid.h>
    24 #include <qwt_plot_curve.h>
    25 #include <qwt_plot_zoomer.h>
    26 #include <qwt_legend.h>
    27 #include <qwt_scale_draw.h>
    28 
    29 #endif
    30 #include "Time.h"
    31 
    3220using namespace std;
    3321
    34 
    35 #ifdef PLOTTING_PLEASE
    36 class TimeScaleDraw: public QwtScaleDraw
    37 {
    38 public:
    39     virtual QwtText label(double v) const
    40     {
    41         Time t(v);
    42         string time = t.GetAsStr("%H:%M:%S%F");
    43         while (time[time.size()-1] == '0' && time.size() > 2)
    44         {
    45             time = time.substr(0, time.size()-1);
    46         }
    47         return QwtText(time.c_str());
    48     }
     22struct MyColumn
     23{
     24    string name;
     25
     26    fits::Table::Column col;
     27
     28    uint32_t first;
     29    uint32_t last;
     30
     31    void *ptr;
    4932};
    50 #endif
    51 
    52 class FitsDumper
    53 {
    54 public:
    55     FitsDumper();
    56     ~FitsDumper();
    57 
     33
     34struct minMaxStruct
     35{
     36    double min;
     37    double max;
     38    long double average;
     39    long double squared;
     40    long numValues;
     41    minMaxStruct() : min(FLT_MAX), max(-FLT_MAX), average(0), squared(0), numValues(0) { }
     42
     43    void add(double val)
     44    {
     45        average += val;
     46        squared += val*val;
     47
     48        if (val<min)
     49            min = val;
     50
     51        if (val>max)
     52            max = val;
     53
     54        numValues++;
     55    }
     56};
     57
     58
     59class FitsDumper : public fits
     60{
    5861private:
    59     fits* fFile;
    60     bool fDotsPlease;
    61     bool fNoZeroPlease;
    6262    string fFilename;
    63 
    64     fits::Table::Columns fColMap;
    65     fits::Table::Keys fKeyMap;
    6663
    6764    // Convert CCfits::ValueType into a human readable string
    6865    string ValueTypeToStr(char type) const;
    69 
    70     // Convert CCfits::ValueType into a number of associated bytes
    71     int    ValueTypeToSize(char type) const;
    72 
    73     /// Calculate the buffer size required to read a row of the fits table, as well as the offsets to each column
    74 //    vector<int> CalculateOffsets() const;
    75 
    76     template<class T>
    77         T PtrToValue(const unsigned char* &ptr) const;
    78 //    template<class T>
    79 //        double PtrToDouble(const unsigned char *ptr) const;
    80 //    double PtrToDouble(const unsigned char *ptr, CCfits::ValueType type) const;
    81 
    82     /// Write a single row of the selected data
    83  //   int  WriteRow(ostream &, const vector<MyColumn*> &, const vector<int> &, unsigned char *, const vector<pair<int, int> >&) const;
    84 
    85     bool OpenFile(const string &, bool);        /// Open a file
    86     bool OpenTable(const string &);       /// Open a table
    8766
    8867    /// Lists all columns of an open file
    … …  
    9170    void ListKeywords(ostream &);
    9271
    93     bool separateColumnsFromRanges(const vector<string>& list,
    94                                    vector<pair<int, int> >& ranges,
    95                                    vector<string>& listNamesOnly);
    96     /// Perform the dumping, based on the current dump list
    97     bool Dump(const string &, const vector<string> &list, int);
     72    vector<MyColumn> InitColumns(const vector<string>& list);
     73
    9874    ///Display the selected columns values VS time
    99 #ifdef PLOTTING_PLEASE
    100     int doCurvesDisplay( const vector<string> &list, const string& tableName);
    101 #endif
    102     int doMinMaxPlease(const string& filename, const vector<string>& list, int precision);
    103     int doStatsPlease(const string &filename, const vector<string>& list, int precision);
    104 //    void doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), true);
    105     //    bool Plot(const vector<string> &list);
     75    void Dump(ofstream &, const vector<MyColumn> &, const string &);
     76    void DumpMinMax(ofstream &, const vector<MyColumn> &, bool);
     77    void DumpStats(ofstream &, const vector<MyColumn> &);
    10678
    10779public:
     80    FitsDumper(const string &fname);
     81
    10882    ///Configures the fitsLoader from the config file and/or command arguments.
    109     int ExecConfig(Configuration& conf);
     83    int Exec(Configuration& conf);
    11084};
    11185
    … …  
    11690//!        the ostream where to redirect the outputs
    11791//
    118 FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false), fNoZeroPlease(false)
    119 {
    120 }
    121 
    122 // --------------------------------------------------------------------------
    123 //
    124 //! Destructor
    125 //
    126 FitsDumper::~FitsDumper()
    127 {
    128     if (fFile)
    129         delete fFile;
    130 }
    131 
     92FitsDumper::FitsDumper(const string &fname) : fits(fname), fFilename(fname)
     93{
     94}
    13295
    13396string FitsDumper::ValueTypeToStr(char type) const
    … …  
    148111}
    149112
    150 int FitsDumper::ValueTypeToSize(char type) const
    151 {
    152     switch (type)
    153     {
    154         case 'A': return sizeof(int8_t);
    155         case 'L': return sizeof(uint8_t);
    156         case 'B': return sizeof(int8_t);
    157         case 'I': return sizeof(int16_t);
    158         case 'J': return sizeof(int32_t);
    159         case 'K': return sizeof(int64_t);
    160         case 'E': return sizeof(float);
    161         case 'D': return sizeof(double);
    162     default:
    163         return 0;
    164     }
    165 }
    166 
    167 template<class T>
    168 T FitsDumper::PtrToValue(const unsigned char* &ptr) const
    169 {
    170     T t;
    171     reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast<unsigned char*>(&t));
    172     ptr += sizeof(T);
    173 
    174     return t;
    175 }
    176 // --------------------------------------------------------------------------
    177 //
    178 //! Loads the fits file based on the current parameters
    179 //
    180 bool FitsDumper::OpenFile(const string &filename, bool force)
    181 {
    182     if (fFile)
    183     {
    184         fFile->close();
    185         delete fFile;
    186     }
    187 
    188     try {
    189         fFile = new fits(filename, force);
    190     }
    191     catch (std::runtime_error e)
    192     {
    193         cout << "Something went wrong while trying to open " << filename;
    194         cout << ": " << e.what() << " Aborting dump." << endl;
    195         return false;
    196     }
    197     fFilename = filename;
    198 
    199     const fits::Table::Columns& tCols = fFile->GetColumns();
    200 
    201     for (auto it=tCols.begin(); it != tCols.end(); it++)
    202         fColMap.insert(*it);
    203 
    204     const fits::Table::Keys& tkeys = fFile->GetKeys();
    205 
    206     for (auto it=tkeys.begin(); it != tkeys.end(); it++)
    207         fKeyMap.insert(*it);
    208 
    209     return true;
    210 }
    211 
    212 bool FitsDumper::OpenTable(const string &)
    213 {
    214     if (!fFile)
    215     {
    216         cerr << "No file open." << endl;
    217         return false;
    218     }
    219 
    220 
    221     return true;
    222 }
    223 
    224 
    225113void FitsDumper::List()
    226114{
    227     if (!fFile)
    228     {
    229         cerr << "No file open." << endl;
    230         return;
    231     }
     115    const fits::Table::Keys    &fKeyMap = GetKeys();
     116    const fits::Table::Columns &fColMap = GetColumns();
    232117
    233118    cout << "\nFile: " << fFilename << "\n";
    … …  
    245130
    246131    cout << endl;
    247     cout << flush;
    248132}
    249133
    250134void FitsDumper::ListKeywords(ostream &out)
    251135{
     136    const fits::Table::Keys &fKeyMap = GetKeys();
     137
    252138    for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++) {
    253139        out << "## " << setw(8) << it->first << " = " << setw(10);
    … …  
    258144void FitsDumper::ListHeader(const string& filename)
    259145{
    260     if (!fFile)
    261     {
    262         cerr << "No table open." << endl;
    263         return;
    264     }
    265146    ofstream out(filename=="-"?"/dev/stdout":filename);
    266147    if (!out)
    … …  
    269150        return;
    270151    }
     152
     153    const fits::Table::Keys &fKeyMap = GetKeys();
    271154
    272155    out << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n";
    … …  
    279162}
    280163
    281 bool FitsDumper::separateColumnsFromRanges(const vector<string>& list,
    282                                vector<pair<int, int> >& ranges,
    283                                vector<string>& listNamesOnly)
    284 {
    285     for (vector<string>::const_iterator it=list.begin(); it!=list.end(); it++)
    286     {
    287         string columnNameOnly = *it;
    288         unsigned long bracketIndex0 = columnNameOnly.find_first_of('[');
    289         unsigned long bracketIndex1 = columnNameOnly.find_first_of(']');
    290         unsigned long colonIndex = columnNameOnly.find_first_of(':');
    291         int columnStart = -1;
    292         int columnEnd = -1;
    293         if (bracketIndex0 != string::npos)
    294         {//there is a range given. Extract the range
    295             if (colonIndex != string::npos)
    296             {//we have a range here
    297                 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, colonIndex-(bracketIndex0+1)).c_str());
    298                 columnEnd = atoi(columnNameOnly.substr(colonIndex+1, bracketIndex1-(colonIndex+1)).c_str());
    299                 columnEnd++;
    300             }
    301             else
    302             {//only a single index there
    303                 columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, bracketIndex1 - (bracketIndex0+1)).c_str());
    304                 columnEnd = columnStart+1;
    305             }
    306             columnNameOnly = columnNameOnly.substr(0, bracketIndex0);
    307         }
    308 
    309         if (fColMap.find(columnNameOnly) == fColMap.end())
    310         {
    311             cerr << "ERROR - Column '" << columnNameOnly << "' not found in table." << endl;
    312             return false;
    313         }
    314         fits::Table::Column& cCol = fColMap.find(columnNameOnly)->second;
    315         if (bracketIndex0 == string::npos)
    316         {//no range given: use the full range
    317             ranges.push_back(make_pair(0, cCol.num));
    318             columnStart = 0;
    319             columnEnd = 1;
    320         }
    321         else
    322         {//use the range extracted earlier
    323             if (columnStart < 0)
    324             {
    325                 cerr << "ERROR - Start range for column " << columnNameOnly << " is less than zero (" << columnStart << "). Aborting" << endl;
    326                 return false;
    327             }
    328             if (columnEnd>1 && columnEnd > (int)(cCol.num))
    329             {
    330                 cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol.num << " vs " << columnEnd << "). Aborting" << endl;
    331                 return false;
    332             }
    333             ranges.push_back(make_pair(columnStart, columnEnd));
    334         }
    335 //        cout << "Will be exporting from " << columnStart << " to " << columnEnd-1 << " for column " << columnNameOnly << endl;
    336         listNamesOnly.push_back(columnNameOnly);
    337     }
    338     return true;
    339 }
     164vector<MyColumn> FitsDumper::InitColumns(const vector<string> &names)
     165{
     166    static const boost::regex expr("([[:word:].]+)(\\[([[:digit:]]+)?(:)?([[:digit:]]+)?\\])?");
     167
     168    const fits::Table::Columns &fColMap = GetColumns();
     169
     170    vector<MyColumn> vec;
     171
     172    for (auto it=names.begin(); it!=names.end(); it++)
     173    {
     174        boost::smatch what;
     175        if (!boost::regex_match(*it, what, expr, boost::match_extra))
     176        {
     177            cerr << "Couldn't parse expression '" << *it << "' " << endl;
     178            return vector<MyColumn>();
     179        }
     180
     181        const string name = what[1];
     182
     183        const auto iter = fColMap.find(name);
     184        if (iter==fColMap.end())
     185        {
     186            cerr << "ERROR - Column '" << name << "' not found in table." << endl;
     187            return vector<MyColumn>();
     188        }
     189
     190        const fits::Table::Column &col = iter->second;
     191
     192        const string val0  = what[3];
     193        const string delim = what[4];
     194        const string val1  = what[5];
     195
     196        const uint32_t first = val0.empty() ? 0 : atoi(val0.c_str());
     197        const uint32_t last  = val0.empty()==delim.empty() ? col.num-1 : (val1.empty() ? first : atoi(val1.c_str()));
     198
     199        if (first>=col.num)
     200        {
     201            cerr << "ERROR - First index " << first << " for column " << name << " exceeds number of elements " << col.num << endl;
     202            return vector<MyColumn>();
     203        }
     204
     205        if (last>=col.num)
     206        {
     207            cerr << "ERROR - Last index " << last << " for column " << name << " exceeds number of elements " << col.num << endl;
     208            return vector<MyColumn>();
     209        }
     210
     211        if (first>last)
     212        {
     213            cerr << "ERROR - Last index " << last << " for column " << name << " exceeds first index " << first << endl;
     214            return vector<MyColumn>();
     215        }
     216
     217        MyColumn mycol;
     218
     219        mycol.name  = name;
     220        mycol.col   = col;
     221        mycol.first = first;
     222        mycol.last  = last;
     223
     224        vec.push_back(mycol);
     225    }
     226
     227    for (auto it=vec.begin(); it!=vec.end(); it++)
     228        it->ptr = SetPtrAddress(it->name);
     229
     230    return vec;
     231}
     232
    340233// --------------------------------------------------------------------------
    341234//
    342235//! Perform the actual dump, based on the current parameters
    343236//
    344 bool FitsDumper::Dump(const string &filename, const vector<string> &list, int precision)
    345 {
    346 
    347     //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    348     vector<pair<int, int> > ranges;
    349     vector<string> listNamesOnly;
    350 
    351     if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
    352     {
    353         cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
    354         return false;
    355     }
    356 
    357     ofstream out(filename=="-"?"/dev/stdout":filename);
    358     if (!out)
    359     {
    360         cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
    361         return false;
    362     }
    363 
    364     out.precision(precision);
     237void FitsDumper::Dump(ofstream &out, const vector<MyColumn> &cols, const string &filename)
     238{
     239    const fits::Table::Keys &fKeyMap = GetKeys();
    365240
    366241    out << "## --------------------------------------------------------------------------\n";
     242    out << "## Fits file:\t" << fFilename << '\n';
    367243    if (filename!="-")
    368244        out << "## File:    \t" << filename << '\n';
    369245    out << "## Table:   \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
    370         out << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
    371 
    372     out << "## NumRows: \t" << fFile->GetInt("NAXIS2") << '\n';
     246    out << "## NumRows: \t" << GetInt("NAXIS2") << '\n';
     247    out << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
    373248    out << "## --------------------------------------------------------------------------\n";
    374249    ListKeywords(out);
    … …  
    376251    out << "#\n";
    377252
    378     auto rangesIt = ranges.begin();
    379     for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
    380     {
    381         const fits::Table::Column& col = fColMap[*it];
    382 //        const MyColumn *col = static_cast<MyColumn*>(fTable->column()[*it]);
    383         if (rangesIt->first != 0 || rangesIt->second != (int)(col.num))
    384         {
    385             out << "#";
    386             for (int i=rangesIt->first; i<rangesIt->second; i++)
    387                 out << " " << *it << "[" << i << "]";
    388             out << ": " << col.unit;
     253    for (auto it=cols.begin(); it!=cols.end(); it++)
     254    {
     255        out << "# " << it->name;
     256
     257        if (it->first==it->last)
     258        {
     259            if (it->first!=0)
     260                out << "[" << it->first << "]";
    389261        }
    390262        else
    391             out << "# " << *it << "[" << col.num << "]: " << col.unit;
    392 
    393 //        FIXME: retrive the column comment
    394 //        if (!col->comment().empty())
    395 //            out << " (" <<col->comment() << ")";
    396         out << '\n';
     263            out << "[" << it->first << ":" << it->last << "]";
     264
     265        out << ": " << it->col.unit << '\n';
    397266    }
    398267    out << "#" << endl;
    399268
    400 
    401     vector<const void*> ptrs;
    402     for (auto namesIt=listNamesOnly.begin(); namesIt!=listNamesOnly.end(); namesIt++)
    403     {
    404         ptrs.push_back(fFile->SetPtrAddress(*namesIt));
    405     }
    406 
    407     while (fFile->GetNextRow())
    408     {
    409         const size_t row = fFile->GetRow();
    410         if (row==fFile->GetNumRows())
    411             break;
    412 
    413         rangesIt = ranges.begin();
    414         auto ptrsIt=ptrs.begin();
    415 
    416         for (vector<string>::const_iterator in=listNamesOnly.begin(); in!=listNamesOnly.end(); in++, rangesIt++, ptrsIt++)
    417         {
    418             const void *ptr = *ptrsIt;
    419 
    420             const fits::Table::Column &cCol = fColMap.find(*in)->second;
    421 
     269    // -----------------------------------------------------------------
     270
     271    while (GetNextRow())
     272    {
     273        const size_t row = GetRow();
     274        if (row==GetNumRows())
     275            break;
     276
     277        for (auto it=cols.begin(); it!=cols.end(); it++)
     278        {
    422279            string msg;
    423             for (int i=rangesIt->first; i<rangesIt->second; i++)
     280            for (uint32_t i=it->first; i<=it->last; i++)
    424281            {
    425                 switch (cCol.type)
     282                switch (it->col.type)
    426283                {
    427284                case 'A':
    428                     msg += reinterpret_cast<const char*>(ptr)[i];
     285                    msg += reinterpret_cast<const char*>(it->ptr)[i];
    429286                    break;
    430287                case 'B':
    431                     out << (unsigned int)reinterpret_cast<const unsigned char*>(ptr)[i] << " ";
     288                    out << (unsigned int)reinterpret_cast<const unsigned char*>(it->ptr)[i] << " ";
    432289                    break;
    433290                case 'L':
    434                     out << reinterpret_cast<const bool*>(ptr)[i] << " ";
     291                    out << reinterpret_cast<const bool*>(it->ptr)[i] << " ";
    435292                    break;
    436293                case 'I':
    437                     out << reinterpret_cast<const int16_t*>(ptr)[i] << " ";
     294                    out << reinterpret_cast<const int16_t*>(it->ptr)[i] << " ";
    438295                    break;
    439296                case 'J':
    440                     out << reinterpret_cast<const int32_t*>(ptr)[i] << " ";
     297                    out << reinterpret_cast<const int32_t*>(it->ptr)[i] << " ";
    441298                    break;
    442299                case 'K':
    443                     out << reinterpret_cast<const int64_t*>(ptr)[i] << " ";
     300                    out << reinterpret_cast<const int64_t*>(it->ptr)[i] << " ";
    444301                    break;
    445302                case 'E':
    446                     out << reinterpret_cast<const float*>(ptr)[i] << " ";
     303                    out << reinterpret_cast<const float*>(it->ptr)[i] << " ";
    447304                    break;
    448305                case 'D':
    449                     out << reinterpret_cast<const double*>(ptr)[i] << " ";
     306                    out << reinterpret_cast<const double*>(it->ptr)[i] << " ";
    450307                    break;
    451308                default:
    … …  
    454311            }
    455312
    456             if (cCol.type=='A')
     313            if (it->col.type=='A')
    457314                out << "'" << msg << "' ";
    458315        }
    459316        out << endl;
    460317    }
    461     return true;
     318}
     319
     320void FitsDumper::DumpMinMax(ofstream &out, const vector<MyColumn> &cols, bool fNoZeroPlease)
     321{
     322    vector<minMaxStruct> statData(cols.size());
     323
     324    // Loop over all columns in our list of requested columns
     325    while (GetNextRow())
     326    {
     327        const size_t row = GetRow();
     328        if (row==GetNumRows())
     329            break;
     330
     331        auto statsIt = statData.begin();
     332
     333        for (auto in=cols.begin(); in!=cols.end(); in++, statsIt++)
     334        {
     335            if ((in->name=="UnixTimeUTC" || in->name=="PCTime") && in->first==0 && in->last==1)
     336            {
     337                const uint32_t *val = reinterpret_cast<const uint32_t*>(in->ptr);
     338                if (fNoZeroPlease && val[0]==0 && val[1]==0)
     339                    continue;
     340
     341                statsIt->add(Time(val[0], val[1]).Mjd());
     342                continue;
     343            }
     344
     345            for (uint32_t i=in->first; i<=in->last; i++)
     346            {
     347                double cValue = 0;
     348                switch (in->col.type)
     349                {
     350                case 'L':
     351                        cValue = reinterpret_cast<const bool*>(in->ptr)[i];
     352                        break;
     353                case 'B':
     354                        cValue = reinterpret_cast<const int8_t*>(in->ptr)[i];
     355                        break;
     356                case 'I':
     357                        cValue = reinterpret_cast<const int16_t*>(in->ptr)[i];
     358                        break;
     359                case 'J':
     360                        cValue = reinterpret_cast<const int32_t*>(in->ptr)[i];
     361                        break;
     362                case 'K':
     363                        cValue = reinterpret_cast<const int64_t*>(in->ptr)[i];
     364                        break;
     365                case 'E':
     366                        cValue = reinterpret_cast<const float*>(in->ptr)[i];
     367                        break;
     368                case 'D':
     369                        cValue = reinterpret_cast<const double*>(in->ptr)[i];
     370                        break;
     371                default:
     372                    ;
     373                }
     374
     375                if (fNoZeroPlease && cValue == 0)
     376                    continue;
     377
     378                statsIt->add(cValue);
     379            }
     380        }
     381    }
     382
     383    // okay. So now I've got ALL the data, loaded.
     384    // let's do the summing and averaging in a safe way (i.e. avoid overflow
     385    // of variables as much as possible)
     386    auto statsIt = statData.begin();
     387    for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
     388    {
     389        cout << "\n[" << it->name << ':' << it->first;
     390        if (it->first!=it->last)
     391            cout << ':' << it->last;
     392        cout << "]\n";
     393
     394        if (statsIt->numValues == 0)
     395        {
     396            out << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
     397            continue;
     398        }
     399
     400        const long &num = statsIt->numValues;
     401
     402        long double &avg = statsIt->average;
     403        long double &rms = statsIt->squared;
     404
     405        avg /= num;
     406        rms  = sqrt(rms/num - avg*avg);
     407
     408        out << "Min: " << statsIt->min << '\n';
     409        out << "Max: " << statsIt->max << '\n';
     410        out << "Avg: " << avg << '\n';
     411        out << "Rms: " << rms << endl;
     412    }
     413}
     414
     415template<typename T>
     416void displayStats(vector<char> &array, ofstream& out)
     417{
     418    const size_t numElems = array.size()/sizeof(T);
     419    if (numElems == 0)
     420    {
     421        out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
     422        return;
     423    }
     424
     425    T *val = reinterpret_cast<T*>(array.data());
     426
     427    sort(val, val+numElems);
     428
     429    out << "Min: " << double(val[0]) << '\n';
     430    out << "Max: " << double(val[numElems-1]) << '\n';
     431
     432    if (numElems>2)
     433    {
     434        if (numElems%2 == 0)
     435            out << "Med: " << (double(val[numElems/2]) + double(val[numElems/2+1]))/2 << '\n';
     436        else
     437            out << "Med: " << double(val[numElems/2+1]) << '\n';
     438    }
     439
     440    long double avg = 0;
     441    long double rms = 0;
     442    for (uint32_t i=0;i<numElems;i++)
     443    {
     444        avg += double(val[i]);
     445        rms += double(val[i])*double(val[i]);
     446    }
     447
     448    avg /= numElems;
     449    rms  = sqrt(rms/numElems - avg*avg);
     450
     451    out << "Avg: " << avg << '\n';
     452    out << "Rms: " << rms << endl;
     453
     454}
     455
     456void FitsDumper::DumpStats(ofstream &out, const vector<MyColumn> &cols)
     457{
     458    // Loop over all columns in our list of requested columns
     459    vector<vector<char>> statData;
     460
     461    for (auto in=cols.begin(); in!=cols.end(); in++)
     462        statData.push_back(vector<char>(in->col.size*GetNumRows()*(in->last-in->first+1)));
     463
     464    while (GetNextRow())
     465    {
     466        const size_t row = GetRow();
     467        if (row==GetNumRows())
     468            break;
     469
     470        auto statsIt = statData.begin();
     471        for (auto in=cols.begin(); in!=cols.end(); in++, statsIt++)
     472        {
     473            const char *src = reinterpret_cast<const char*>(in->ptr);
     474            const size_t sz = (in->last-in->first+1)*in->col.size;
     475            memcpy(statsIt->data()+row*sz, src+in->first*in->col.size, sz);
     476        }
     477    }
     478
     479    auto statsIt = statData.begin();
     480    for (auto in=cols.begin(); in!=cols.end(); in++, statsIt++)
     481    {
     482        out << "\n[" << in->name << ':' << in->first;
     483        if (in->last!=in->first)
     484            out << ':' << in->last;
     485        out << "]\n";
     486
     487        switch (in->col.type)
     488        {
     489        case 'L':
     490            displayStats<bool>(*statsIt, out);
     491            break;
     492        case 'B':
     493            displayStats<char>(*statsIt, out);
     494            break;
     495        case 'I':
     496            displayStats<int16_t>(*statsIt, out);
     497            break;
     498        case 'J':
     499            displayStats<int32_t>(*statsIt, out);
     500            break;
     501        case 'K':
     502            displayStats<int64_t>(*statsIt, out);
     503            break;
     504        case 'E':
     505            displayStats<float>(*statsIt, out);
     506            break;
     507        case 'D':
     508            displayStats<double>(*statsIt, out);
     509            break;
     510        default:
     511            ;
     512        }
     513    }
    462514}
    463515
    … …  
    468520//!             the configuration object
    469521//
    470 int FitsDumper::ExecConfig(Configuration& conf)
    471 {
    472     if (conf.Has("fitsfile"))
    473     {
    474         if (!OpenFile(conf.Get<string>("fitsfile"), conf.Get<bool>("force")))
    475             return -1;
    476     }
    477 #ifdef PLOTTING_PLEASE
    478     if (conf.Get<bool>("stat") && conf.Get<bool>("graph"))
    479     {
    480         cout << "Invalid conbination of options: cannot graph stats. Aborting" << endl;
    481         return -1;
    482     }
    483 #endif
    484     if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
    485     {
    486         cout << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl;
    487         return -1;
    488     }
    489 #ifdef PLOTTING_PLEASE
    490     if (conf.Get<bool>("minmax") && conf.Get<bool>("graph"))
    491     {
    492         cout << "Invalid combination of options: cannot graph minmax. Aborting" << endl;
    493         return -1;
    494     }
    495 #endif
    496     if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
    497     {
    498         cout << "Invalid combination of options: nozero only works with minmax. Aborting" << endl;
    499         return -1;
    500     }
    501 
    502     if (conf.Get<bool>("nozero"))
    503     {
    504         fNoZeroPlease = true;
    505     }
    506 
    507     if (conf.Has("tablename"))
    508     {
    509         if (!OpenTable(conf.Get<string>("tablename")))
    510             return -1;
    511     }
    512 
     522int FitsDumper::Exec(Configuration& conf)
     523{
    513524    if (conf.Get<bool>("list"))
    514525        List();
    515     /*
    516     if (conf.Get<bool>("tstart"))
    517     {
    518         doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), true);
    519     }
    520     if (conf.Get<bool>("tstop"))
    521     {
    522         doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), false);
    523     }
    524     */
    525     if (conf.Get<bool>("minmax"))
    526     {
    527         if (!conf.Has("col"))
    528         {
    529             cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
    530             return 0;
    531         }
    532         doMinMaxPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
    533         return 0;
    534     }
    535 
    536     if (conf.Get<bool>("stat"))
    537     {
    538         if (!conf.Has("col"))
    539         {
    540             cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
    541             return 0;
    542         }
    543         doStatsPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
    544         return 0;
    545     }
    546 
    547 #ifdef PLOTTING_PLEASE
    548     if (conf.Get<bool>("graph"))
    549     {
    550         if (!conf.Has("col"))
    551         {
    552             cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
    553             return 0;
    554         }
    555         if (conf.Get<bool>("dots"))
    556             fDotsPlease = true;
    557         doCurvesDisplay(conf.Get<vector<string>>("col"),
    558                         conf.Get<string>("tablename"));
    559         return 1;
    560     }
    561 #endif
    562526
    563527    if (conf.Get<bool>("header"))
    564528        ListHeader(conf.Get<string>("outfile"));
    565529
     530
    566531    if (conf.Get<bool>("header") || conf.Get<bool>("list"))
    567532        return 1;
    568533
    569     if (conf.Has("outfile"))
    570     {
    571         if (!conf.Has("col"))
    572         {
    573             cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
    574             return 0;
    575         }
    576         if (!Dump(conf.Get<string>("outfile"),
    577                   conf.Get<vector<string>>("col"),
    578                   conf.Get<int>("precision")))
    579             return -1;
    580     }
    581 
     534    // ------------------------------------------------------------
     535
     536    if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
     537    {
     538        cerr << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl;
     539        return -1;
     540    }
     541    if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
     542    {
     543        cerr << "Invalid combination of options: nozero only works with minmax. Aborting" << endl;
     544        return -1;
     545    }
     546
     547    // ------------------------------------------------------------
     548
     549    if (conf.Vec<string>("col").size()==0)
     550    {
     551        cerr << "No columns specifiec." << endl;
     552        return 0;
     553    }
     554
     555    const string filename = conf.Get<string>("outfile");
     556
     557    ofstream out(filename=="-"?"/dev/stdout":filename);
     558    if (!out)
     559    {
     560        cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
     561        return false;
     562    }
     563    out.precision(conf.Get<int>("precision"));
     564
     565    const vector<MyColumn> cols = InitColumns(conf.Vec<string>("col"));
     566    if (cols.size()==0)
     567        return false;
     568
     569    if (conf.Get<bool>("minmax"))
     570    {
     571        DumpMinMax(out, cols, conf.Get<bool>("nozero"));
     572        return 0;
     573    }
     574
     575    if (conf.Get<bool>("stat"))
     576    {
     577        DumpStats(out, cols);
     578        return 0;
     579    }
     580
     581    Dump(out, cols, filename);
    582582
    583583    return 0;
    … …  
    598598    //
    599599}
    600 
    601 struct minMaxStruct {
    602     double min;
    603     double max;
    604     double average;
    605     double squared;
    606     long numValues;
    607     minMaxStruct() : min(1e10), max(-1e10), average(0), numValues(0) {}
    608 };
    609 
    610 int FitsDumper::doMinMaxPlease(const string& filename, const vector<string>& list, int precision)
    611 {
    612     //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    613     vector<pair<int, int> > ranges;
    614     vector<string> names;
    615 
    616     if (!separateColumnsFromRanges(list, ranges, names))
    617     {
    618         cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
    619         return false;
    620     }
    621 
    622     ofstream out(filename=="-"?"/dev/stdout":filename);
    623     if (!out)
    624     {
    625         cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
    626         return false;
    627     }
    628 
    629     out.precision(precision);
    630 
    631     vector<minMaxStruct> statData(names.size());
    632 
    633     vector<const void*> ptrs;
    634     for (auto namesIt=names.begin(); namesIt!=names.end(); namesIt++)
    635     {
    636         ptrs.push_back(fFile->SetPtrAddress(*namesIt));
    637     }
    638 
    639     // Loop over all columns in our list of requested columns
    640     while (fFile->GetNextRow())
    641     {
    642         const size_t row = fFile->GetRow();
    643         if (row==fFile->GetNumRows())
    644             break;
    645 
    646         auto rangesIt = ranges.begin();
    647         auto statsIt  = statData.begin();
    648         auto ptrsIt = ptrs.begin();
    649 
    650         for (auto in=names.begin(); in!=names.end(); in++, rangesIt++, statsIt++, ptrsIt++)
    651         {
    652             const void *ptr = *ptrsIt;
    653 
    654             const fits::Table::Column &cCol = fColMap.find(*in)->second;
    655 
    656             if (*in=="UnixTimeUTC" || *in=="PCTime")
    657             {
    658                 const uint32_t *val = reinterpret_cast<const uint32_t*>(ptr);
    659                 if (fNoZeroPlease && val[0]==0 && val[1]==0)
    660                     continue;
    661 
    662                 const double cValue = Time(val[0], val[1]).Mjd();
    663 
    664                 statsIt->average += cValue;
    665                 statsIt->squared += cValue*cValue;
    666 
    667                 if (cValue < statsIt->min)
    668                     statsIt->min = cValue;
    669 
    670                 if (cValue > statsIt->max)
    671                     statsIt->max = cValue;
    672 
    673                 statsIt->numValues++;
    674 
    675                 continue;
    676             }
    677 
    678             for (int i=rangesIt->first; i<rangesIt->second; i++)
    679             {
    680                 double cValue = 0;
    681                 switch (cCol.type)
    682                 {
    683                 case 'L':
    684                         cValue = reinterpret_cast<const bool*>(ptr)[i];
    685                         break;
    686                 case 'B':
    687                         cValue = reinterpret_cast<const int8_t*>(ptr)[i];
    688                         break;
    689                 case 'I':
    690                         cValue = reinterpret_cast<const int16_t*>(ptr)[i];
    691                         break;
    692                 case 'J':
    693                         cValue = reinterpret_cast<const int32_t*>(ptr)[i];
    694                         break;
    695                 case 'K':
    696                         cValue = reinterpret_cast<const int64_t*>(ptr)[i];
    697                         break;
    698                 case 'E':
    699                         cValue = reinterpret_cast<const float*>(ptr)[i];
    700                         break;
    701                 case 'D':
    702                         cValue = reinterpret_cast<const double*>(ptr)[i];
    703                         break;
    704                 default:
    705                     ;
    706                 }
    707 
    708                 if (fNoZeroPlease && cValue == 0)
    709                     continue;
    710 
    711                 statsIt->average += cValue;
    712                 statsIt->squared += cValue*cValue;
    713 
    714                 if (cValue < statsIt->min)
    715                     statsIt->min = cValue;
    716 
    717                 if (cValue > statsIt->max)
    718                     statsIt->max = cValue;
    719 
    720                 statsIt->numValues++;
    721             }
    722         }
    723     }
    724 
    725     // okay. So now I've got ALL the data, loaded.
    726     // let's do the summing and averaging in a safe way (i.e. avoid overflow
    727     // of variables as much as possible)
    728     auto statsIt = statData.begin();
    729 
    730     for (auto it=names.begin(); it!=names.end(); it++, statsIt++)
    731     {
    732         cout << "[" << *it << "]" << endl;
    733         if (statsIt->numValues == 0)
    734         {
    735             out << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
    736             continue;
    737         }
    738 
    739         double &avg = statsIt->average;
    740         double &rms = statsIt->squared;
    741         long   &num = statsIt->numValues;
    742 
    743         avg /= num;
    744         rms  = sqrt(rms/num - avg*avg);
    745 
    746         out << "Min: " << statsIt->min << '\n';
    747         out << "Max: " << statsIt->max << '\n';
    748         out << "Avg: " << avg << '\n';
    749         out << "Rms: " << rms << endl;
    750     }
    751 
    752     /*
    753     vector<pair<char, char*> > columnsData;
    754     vector<minMaxStruct> statData;
    755     int numRows = fFile->GetInt("NAXIS2");
    756     auto rangesIt = ranges.begin();
    757     for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
    758     {
    759         fits::Table::Column& cCol = fColMap.find(*it)->second;
    760         columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
    761 //        minMaxStuct initData;
    762         statData.push_back(minMaxStruct());
    763         fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
    764     }
    765 
    766     int row = 0;
    767     long UTCvalue0=0;
    768     long UTCvalue1=0;
    769     while (fFile->GetNextRow() && row < numRows)
    770     {
    771         rangesIt = ranges.begin();
    772         auto statsIt = statData.begin();
    773         for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
    774         {
    775             double cValue = 0;
    776             for (int i=rangesIt->first; i<rangesIt->second; i++)
    777             {
    778             switch (it->first) {
    779                 case 'L':
    780                         cValue = reinterpret_cast<unsigned char*>(it->second)[i];
    781                         break;
    782                 case 'B':
    783                         cValue = reinterpret_cast<bool*>(it->second)[i];
    784                         break;
    785                 case 'I':
    786                         cValue = reinterpret_cast<int16_t*>(it->second)[i];
    787                         break;
    788                 case 'J':
    789                         cValue = reinterpret_cast<int32_t*>(it->second)[i];
    790                         break;
    791                 case 'K':
    792                         cValue = reinterpret_cast<int64_t*>(it->second)[i];
    793                         break;
    794                 case 'E':
    795                         cValue = reinterpret_cast<float*>(it->second)[i];
    796                         break;
    797                 case 'D':
    798                         cValue = reinterpret_cast<double*>(it->second)[i];
    799                         break;
    800                 default:
    801                     ;
    802             }
    803             if (list.size() == 1 && (list[0] == "UnixTimeUTC" || list[0] == "PCTime"))
    804             {
    805                 if (i==0)
    806                 {
    807                     UTCvalue0 = cValue;
    808                 }
    809                 else
    810                 {
    811                     UTCvalue1 = cValue;
    812                     boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
    813                             boost::posix_time::seconds(UTCvalue0) +  boost::posix_time::microsec(UTCvalue1));
    814 
    815                     Time mjdTime(unixTimeT);
    816                     cValue = mjdTime.Mjd();
    817                     if (!fNoZeroPlease || cValue != 0)
    818                     {
    819                         statsIt->average += cValue;
    820                         if (cValue < statsIt->min)
    821                             statsIt->min = cValue;
    822                         if (cValue > statsIt->max)
    823                             statsIt->max = cValue;
    824                         statsIt->numValues++;
    825                     }
    826                 }
    827             }
    828             else {
    829                 if (!fNoZeroPlease || cValue != 0)
    830                 {
    831                     statsIt->average += cValue;
    832                     if (cValue < statsIt->min)
    833                         statsIt->min = cValue;
    834                     if (cValue > statsIt->max)
    835                         statsIt->max = cValue;
    836                     statsIt->numValues++;
    837                 }
    838             }
    839             }
    840         }
    841         row++;
    842     }
    843     for (auto it = columnsData.begin(); it != columnsData.end(); it++)
    844         delete[] it->second;
    845 */
    846     return true;
    847 }
    848 
    849 /*
    850 void FitsDumper::doTBoundary(const string& filename, int precision, bool tStop)
    851 {
    852 
    853     //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    854      vector<pair<int, int> > ranges;
    855      vector<string> listNamesOnly;
    856 
    857      if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
    858      {
    859          cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
    860          return false;
    861      }
    862 
    863      ofstream out(filename=="-"?"/dev/stdout":filename);
    864      if (!out)
    865      {
    866          cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
    867          return false;
    868      }
    869 
    870      out.precision(precision);
    871 
    872      // Loop over all columns in our list of requested columns
    873      vector<pair<char, char*> > columnsData;
    874      vector<minMaxStruct> statData;
    875      int numRows = fFile->GetInt("NAXIS2");
    876      auto rangesIt = ranges.begin();
    877      for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
    878      {
    879          fits::Table::Column& cCol = fColMap.find(*it)->second;
    880          columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
    881  //        minMaxStuct initData;
    882          statData.push_back(minMaxStruct());
    883          fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
    884      }
    885 
    886      int row = 0;
    887      while (fFile->GetNextRow() && row < numRows)
    888      {
    889          rangesIt = ranges.begin();
    890          auto statsIt = statData.begin();
    891          for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
    892          {
    893              double cValue = 0;
    894              for (int i=rangesIt->first; i<rangesIt->second; i++)
    895              {
    896              switch (it->first) {
    897                  case 'L':
    898                          cValue = reinterpret_cast<bool*>(it->second)[i];
    899                          break;
    900                  case 'B':
    901                          cValue = reinterpret_cast<bool*>(it->second)[i];
    902                          break;
    903                  case 'I':
    904                          cValue = reinterpret_cast<int16_t*>(it->second)[i];
    905                          break;
    906                  case 'J':
    907                          cValue = reinterpret_cast<int32_t*>(it->second)[i];
    908                          break;
    909                  case 'K':
    910                          cValue = reinterpret_cast<int64_t*>(it->second)[i];
    911                          break;
    912                  case 'E':
    913                          cValue = reinterpret_cast<float*>(it->second)[i];
    914                          break;
    915                  case 'D':
    916                          cValue = reinterpret_cast<double*>(it->second)[i];
    917                          break;
    918                  default:
    919                      ;
    920              }
    921              if (!fNoZeroPlease || cValue != 0)
    922              {
    923                  statsIt->average += cValue;
    924                  if (cValue < statsIt->min)
    925                      statsIt->min = cValue;
    926                  if (cValue > statsIt->max)
    927                      statsIt->max = cValue;
    928                  statsIt->numValues++;
    929              }
    930              }
    931          }
    932          row++;
    933      }
    934      for (auto it = columnsData.begin(); it != columnsData.end(); it++)
    935          delete[] it->second;
    936 
    937      //okay. So now I've got ALL the data, loaded.
    938      //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
    939      rangesIt = ranges.begin();
    940      auto statsIt = statData.begin();
    941 
    942      auto nameIt = listNamesOnly.begin();
    943      for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
    944      {
    945          int span = rangesIt->second - rangesIt->first;
    946          cout << *nameIt << ": " << endl;
    947          if (statsIt->numValues != 0)
    948          {
    949              statsIt->average /= statsIt->numValues;
    950              out << "min: " << statsIt->min << endl;
    951              out << "max: " << statsIt->max << endl;
    952              out << "mea: " << statsIt->average << endl;
    953          }
    954          else
    955          {
    956              out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
    957          }
    958 
    959      }
    960      return true;
    961 
    962 }
    963 */
    964 template<typename T>
    965 void displayStats(char* array, int numElems, ofstream& out)
    966 {
    967     if (numElems == 0)
    968     {
    969         out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
    970         return;
    971     }
    972 
    973     T *val = reinterpret_cast<T*>(array);
    974 
    975     sort(val, val+numElems);
    976 
    977     out << "Min: " << double(val[0]) << '\n';
    978     out << "Max: " << double(val[numElems-1]) << '\n';
    979 
    980     if (numElems>2)
    981     {
    982         if (numElems%2 == 0)
    983             out << "Med: " << double(val[(numElems-1)/2] + val[(numElems-1)/2+1])/2 << '\n';
    984         else
    985             out << "Med: " << double(val[numElems/2+1]) << '\n';
    986     }
    987 
    988     double avg = 0;
    989     double rms = 0;
    990     for (int i=0;i<numElems;i++)
    991     {
    992         avg += val[i];
    993         rms += val[i]*val[i];
    994     }
    995 
    996     avg /= numElems;
    997     rms  = sqrt(rms/numElems - avg*avg);
    998 
    999     out << "Avg: " << avg << '\n';
    1000     out << "Rms: " << rms << endl;
    1001 
    1002 }
    1003 int FitsDumper::doStatsPlease(const string &filename, const vector<string>& list, int precision)
    1004 {
    1005 
    1006     //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    1007     vector<pair<int, int> > ranges;
    1008     vector<string> names;
    1009 
    1010     if (!separateColumnsFromRanges(list, ranges, names))
    1011     {
    1012         cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
    1013         return false;
    1014     }
    1015 
    1016     ofstream out(filename=="-"?"/dev/stdout":filename);
    1017     if (!out)
    1018     {
    1019         cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
    1020         return false;
    1021     }
    1022 
    1023     out.precision(precision);
    1024 
    1025     // Loop over all columns in our list of requested columns
    1026     vector<char*> statData;
    1027 
    1028     auto rangesIt = ranges.begin();
    1029     for (auto in=names.begin(); in!=names.end(); in++, rangesIt++)
    1030     {
    1031         cout << *in << endl;
    1032         const fits::Table::Column &cCol = fColMap.find(*in)->second;
    1033         statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*fFile->GetNumRows()]);
    1034     }
    1035     vector<const void*> ptrs;
    1036     for (auto namesIt=names.begin(); namesIt!=names.end(); namesIt++)
    1037     {
    1038         ptrs.push_back(fFile->SetPtrAddress(*namesIt));
    1039     }
    1040 
    1041     while (fFile->GetNextRow())
    1042     {
    1043         const size_t row = fFile->GetRow();
    1044         if (row==fFile->GetNumRows())
    1045             break;
    1046 
    1047         rangesIt = ranges.begin();
    1048         auto statsIt = statData.begin();
    1049         auto ptrsIt = ptrs.begin();
    1050 
    1051         for (auto in=names.begin(); in!=names.end(); in++, rangesIt++, statsIt++, ptrsIt++)
    1052         {
    1053             const void *ptr = *ptrsIt;
    1054 
    1055             const fits::Table::Column &cCol = fColMap.find(*in)->second;
    1056 
    1057             const int span = rangesIt->second - rangesIt->first;
    1058 
    1059             for (int i=rangesIt->first; i<rangesIt->second; i++)
    1060             {
    1061                 switch (cCol.type)
    1062                 {
    1063                 case 'L':
    1064                         reinterpret_cast<bool*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const bool*>(ptr)[i];
    1065                         break;
    1066                 case 'B':
    1067                         reinterpret_cast<char*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const char*>(ptr)[i];
    1068                         break;
    1069                 case 'I':
    1070                         reinterpret_cast<int16_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int16_t*>(ptr)[i];
    1071                         break;
    1072                 case 'J':
    1073                         reinterpret_cast<int32_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int32_t*>(ptr)[i];
    1074                         break;
    1075                 case 'K':
    1076                         reinterpret_cast<int64_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const int64_t*>(ptr)[i];
    1077                         break;
    1078                 case 'E':
    1079                         reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const float*>(ptr)[i];
    1080                         break;
    1081                 case 'D':
    1082                         reinterpret_cast<double*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<const double*>(ptr)[i];
    1083                         break;
    1084                 default:
    1085                     ;
    1086                 }
    1087             }
    1088         }
    1089     }
    1090 
    1091     //okay. So now I've got ALL the data, loaded.
    1092     //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
    1093     rangesIt = ranges.begin();
    1094     auto statsIt = statData.begin();
    1095 
    1096     for (auto it=list.begin(); it!=list.end(); it++, rangesIt++, statsIt++)
    1097     {
    1098         const fits::Table::Column &cCol = fColMap.find(*it)->second;
    1099 
    1100         const int span = rangesIt->second - rangesIt->first;
    1101 
    1102         out << "[" << *it << "]" << endl;
    1103         switch (cCol.type)
    1104         {
    1105             case 'L':
    1106                     displayStats<bool>(*statsIt, fFile->GetNumRows()*span, out);
    1107                     break;
    1108             case 'B':
    1109                     displayStats<char>(*statsIt, fFile->GetNumRows()*span, out);
    1110                     break;
    1111             case 'I':
    1112                     displayStats<int16_t>(*statsIt, fFile->GetNumRows()*span, out);
    1113                     break;
    1114             case 'J':
    1115                     displayStats<int32_t>(*statsIt, fFile->GetNumRows()*span, out);
    1116                     break;
    1117             case 'K':
    1118                     displayStats<int64_t>(*statsIt, fFile->GetNumRows()*span, out);
    1119                     break;
    1120             case 'E':
    1121                     displayStats<float>(*statsIt, fFile->GetNumRows()*span, out);
    1122                     break;
    1123             case 'D':
    1124                     displayStats<double>(*statsIt, fFile->GetNumRows()*span, out);
    1125                     break;
    1126             default:
    1127                 ;
    1128         }
    1129     }
    1130 
    1131     return true;
    1132 }
    1133 #ifdef PLOTTING_PLEASE
    1134 int FitsDumper::doCurvesDisplay( const vector<string> &list, const string& tableName)
    1135 {
    1136     //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    1137      vector<pair<int, int> > ranges;
    1138      vector<string> listNamesOnly;
    1139      if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
    1140      {
    1141          cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
    1142          return false;
    1143      }
    1144      vector<string> curvesNames;
    1145      stringstream str;
    1146      for (auto it=ranges.begin(), jt=listNamesOnly.begin(); it != ranges.end(); it++, jt++)
    1147      {
    1148          for (int i=it->first; i<it->second;i++)
    1149          {
    1150                  str.str("");
    1151                  str << *jt << "[" << i << "]";
    1152                  curvesNames.push_back(str.str());
    1153          }
    1154      }
    1155      char* handle = new char[17];
    1156      sprintf(handle,"FitsDump Display");
    1157 //    Qt::HANDLE h = *handle;//NULL
    1158     int argc = 1;
    1159     char ** argv = &handle;
    1160     QApplication a(argc, argv);
    1161 
    1162 
    1163 
    1164     QwtPlot* plot = new QwtPlot();
    1165     QwtPlotGrid* grid = new QwtPlotGrid;
    1166     grid->enableX(false);
    1167     grid->enableY(true);
    1168     grid->enableXMin(false);
    1169     grid->enableYMin(false);
    1170     grid->setMajPen(QPen(Qt::black, 0, Qt::DotLine));
    1171     grid->attach(plot);
    1172     plot->setAutoReplot(true);
    1173     string title = tableName;
    1174     plot->setAxisScaleDraw(  QwtPlot::xBottom, new TimeScaleDraw());
    1175 
    1176     QWidget window;
    1177     QHBoxLayout* layout = new QHBoxLayout(&window);
    1178     layout->setContentsMargins(0,0,0,0);
    1179     layout->addWidget(plot);
    1180 
    1181     QwtPlotZoomer zoom(plot->canvas());
    1182     zoom.setRubberBandPen(QPen(Qt::gray, 2, Qt::DotLine));
    1183     zoom.setTrackerPen(QPen(Qt::gray));
    1184     int totalSize = 0;
    1185     for (unsigned int i=0;i<list.size();i++)
    1186         totalSize += ranges[i].second - ranges[i].first;
    1187 
    1188     vector<QwtPlotCurve*> curves(totalSize);
    1189     int ii=0;
    1190     for (auto it = curves.begin(), jt=curvesNames.begin(); it != curves.end(); it++, jt++)
    1191     {
    1192         *it = new QwtPlotCurve(jt->c_str());
    1193         switch (ii%6)
    1194         {
    1195         case 0:
    1196             (*it)->setPen(QColor(255,0,0));
    1197             break;
    1198         case 1:
    1199             (*it)->setPen(QColor(0,255,0));
    1200             break;
    1201         case 2:
    1202             (*it)->setPen(QColor(0,0,255));
    1203             break;
    1204         case 3:
    1205             (*it)->setPen(QColor(255,255,0));
    1206             break;
    1207         case 4:
    1208             (*it)->setPen(QColor(0,255,255));
    1209             break;
    1210         case 5:
    1211             (*it)->setPen(QColor(255,0,255));
    1212             break;
    1213         default:
    1214             (*it)->setPen(QColor(0,0,0));
    1215         };
    1216         ii++;
    1217         if (fDotsPlease)
    1218             (*it)->setStyle(QwtPlotCurve::Dots);
    1219         else
    1220             (*it)->setStyle(QwtPlotCurve::Lines);
    1221         (*it)->attach(plot);
    1222     }
    1223     plot->insertLegend(new QwtLegend(), QwtPlot::RightLegend);
    1224 
    1225 
    1226       vector<pair<char, char*> > columnsData;
    1227 
    1228       for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
    1229       {
    1230           fits::Table::Column& cCol = fColMap.find(*it)->second;
    1231           columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
    1232           fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
    1233       }
    1234 
    1235       //add the time column to the given columns
    1236       if (fColMap.find("Time") == fColMap.end() &&
    1237           fColMap.find("UnixTimeUTC") == fColMap.end())
    1238       {
    1239           cerr << "Error: time column could not be found in given table. Aborting" << endl;
    1240           return false;
    1241       }
    1242       const fits::Table::Column& timeCol = (fColMap.find("Time") != fColMap.end()) ? fColMap.find("Time")->second : fColMap.find("UnixTimeUTC")->second;
    1243       bool unixTime = (fColMap.find("Time") == fColMap.end());
    1244       if (unixTime)
    1245           ranges.push_back(make_pair(0,2));
    1246       else
    1247           ranges.push_back(make_pair(0,1));
    1248       columnsData.push_back(make_pair(timeCol.type, new char[timeCol.num*timeCol.size]));
    1249       fFile->SetPtrAddress(unixTime ? "UnixTimeUTC" : "Time", columnsData[columnsData.size()-1].second);
    1250 
    1251 //      stringstream str;
    1252       str.str("");
    1253 
    1254 
    1255       vector<double*> xValues(totalSize);
    1256       double* yValues;
    1257       cout.precision(10);
    1258       str.precision(20);
    1259       for (auto it=xValues.begin(); it!=xValues.end(); it++)
    1260           *it = new double[fFile->GetInt("NAXIS2")];
    1261 
    1262       yValues = new double[fFile->GetInt("NAXIS2")];
    1263 
    1264       cout.precision(3);
    1265       int endIndex = 0;
    1266       int numRows = fFile->GetInt("NAXIS2");
    1267       for (int i=1;i<numRows;i++)
    1268       {
    1269           fFile->GetNextRow();
    1270           cout << "\r" << "Constructing graph " << ((float)(endIndex)/(float)(fFile->GetInt("NAXIS2")))*100.0 << "%";
    1271           endIndex++;
    1272           auto rangesIt = ranges.begin();
    1273           for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
    1274           {
    1275               for (int j=rangesIt->first; j<rangesIt->second; j++)
    1276               {
    1277               switch (it->first) {
    1278                   case 'L':
    1279                           str << reinterpret_cast<bool*>(it->second)[j] << " ";
    1280                           break;
    1281                   case 'B':
    1282                           str << reinterpret_cast<char*>(it->second)[j] << " ";
    1283                           break;
    1284                   case 'I':
    1285                           str << reinterpret_cast<int16_t*>(it->second)[j] << " ";
    1286                           break;
    1287                   case 'J':
    1288                           str << reinterpret_cast<int32_t*>(it->second)[j] << " ";
    1289                           break;
    1290                   case 'K':
    1291                           str << reinterpret_cast<int64_t*>(it->second)[j] << " ";
    1292                           break;
    1293                   case 'E':
    1294                           str << reinterpret_cast<float*>(it->second)[j] << " ";
    1295                           break;
    1296                   case 'D':
    1297                           str << reinterpret_cast<double*>(it->second)[j] << " ";
    1298                           break;
    1299                   default:
    1300                       ;
    1301               }
    1302               }
    1303           }
    1304 
    1305           for (auto it=xValues.begin(); it!= xValues.end(); it++)
    1306           {
    1307               str >> (*it)[i-1];
    1308           }
    1309           if (unixTime)
    1310           {
    1311               long u1, u2;
    1312               str >> u1 >> u2;
    1313 
    1314               boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
    1315                       boost::posix_time::seconds(u1) +  boost::posix_time::microsec(u2));
    1316 
    1317               Time mjdTime(unixTimeT);
    1318               yValues[i-1] = mjdTime.Mjd();
    1319               if (yValues[i-1] < 40587)
    1320                 yValues[i-1] += 40587;
    1321           }
    1322           else
    1323           {
    1324               str >> yValues[i-1];
    1325               if (yValues[i-1] < 40587)
    1326                   yValues[i-1] += 40587;
    1327 
    1328               Time t(yValues[i-1]);
    1329               string time = t.GetAsStr("%H:%M:%S%F");
    1330               while (time[time.size()-1] == '0' && time.size() > 2)
    1331               {
    1332                   time = time.substr(0, time.size()-1);
    1333               }
    1334           }
    1335           if (i==1)
    1336           {
    1337               Time t(yValues[0]);
    1338               title += " - " + t.GetAsStr("%Y-%m-%d");
    1339               plot->setTitle(title.c_str());
    1340           }
    1341       }
    1342       //set the actual data.
    1343       auto jt = xValues.begin();
    1344       for (auto it=curves.begin(); it != curves.end(); it++, jt++)
    1345           (*it)->setRawData(yValues, *jt, endIndex-1);
    1346 
    1347       QStack<QRectF> stack;
    1348       double minX, minY, maxX, maxY;
    1349       minX = minY = 1e10;
    1350       maxX = maxY = -1e10;
    1351       QRectF rect;
    1352       QPointF point;
    1353       for (auto it=curves.begin(); it!= curves.end(); it++)
    1354       {
    1355           rect = (*it)->boundingRect();
    1356           point = rect.bottomRight();
    1357           if (point.x() < minX) minX = point.x();
    1358           if (point.y() < minY) minY = point.y();
    1359           if (point.x() > maxX) maxX = point.x();
    1360           if (point.y() > maxY) maxY = point.y();
    1361           point = rect.topLeft();
    1362           if (point.x() < minX) minX = point.x();
    1363           if (point.y() < minY) minY = point.y();
    1364           if (point.x() > maxX) maxX = point.x();
    1365           if (point.y() > maxY) maxY = point.y();
    1366       }
    1367       QPointF bottomRight(maxX, minY);
    1368       QPointF topLeft(minX, maxY);
    1369       QPointF center((bottomRight+topLeft)/2.f);
    1370       stack.push(QRectF(topLeft + (topLeft-center)*(.5f),bottomRight + (bottomRight-center)*(.5f)));
    1371       zoom.setZoomStack(stack);
    1372 
    1373 //      delete[] fitsBuffer;
    1374       for (auto it = columnsData.begin(); it != columnsData.end(); it++)
    1375           delete[] it->second;
    1376     window.resize(600, 400);
    1377     window.show();
    1378 
    1379     a.exec();
    1380 
    1381 
    1382     for (auto it = curves.begin(); it != curves.end(); it++)
    1383     {
    1384        (*it)->detach();
    1385         delete *it;
    1386     }
    1387     grid->detach();
    1388     for (auto it = xValues.begin(); it != xValues.end(); it++)
    1389         delete[] *it;
    1390 
    1391     delete[] yValues;
    1392 
    1393     delete[] handle;
    1394 
    1395     return 0;
    1396 }
    1397 #endif
    1398600
    1399601void SetupConfiguration(Configuration& conf)
    … …  
    1406608#endif
    1407609                                                  , "Name of FITS file")
    1408         ("tablename,t", var<string>("DATA")
    1409 #if BOOST_VERSION >= 104200
    1410          ->required()
    1411 #endif
    1412                                                   , "Name of input table")
    1413610        ("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")
    1414611        ("outfile,o",   var<string>("/dev/stdout"), "Name of output file (-:/dev/stdout)")
    … …  
    1419616        ("minmax,m",    po_switch(),                "Calculates min and max of data")
    1420617        ("nozero,z",    po_switch(),                "skip 0 values for stats")
    1421         ("tstart,a",    po_switch(),                "Give the mjdStart from reading the file data")
    1422         ("tstop,b",     po_switch(),                "Give the mjdStop from reading the file data")
    1423618        ("force",       po_switch(),                "Force reading the fits file even if END key is missing")
    1424 #ifdef PLOTTING_PLEASE
    1425         ("graph,g",     po_switch(),                "Plot the columns instead of dumping them")
    1426         ("dots,d",      po_switch(),                "Plot using dots instead of lines")
    1427 #endif
    1428619        ;
    1429620
    1430621    po::positional_options_description p;
    1431     p.add("fitsfile",   1); // The first positional options
    1432     p.add("col",       -1); // All others
     622    p.add("fitsfile",  1); // The first positional options
     623    p.add("col",      -1); // All others
    1433624
    1434625    conf.AddOptions(configs);
    … …  
    1445636        return -1;
    1446637
    1447     FitsDumper loader;
    1448     return loader.ExecConfig(conf);
    1449 }
     638    if (!conf.Has("fitsfile"))
     639    {
     640        cerr << "Filename required." << endl;
     641        return -1;
     642    }
     643
     644    FitsDumper loader(conf.Get<string>("fitsfile"));
     645    if (!loader)
     646    {
     647        cerr << "ERROR - Opening " << conf.Get<string>("fitsfile");
     648        cerr << " failed: " << strerror(errno) << endl;
     649        return -1;
     650    }
     651
     652    return loader.Exec(conf);
     653}
Note: See TracChangeset for help on using the changeset viewer.