//****************************************************************
/** @class FitsDumper

  @brief Dumps contents of fits tables to stdout or a file

 */
 //****************************************************************
#include "Configuration.h"

#include <map>
#include <fstream>

#include "externals/fits.h"

//#define PLOTTING_PLEASE

#ifdef PLOTTING_PLEASE
#include <QPen>
#include <QtGui>
#include <QApplication>

#include <qwt_plot.h>
#include <qwt_plot_grid.h>
#include <qwt_plot_curve.h>
#include <qwt_plot_zoomer.h>
#include <qwt_legend.h>
#include <qwt_scale_draw.h>

#endif
#include "Time.h"

using namespace std;


#ifdef PLOTTING_PLEASE
class TimeScaleDraw: public QwtScaleDraw
{
public:
    virtual QwtText label(double v) const
    {
        Time t(v);
        string time = t.GetAsStr("%H:%M:%S%F");
        while (time[time.size()-1] == '0' && time.size() > 2)
        {
            time = time.substr(0, time.size()-1);
        }
        return QwtText(time.c_str());
    }
};
#endif

class FitsDumper
{
public:
    FitsDumper();
    ~FitsDumper();

private:
    fits* fFile;
    bool fDotsPlease;
    bool fNoZeroPlease;
    string fFilename;

    fits::Table::Columns fColMap;
    fits::Table::Keys fKeyMap;

    // Convert CCfits::ValueType into a human readable string
    string ValueTypeToStr(char type) const;

    // Convert CCfits::ValueType into a number of associated bytes
    int    ValueTypeToSize(char type) const;

    /// Calculate the buffer size required to read a row of the fits table, as well as the offsets to each column
//    vector<int> CalculateOffsets() const;

    template<class T>
        T PtrToValue(const unsigned char* &ptr) const;
//    template<class T>
//        double PtrToDouble(const unsigned char *ptr) const;
//    double PtrToDouble(const unsigned char *ptr, CCfits::ValueType type) const;

    /// Write a single row of the selected data
 //   int  WriteRow(ostream &, const vector<MyColumn*> &, const vector<int> &, unsigned char *, const vector<pair<int, int> >&) const;

    bool OpenFile(const string &);        /// Open a file
    bool OpenTable(const string &);       /// Open a table

    /// Lists all columns of an open file
    void List();                          
    void ListHeader(const string& filename);
    void ListKeywords(ostream &);

    bool separateColumnsFromRanges(const vector<string>& list,
                                   vector<pair<int, int> >& ranges,
                                   vector<string>& listNamesOnly);
    /// Perform the dumping, based on the current dump list
    bool Dump(const string &, const vector<string> &list, int);
    ///Display the selected columns values VS time
#ifdef PLOTTING_PLEASE
    int doCurvesDisplay( const vector<string> &list, const string& tableName);
#endif
    int doMinMaxPlease(const string& filename, const vector<string>& list, int precision);
    int doStatsPlease(const string &filename, const vector<string>& list, int precision);
//    void doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), true);
    //    bool Plot(const vector<string> &list);

public:
    ///Configures the fitsLoader from the config file and/or command arguments.
    int ExecConfig(Configuration& conf);
};

// --------------------------------------------------------------------------
//
//! Constructor
//! @param out
//!        the ostream where to redirect the outputs
//
FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false), fNoZeroPlease(false)
{
}

// --------------------------------------------------------------------------
//
//! Destructor
//
FitsDumper::~FitsDumper()
{
    if (fFile)
        delete fFile;
}


string FitsDumper::ValueTypeToStr(char type) const
{
    switch (type)
    {
        case 'L': return "bool(8)";
        case 'B': return "byte(8)";
        case 'I': return "short(16)";
        case 'J': return "int(32)";
        case 'K': return "int(64)";
        case 'E': return "float(32)";
        case 'D': return "double(64)";
    default:
        return "unknwown";
    }
}

int FitsDumper::ValueTypeToSize(char type) const
{
    switch (type)
    {
        case 'L': return sizeof(uint8_t);
        case 'B': return sizeof(int8_t);
        case 'I': return sizeof(int16_t);
        case 'J': return sizeof(int32_t);
        case 'K': return sizeof(int64_t);
        case 'E': return sizeof(float);
        case 'D': return sizeof(double);
    default:
        return 0;
    }
}

template<class T>
T FitsDumper::PtrToValue(const unsigned char* &ptr) const
{
    T t;
    reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast<unsigned char*>(&t));
    ptr += sizeof(T);

    return t;
}
// --------------------------------------------------------------------------
//
//! Loads the fits file based on the current parameters
//
bool FitsDumper::OpenFile(const string &filename)
{
    if (fFile)
    {
        fFile->close();
        delete fFile;
    }

    try {
        fFile = new fits(filename);
    }
    catch (std::runtime_error e)
    {
        cout << "Something went wrong while trying to open " << filename;
        cout << ": " << e.what() << " Aborting dump." << endl;
        return false;
    }
    fFilename = filename;

    const fits::Table::Columns& tCols = fFile->GetColumns();

    for (auto it=tCols.begin(); it != tCols.end(); it++)
        fColMap.insert(*it);

    const fits::Table::Keys& tkeys = fFile->GetKeys();

    for (auto it=tkeys.begin(); it != tkeys.end(); it++)
        fKeyMap.insert(*it);

    return true;
}

bool FitsDumper::OpenTable(const string &)
{
    if (!fFile)
    {
        cerr << "No file open." << endl;
        return false;
    }


    return true;
}


void FitsDumper::List()
{
    if (!fFile)
    {
        cerr << "No file open." << endl;
        return;
    }

    cout << "\nFile: " << fFilename << "\n";

    cout << " " << fKeyMap.find("EXTNAME")->second.value << " [";
    cout << fKeyMap.find("NAXIS2")->second.value << "]\n";

    for (auto it = fColMap.begin(); it != fColMap.end(); it++)
    {
        cout << "   " << it->first << "[" << it->second.num << "] (" << it->second.unit << ":" << ValueTypeToStr(it->second.type) << ") ";
        for (auto jt = fKeyMap.begin(); jt != fKeyMap.end(); jt++)
            if (jt->second.value == it->first)
                cout << jt->second.comment << endl;
    }

    cout << endl;
    cout << flush;
}

void FitsDumper::ListKeywords(ostream &out)
{
    for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++) {
        out << "## " << setw(8) << it->first << " = " << setw(10);
        out << "'" << it->second.value << "'" << " / " << it->second.comment << endl;
    }
}

void FitsDumper::ListHeader(const string& filename)
{
    if (!fFile)
    {
        cerr << "No table open." << endl;
        return;
    }
    ofstream out(filename=="-"?"/dev/stdout":filename);
    if (!out)
    {
        cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
        return;
    }

    out << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n";
    if (fKeyMap.find("COMMENT") != fKeyMap.end())
        out << "Comment: \t" << fKeyMap.find("COMMENT")->second.value << "\n";

    ListKeywords(out);
    out << endl;

}

bool FitsDumper::separateColumnsFromRanges(const vector<string>& list,
                               vector<pair<int, int> >& ranges,
                               vector<string>& listNamesOnly)
{
    for (vector<string>::const_iterator it=list.begin(); it!=list.end(); it++)
    {
        string columnNameOnly = *it;
        unsigned long bracketIndex0 = columnNameOnly.find_first_of('[');
        unsigned long bracketIndex1 = columnNameOnly.find_first_of(']');
        unsigned long colonIndex = columnNameOnly.find_first_of(':');
        int columnStart = -1;
        int columnEnd = -1;
        if (bracketIndex0 != string::npos)
        {//there is a range given. Extract the range
            if (colonIndex != string::npos)
            {//we have a range here
                columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, colonIndex-(bracketIndex0+1)).c_str());
                columnEnd = atoi(columnNameOnly.substr(colonIndex+1, bracketIndex1-(colonIndex+1)).c_str());
                columnEnd++;
            }
            else
            {//only a single index there
                columnStart = atoi(columnNameOnly.substr(bracketIndex0+1, bracketIndex1 - (bracketIndex0+1)).c_str());
                columnEnd = columnStart+1;
            }
            columnNameOnly = columnNameOnly.substr(0, bracketIndex0);
        }

        if (fColMap.find(columnNameOnly) == fColMap.end())
        {
            cerr << "ERROR - Column '" << columnNameOnly << "' not found in table." << endl;
            return false;
        }
        fits::Table::Column& cCol = fColMap.find(columnNameOnly)->second;
        if (bracketIndex0 == string::npos)
        {//no range given: use the full range
            ranges.push_back(make_pair(0, cCol.num));
            columnStart = 0;
            columnEnd = 1;
        }
        else
        {//use the range extracted earlier
            if (columnStart < 0)
            {
                cerr << "ERROR - Start range for column " << columnNameOnly << " is less than zero (" << columnStart << "). Aborting" << endl;
                return false;
            }
            if (columnEnd>1 && columnEnd > (int)(cCol.num))
            {
                cerr << "ERROR - End range for column " << columnNameOnly << " is greater than the last element (" << cCol.num << " vs " << columnEnd << "). Aborting" << endl;
                return false;
            }
            ranges.push_back(make_pair(columnStart, columnEnd));
        }
//        cout << "Will be exporting from " << columnStart << " to " << columnEnd-1 << " for column " << columnNameOnly << endl;
        listNamesOnly.push_back(columnNameOnly);
    }
    return true;
}
// --------------------------------------------------------------------------
//
//! Perform the actual dump, based on the current parameters
//
bool FitsDumper::Dump(const string &filename, const vector<string> &list, int precision)
{

    //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    vector<pair<int, int> > ranges;
    vector<string> listNamesOnly;

    if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
    {
        cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
        return false;
    }

    ofstream out(filename=="-"?"/dev/stdout":filename);
    if (!out)
    {
        cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
        return false;
    }

    out.precision(precision);

    out << "## --------------------------------------------------------------------------\n";
    if (filename!="-")
        out << "## File:    \t" << filename << '\n';
    out << "## Table:   \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
        out << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';

    out << "## NumRows: \t" << fFile->GetInt("NAXIS2") << '\n';
    out << "## --------------------------------------------------------------------------\n";
    ListKeywords(out);
    out << "## --------------------------------------------------------------------------\n";
    out << "#\n";

    auto rangesIt = ranges.begin();
    for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
    {
        const fits::Table::Column& col = fColMap[*it];
//        const MyColumn *col = static_cast<MyColumn*>(fTable->column()[*it]);
        if (rangesIt->first != 0 || rangesIt->second != (int)(col.num))
        {
            out << "#";
            for (int i=rangesIt->first; i<rangesIt->second; i++)
                out << " " << *it << "[" << i << "]";
            out << ": " << col.unit;
        }
        else
            out << "# " << *it << "[" << col.num << "]: " << col.unit;

//        FIXME: retrive the column comment
//        if (!col->comment().empty())
//            out << " (" <<col->comment() << ")";
        out << '\n';
    }
    out << "#" << endl;


    // Loop over all columns in our list of requested columns
    vector<pair<char, char*> > columnsData;

    for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
    {
        fits::Table::Column& cCol = fColMap.find(*it)->second;
        columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
        fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
    }
    int numRows = fFile->GetInt("NAXIS2");
    int row = 0;
    while (fFile->GetNextRow() && row < numRows)
    {
        row++;
        rangesIt = ranges.begin();
        for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
        {
            for (int i=rangesIt->first; i<rangesIt->second; i++)
            {
            switch (it->first) {
                case 'L':
                        out << reinterpret_cast<bool*>(it->second)[i] << " ";
                        break;
                case 'B':
                        out << reinterpret_cast<char*>(it->second)[i] << " ";
                        break;
                case 'I':
                        out << reinterpret_cast<int16_t*>(it->second)[i] << " ";
                        break;
                case 'J':
                        out << reinterpret_cast<int32_t*>(it->second)[i] << " ";
                        break;
                case 'K':
                        out << reinterpret_cast<int64_t*>(it->second)[i] << " ";
                        break;
                case 'E':
                        out << reinterpret_cast<float*>(it->second)[i] << " ";
                        break;
                case 'D':
                        out << reinterpret_cast<double*>(it->second)[i] << " ";
                        break;
                default:
                    ;
            }
            }
        }
    }
    for (auto it = columnsData.begin(); it != columnsData.end(); it++)
        delete[] it->second;
    return true;
}

// --------------------------------------------------------------------------
//
//! Retrieves the configuration parameters
//! @param conf
//!             the configuration object
//
int FitsDumper::ExecConfig(Configuration& conf)
{
    if (conf.Has("fitsfile"))
    {
        if (!OpenFile(conf.Get<string>("fitsfile")))
            return -1;
    }
#ifdef PLOTTING_PLEASE
    if (conf.Get<bool>("stat") && conf.Get<bool>("graph"))
    {
        cout << "Invalid conbination of options: cannot graph stats. Aborting" << endl;
        return -1;
    }
#endif
    if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
    {
        cout << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl;
        return -1;
    }
#ifdef PLOTTING_PLEASE
    if (conf.Get<bool>("minmax") && conf.Get<bool>("graph"))
    {
        cout << "Invalid combination of options: cannot graph minmax. Aborting" << endl;
        return -1;
    }
#endif
    if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
    {
        cout << "Invalid combination of options: nozero only works with minmax. Aborting" << endl;
        return -1;
    }

    if (conf.Get<bool>("nozero"))
    {
        fNoZeroPlease = true;
    }

    if (conf.Has("tablename"))
    {
        if (!OpenTable(conf.Get<string>("tablename")))
            return -1;
    }

    if (conf.Get<bool>("list"))
        List();
    /*
    if (conf.Get<bool>("tstart"))
    {
        doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), true);
    }
    if (conf.Get<bool>("tstop"))
    {
        doTBoundary(conf.Get<string>("outfile"), conf.Get<int>("precision"), false);
    }
    */
    if (conf.Get<bool>("minmax"))
    {
        if (!conf.Has("col"))
        {
            cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
            return 0;
        }
        doMinMaxPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
        return 0;
    }

    if (conf.Get<bool>("stat"))
    {
        if (!conf.Has("col"))
        {
            cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
            return 0;
        }
        doStatsPlease(conf.Get<string>("outfile"), conf.Get<vector<string>>("col"), conf.Get<int>("precision"));
        return 0;
    }

#ifdef PLOTTING_PLEASE
    if (conf.Get<bool>("graph"))
    {
        if (!conf.Has("col"))
        {
            cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
            return 0;
        }
        if (conf.Get<bool>("dots"))
            fDotsPlease = true;
        doCurvesDisplay(conf.Get<vector<string>>("col"),
                        conf.Get<string>("tablename"));
        return 1;
    }
#endif

    if (conf.Get<bool>("header"))
        ListHeader(conf.Get<string>("outfile"));

    if (conf.Get<bool>("header") || conf.Get<bool>("list"))
        return 1;

    if (conf.Has("outfile"))
    {
        if (!conf.Has("col"))
        {
            cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl;
            return 0;
        }
        if (!Dump(conf.Get<string>("outfile"),
                  conf.Get<vector<string>>("col"),
                  conf.Get<int>("precision")))
            return -1;
    }


    return 0;
}

void PrintUsage()
{
    cout <<
        "fitsdump is a tool to dump data from a FITS table as ascii.\n"
        "\n"
        "Usage: fitsdump [OPTIONS] fitsfile col col ... \n"
        "  or:  fitsdump [OPTIONS]\n";
    cout << endl;
}

void PrintHelp()
{
    // 
}

struct minMaxStruct {
    double min;
    double max;
    double average;
    long numValues;
    minMaxStruct() : min(1e10), max(-1e10), average(0), numValues(0) {}
};

int FitsDumper::doMinMaxPlease(const string& filename, const vector<string>& list, int precision)
{
    //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    vector<pair<int, int> > ranges;
    vector<string> listNamesOnly;

    if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
    {
        cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
        return false;
    }

    ofstream out(filename=="-"?"/dev/stdout":filename);
    if (!out)
    {
        cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
        return false;
    }

    out.precision(precision);

    // Loop over all columns in our list of requested columns
    vector<pair<char, char*> > columnsData;
    vector<minMaxStruct> statData;
    int numRows = fFile->GetInt("NAXIS2");
    auto rangesIt = ranges.begin();
    for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
    {
        fits::Table::Column& cCol = fColMap.find(*it)->second;
        columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
//        minMaxStuct initData;
        statData.push_back(minMaxStruct());
        fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
    }

    int row = 0;
    long UTCvalue0=0;
    long UTCvalue1=0;
    while (fFile->GetNextRow() && row < numRows)
    {
        rangesIt = ranges.begin();
        auto statsIt = statData.begin();
        for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
        {
            double cValue = 0;
            for (int i=rangesIt->first; i<rangesIt->second; i++)
            {
            switch (it->first) {
                case 'L':
                        cValue = reinterpret_cast<bool*>(it->second)[i];
                        break;
                case 'B':
                        cValue = reinterpret_cast<bool*>(it->second)[i];
                        break;
                case 'I':
                        cValue = reinterpret_cast<int16_t*>(it->second)[i];
                        break;
                case 'J':
                        cValue = reinterpret_cast<int32_t*>(it->second)[i];
                        break;
                case 'K':
                        cValue = reinterpret_cast<int64_t*>(it->second)[i];
                        break;
                case 'E':
                        cValue = reinterpret_cast<float*>(it->second)[i];
                        break;
                case 'D':
                        cValue = reinterpret_cast<double*>(it->second)[i];
                        break;
                default:
                    ;
            }
            if (list.size() == 1 && (list[0] == "UnixTimeUTC" || list[0] == "PCTime"))
            {
                if (i==0)
                {
                    UTCvalue0 = cValue;
                }
                else
                {
                    UTCvalue1 = cValue;
                    boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
                            boost::posix_time::seconds(UTCvalue0) +  boost::posix_time::microsec(UTCvalue1));

                    Time mjdTime(unixTimeT);
                    cValue = mjdTime.Mjd();
                    if (!fNoZeroPlease || cValue != 0)
                    {
                        statsIt->average += cValue;
                        if (cValue < statsIt->min)
                            statsIt->min = cValue;
                        if (cValue > statsIt->max)
                            statsIt->max = cValue;
                        statsIt->numValues++;
                    }
                }
            }
            else {
                if (!fNoZeroPlease || cValue != 0)
                {
                    statsIt->average += cValue;
                    if (cValue < statsIt->min)
                        statsIt->min = cValue;
                    if (cValue > statsIt->max)
                        statsIt->max = cValue;
                    statsIt->numValues++;
                }
            }
            }
        }
        row++;
    }
    for (auto it = columnsData.begin(); it != columnsData.end(); it++)
        delete[] it->second;

    //okay. So now I've got ALL the data, loaded.
    //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
    rangesIt = ranges.begin();
    auto statsIt = statData.begin();

    auto nameIt = listNamesOnly.begin();
    for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
    {
//        int span = rangesIt->second - rangesIt->first;
        cout << *nameIt << ": " << endl;
        if (statsIt->numValues != 0)
        {
            statsIt->average /= statsIt->numValues;
            out << "min: " << statsIt->min << endl;
            out << "max: " << statsIt->max << endl;
            out << "mea: " << statsIt->average << endl;
        }
        else
        {
            out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
        }

    }
    return true;
}

/*
void FitsDumper::doTBoundary(const string& filename, int precision, bool tStop)
{

    //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
     vector<pair<int, int> > ranges;
     vector<string> listNamesOnly;

     if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
     {
         cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
         return false;
     }

     ofstream out(filename=="-"?"/dev/stdout":filename);
     if (!out)
     {
         cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
         return false;
     }

     out.precision(precision);

     // Loop over all columns in our list of requested columns
     vector<pair<char, char*> > columnsData;
     vector<minMaxStruct> statData;
     int numRows = fFile->GetInt("NAXIS2");
     auto rangesIt = ranges.begin();
     for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
     {
         fits::Table::Column& cCol = fColMap.find(*it)->second;
         columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
 //        minMaxStuct initData;
         statData.push_back(minMaxStruct());
         fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
     }

     int row = 0;
     while (fFile->GetNextRow() && row < numRows)
     {
         rangesIt = ranges.begin();
         auto statsIt = statData.begin();
         for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
         {
             double cValue = 0;
             for (int i=rangesIt->first; i<rangesIt->second; i++)
             {
             switch (it->first) {
                 case 'L':
                         cValue = reinterpret_cast<bool*>(it->second)[i];
                         break;
                 case 'B':
                         cValue = reinterpret_cast<bool*>(it->second)[i];
                         break;
                 case 'I':
                         cValue = reinterpret_cast<int16_t*>(it->second)[i];
                         break;
                 case 'J':
                         cValue = reinterpret_cast<int32_t*>(it->second)[i];
                         break;
                 case 'K':
                         cValue = reinterpret_cast<int64_t*>(it->second)[i];
                         break;
                 case 'E':
                         cValue = reinterpret_cast<float*>(it->second)[i];
                         break;
                 case 'D':
                         cValue = reinterpret_cast<double*>(it->second)[i];
                         break;
                 default:
                     ;
             }
             if (!fNoZeroPlease || cValue != 0)
             {
                 statsIt->average += cValue;
                 if (cValue < statsIt->min)
                     statsIt->min = cValue;
                 if (cValue > statsIt->max)
                     statsIt->max = cValue;
                 statsIt->numValues++;
             }
             }
         }
         row++;
     }
     for (auto it = columnsData.begin(); it != columnsData.end(); it++)
         delete[] it->second;

     //okay. So now I've got ALL the data, loaded.
     //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
     rangesIt = ranges.begin();
     auto statsIt = statData.begin();

     auto nameIt = listNamesOnly.begin();
     for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
     {
         int span = rangesIt->second - rangesIt->first;
         cout << *nameIt << ": " << endl;
         if (statsIt->numValues != 0)
         {
             statsIt->average /= statsIt->numValues;
             out << "min: " << statsIt->min << endl;
             out << "max: " << statsIt->max << endl;
             out << "mea: " << statsIt->average << endl;
         }
         else
         {
             out << "min: 0" << endl << "max: 0" << endl << "mea: " << endl;
         }

     }
     return true;

}
*/
template<typename T>
void displayStats(char* array, int numElems, ofstream& out)
{
	if (numElems == 0)
	{
		out << "Min: 0" << endl;
		out << "Max: 0" << endl;
		out << "Med: 0" << endl;
		out << "Mea: 0" << endl;
		return;
	}
    sort(reinterpret_cast<T*>(array), &reinterpret_cast<T*>(array)[numElems]);
    out << "Min: " << reinterpret_cast<T*>(array)[0] << endl;
    out << "Max: " << reinterpret_cast<T*>(array)[numElems-1] << endl;
    if (numElems%2 == 0)
        out << "Med: " << (reinterpret_cast<T*>(array)[numElems/2] + reinterpret_cast<int16_t*>(array)[numElems/2+1])/2.f << endl;
    else
        out << "Med: " << reinterpret_cast<T*>(array)[numElems/2+1] << endl;

    double average = 0;
    for (int i=0;i<numElems;i++)
        average += reinterpret_cast<T*>(array)[i];
    out << "Mea: " << average/(double)numElems << endl;

}
int FitsDumper::doStatsPlease(const string &filename, const vector<string>& list, int precision)
{

    //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
    vector<pair<int, int> > ranges;
    vector<string> listNamesOnly;

    if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
    {
        cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
        return false;
    }

    ofstream out(filename=="-"?"/dev/stdout":filename);
    if (!out)
    {
        cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
        return false;
    }

    out.precision(precision);

    // Loop over all columns in our list of requested columns
    vector<pair<char, char*> > columnsData;
    vector<char*> statData;
    int numRows = fFile->GetInt("NAXIS2");
    auto rangesIt = ranges.begin();
    for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++)
    {
        fits::Table::Column& cCol = fColMap.find(*it)->second;
        columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
        statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*numRows]);
        fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
    }

    int row = 0;
    while (fFile->GetNextRow() && row < numRows)
    {
        rangesIt = ranges.begin();
        auto statsIt = statData.begin();
        for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++)
        {
            int span = rangesIt->second - rangesIt->first;
            for (int i=rangesIt->first; i<rangesIt->second; i++)
            {
            switch (it->first) {
                case 'L':
                        reinterpret_cast<bool*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<bool*>(it->second)[i];
                        break;
                case 'B':
                        reinterpret_cast<char*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<char*>(it->second)[i];
                        break;
                case 'I':
                        reinterpret_cast<int16_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int16_t*>(it->second)[i];
                        break;
                case 'J':
                        reinterpret_cast<int32_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int32_t*>(it->second)[i];
                        break;
                case 'K':
                        reinterpret_cast<int64_t*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<int64_t*>(it->second)[i];
                        break;
                case 'E':
                        reinterpret_cast<float*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<float*>(it->second)[i];
                        break;
                case 'D':
                        reinterpret_cast<double*>(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast<double*>(it->second)[i];
                        break;
                default:
                    ;
            }
            }
        }
        row++;
    }
    for (auto it = columnsData.begin(); it != columnsData.end(); it++)
        delete[] it->second;

    //okay. So now I've got ALL the data, loaded.
    //let's do the summing and averaging in a safe way (i.e. avoid overflow of variables as much as possible)
    rangesIt = ranges.begin();
    auto statsIt = statData.begin();

    auto nameIt = list.begin();
    for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++, statsIt++, nameIt++)
    {
        int span = rangesIt->second - rangesIt->first;
        out << *nameIt << ": " << endl;
        switch (it->first) {
            case 'L':
                    displayStats<bool>(*statsIt, numRows*span, out);
                    break;
            case 'B':
                    displayStats<char>(*statsIt, numRows*span, out);
                    break;
            case 'I':
                    displayStats<int16_t>(*statsIt, numRows*span, out);
                    break;
            case 'J':
                    displayStats<int32_t>(*statsIt, numRows*span, out);
                    break;
            case 'K':
                    displayStats<int64_t>(*statsIt, numRows*span, out);
                    break;
            case 'E':
                    displayStats<float>(*statsIt, numRows*span, out);
                    break;
            case 'D':
                    displayStats<double>(*statsIt, numRows*span, out);
                    break;
            default:
                ;
        }
    }
    return true;
}
#ifdef PLOTTING_PLEASE
int FitsDumper::doCurvesDisplay( const vector<string> &list, const string& tableName)
{
    //first of all, let's separate the columns from their ranges and check that the requested columns are indeed part of the file
     vector<pair<int, int> > ranges;
     vector<string> listNamesOnly;
     if (!separateColumnsFromRanges(list, ranges, listNamesOnly))
     {
         cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl;
         return false;
     }
     vector<string> curvesNames;
     stringstream str;
     for (auto it=ranges.begin(), jt=listNamesOnly.begin(); it != ranges.end(); it++, jt++)
     {
         for (int i=it->first; i<it->second;i++)
         {
                 str.str("");
                 str << *jt << "[" << i << "]";
                 curvesNames.push_back(str.str());
         }
     }
     char* handle = new char[17];
     sprintf(handle,"FitsDump Display");
//    Qt::HANDLE h = *handle;//NULL
    int argc = 1;
    char ** argv = &handle;
    QApplication a(argc, argv);



    QwtPlot* plot = new QwtPlot();
    QwtPlotGrid* grid = new QwtPlotGrid;
    grid->enableX(false);
    grid->enableY(true);
    grid->enableXMin(false);
    grid->enableYMin(false);
    grid->setMajPen(QPen(Qt::black, 0, Qt::DotLine));
    grid->attach(plot);
    plot->setAutoReplot(true);
    string title = tableName;
    plot->setAxisScaleDraw(  QwtPlot::xBottom, new TimeScaleDraw());

    QWidget window;
    QHBoxLayout* layout = new QHBoxLayout(&window);
    layout->setContentsMargins(0,0,0,0);
    layout->addWidget(plot);

    QwtPlotZoomer zoom(plot->canvas());
    zoom.setRubberBandPen(QPen(Qt::gray, 2, Qt::DotLine));
    zoom.setTrackerPen(QPen(Qt::gray));
    int totalSize = 0;
    for (unsigned int i=0;i<list.size();i++)
        totalSize += ranges[i].second - ranges[i].first;

    vector<QwtPlotCurve*> curves(totalSize);
    int ii=0;
    for (auto it = curves.begin(), jt=curvesNames.begin(); it != curves.end(); it++, jt++)
    {
        *it = new QwtPlotCurve(jt->c_str());
        switch (ii%6)
        {
        case 0:
            (*it)->setPen(QColor(255,0,0));
            break;
        case 1:
            (*it)->setPen(QColor(0,255,0));
            break;
        case 2:
            (*it)->setPen(QColor(0,0,255));
            break;
        case 3:
            (*it)->setPen(QColor(255,255,0));
            break;
        case 4:
            (*it)->setPen(QColor(0,255,255));
            break;
        case 5:
            (*it)->setPen(QColor(255,0,255));
            break;
        default:
            (*it)->setPen(QColor(0,0,0));
        };
        ii++;
        if (fDotsPlease)
            (*it)->setStyle(QwtPlotCurve::Dots);
        else
            (*it)->setStyle(QwtPlotCurve::Lines);
        (*it)->attach(plot);
    }
    plot->insertLegend(new QwtLegend(), QwtPlot::RightLegend);


      vector<pair<char, char*> > columnsData;

      for (vector<string>::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++)
      {
          fits::Table::Column& cCol = fColMap.find(*it)->second;
          columnsData.push_back(make_pair(cCol.type, new char[cCol.num*cCol.size]));
          fFile->SetPtrAddress(*it, columnsData[columnsData.size()-1].second);
      }

      //add the time column to the given columns
      if (fColMap.find("Time") == fColMap.end() && 
          fColMap.find("UnixTimeUTC") == fColMap.end())
      {
          cerr << "Error: time column could not be found in given table. Aborting" << endl;
          return false;
      }
      const fits::Table::Column& timeCol = (fColMap.find("Time") != fColMap.end()) ? fColMap.find("Time")->second : fColMap.find("UnixTimeUTC")->second;
      bool unixTime = (fColMap.find("Time") == fColMap.end());
      if (unixTime)
          ranges.push_back(make_pair(0,2));
      else
          ranges.push_back(make_pair(0,1));
      columnsData.push_back(make_pair(timeCol.type, new char[timeCol.num*timeCol.size]));
      fFile->SetPtrAddress(unixTime ? "UnixTimeUTC" : "Time", columnsData[columnsData.size()-1].second);

//      stringstream str;
      str.str("");


      vector<double*> xValues(totalSize);
      double* yValues;
      cout.precision(10);
      str.precision(20);
      for (auto it=xValues.begin(); it!=xValues.end(); it++)
          *it = new double[fFile->GetInt("NAXIS2")];

      yValues = new double[fFile->GetInt("NAXIS2")];

      cout.precision(3);
      int endIndex = 0;
      int numRows = fFile->GetInt("NAXIS2");
      for (int i=1;i<numRows;i++)
      {
          fFile->GetNextRow();
          cout << "\r" << "Constructing graph " << ((float)(endIndex)/(float)(fFile->GetInt("NAXIS2")))*100.0 << "%";
          endIndex++;
          auto rangesIt = ranges.begin();
          for (auto it=columnsData.begin(); it != columnsData.end(); it++, rangesIt++)
          {
              for (int j=rangesIt->first; j<rangesIt->second; j++)
              {
              switch (it->first) {
                  case 'L':
                          str << reinterpret_cast<bool*>(it->second)[j] << " ";
                          break;
                  case 'B':
                          str << reinterpret_cast<char*>(it->second)[j] << " ";
                          break;
                  case 'I':
                          str << reinterpret_cast<int16_t*>(it->second)[j] << " ";
                          break;
                  case 'J':
                          str << reinterpret_cast<int32_t*>(it->second)[j] << " ";
                          break;
                  case 'K':
                          str << reinterpret_cast<int64_t*>(it->second)[j] << " ";
                          break;
                  case 'E':
                          str << reinterpret_cast<float*>(it->second)[j] << " ";
                          break;
                  case 'D':
                          str << reinterpret_cast<double*>(it->second)[j] << " ";
                          break;
                  default:
                      ;
              }
              }
          }

          for (auto it=xValues.begin(); it!= xValues.end(); it++)
          {
              str >> (*it)[i-1];
          }
          if (unixTime)
          {
              long u1, u2;
              str >> u1 >> u2;

              boost::posix_time::ptime unixTimeT( boost::gregorian::date(1970, boost::gregorian::Jan, 1),
                      boost::posix_time::seconds(u1) +  boost::posix_time::microsec(u2));

              Time mjdTime(unixTimeT);
              yValues[i-1] = mjdTime.Mjd();
	      if (yValues[i-1] < 40587)
	      	yValues[i-1] += 40587;
          }
          else
          {
              str >> yValues[i-1];
              if (yValues[i-1] < 40587)
                  yValues[i-1] += 40587;

              Time t(yValues[i-1]);
              string time = t.GetAsStr("%H:%M:%S%F");
              while (time[time.size()-1] == '0' && time.size() > 2)
              {
                  time = time.substr(0, time.size()-1);
              }
          }
          if (i==1)
          {
              Time t(yValues[0]);
              title += " - " + t.GetAsStr("%Y-%m-%d");
              plot->setTitle(title.c_str());
          }
      }
      //set the actual data.
      auto jt = xValues.begin();
      for (auto it=curves.begin(); it != curves.end(); it++, jt++)
          (*it)->setRawData(yValues, *jt, endIndex-1);

      QStack<QRectF> stack;
      double minX, minY, maxX, maxY;
      minX = minY = 1e10;
      maxX = maxY = -1e10;
      QRectF rect;
      QPointF point;
      for (auto it=curves.begin(); it!= curves.end(); it++)
      {
          rect = (*it)->boundingRect();
          point = rect.bottomRight();
          if (point.x() < minX) minX = point.x();
          if (point.y() < minY) minY = point.y();
          if (point.x() > maxX) maxX = point.x();
          if (point.y() > maxY) maxY = point.y();
          point = rect.topLeft();
          if (point.x() < minX) minX = point.x();
          if (point.y() < minY) minY = point.y();
          if (point.x() > maxX) maxX = point.x();
          if (point.y() > maxY) maxY = point.y();
      }
      QPointF bottomRight(maxX, minY);
      QPointF topLeft(minX, maxY);
      QPointF center((bottomRight+topLeft)/2.f);
      stack.push(QRectF(topLeft + (topLeft-center)*(.5f),bottomRight + (bottomRight-center)*(.5f)));
      zoom.setZoomStack(stack);

//      delete[] fitsBuffer;
      for (auto it = columnsData.begin(); it != columnsData.end(); it++)
          delete[] it->second;
    window.resize(600, 400);
    window.show();

    a.exec();


    for (auto it = curves.begin(); it != curves.end(); it++)
    {
       (*it)->detach();
        delete *it;
    }
    grid->detach();
    for (auto it = xValues.begin(); it != xValues.end(); it++)
        delete[] *it;

    delete[] yValues;

    delete[] handle;

    return 0;
}
#endif

void SetupConfiguration(Configuration& conf)
{
    po::options_description configs("Fitsdump options");
    configs.add_options()
        ("fitsfile,f",  var<string>()
#if BOOST_VERSION >= 104200
         ->required()
#endif
                                                  , "Name of FITS file")
        ("tablename,t", var<string>("DATA")
#if BOOST_VERSION >= 104200
         ->required()
#endif
                                                  , "Name of input table")
        ("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")
        ("outfile,o",   var<string>("/dev/stdout"), "Name of output file (-:/dev/stdout)")
        ("precision,p", var<int>(20),               "Precision of ofstream")
        ("list,l",      po_switch(),                "List all tables and columns in file")
        ("header,h",    po_switch(),                "Dump header of given table")
        ("stat,s",      po_switch(),                "Perform statistics instead of dump")
        ("minmax,m",    po_switch(),                "Calculates min and max of data")
        ("nozero,z",    po_switch(),                "skip 0 values for stats")
        ("tstart,a",    po_switch(),                "Give the mjdStart from reading the file data")
        ("tstop,b",     po_switch(),                "Give the mjdStop from reading the file data")
#ifdef PLOTTING_PLEASE
        ("graph,g",     po_switch(),                "Plot the columns instead of dumping them")
        ("dots,d",      po_switch(),                "Plot using dots instead of lines")
#endif
        ;

    po::positional_options_description p;
    p.add("fitsfile",   1); // The first positional options
    p.add("col",       -1); // All others

    conf.AddOptions(configs);
    conf.SetArgumentPositions(p);
}

int main(int argc, const char** argv)
{
    Configuration conf(argv[0]);
    conf.SetPrintUsage(PrintUsage);
    SetupConfiguration(conf);

    if (!conf.DoParse(argc, argv, PrintHelp))
        return -1;

    FitsDumper loader;
    return loader.ExecConfig(conf);
}