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

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

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

#include <float.h>

#include <map>
#include <fstream>

#include <boost/regex.hpp>

#include "Time.h"
#include "externals/fits.h"

using namespace std;

struct MyColumn
{
    string name;

    fits::Table::Column col;

    uint32_t first;
    uint32_t last;

    void *ptr;
};

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

    void add(double val)
    {
        average += val;
        squared += val*val;

        if (val<min)
            min = val;

        if (val>max)
            max = val;

        numValues++;
    }
};


class FitsDumper : public fits
{
private:
    string fFilename;

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

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

    vector<MyColumn> InitColumns(vector<string> list);

    ///Display the selected columns values VS time
    void Dump(ofstream &, const vector<MyColumn> &, const string &);
    void DumpMinMax(ofstream &, const vector<MyColumn> &, bool);
    void DumpStats(ofstream &, const vector<MyColumn> &);

public:
    FitsDumper(const string &fname);

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

// --------------------------------------------------------------------------
//
//! Constructor
//! @param out
//!        the ostream where to redirect the outputs
//
FitsDumper::FitsDumper(const string &fname) : fits(fname), fFilename(fname)
{
}

string FitsDumper::ValueTypeToStr(char type) const
{
    switch (type)
    {
        case 'L': return "bool(8)";
        case 'A': return "char(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 "unknown";
    }
}

void FitsDumper::List()
{
    const fits::Table::Keys    &fKeyMap = GetKeys();
    const fits::Table::Columns &fColMap = GetColumns();

    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;
}

void FitsDumper::ListKeywords(ostream &fout)
{
    const fits::Table::Keys &fKeyMap = GetKeys();

    for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++)
    {
        fout << "## " << ::left << setw(8) << it->first << "= ";

        if (it->second.type=='T')
            fout << ::left  << setw(20) << ("'"+it->second.value+"'");
        else
            fout << ::right << setw(20) << it->second.value;

        if (it->second.comment.size()>0)
            fout << " / " << it->second.comment;
        fout << '\n';
    }

    fout << flush;
}

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

    const fits::Table::Keys &fKeyMap = GetKeys();

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

    ListKeywords(fout);
    fout << endl;
}

vector<MyColumn> FitsDumper::InitColumns(vector<string> names)
{
    static const boost::regex expr("([[:word:].]+)(\\[([[:digit:]]+)?(:)?([[:digit:]]+)?\\])?");

    const fits::Table::Columns &fColMap = GetColumns();

    if (names.size()==0)
        for (auto it=fColMap.begin(); it!=fColMap.end(); it++)
            if (it->second.num>0)
                names.push_back(it->first);

    vector<MyColumn> vec;

    for (auto it=names.begin(); it!=names.end(); it++)
    {
        boost::smatch what;
        if (!boost::regex_match(*it, what, expr, boost::match_extra))
        {
            cerr << "Couldn't parse expression '" << *it << "' " << endl;
            return vector<MyColumn>();
        }

        const string name = what[1];

        const auto iter = fColMap.find(name);
        if (iter==fColMap.end())
        {
            cerr << "ERROR - Column '" << name << "' not found in table." << endl;
            return vector<MyColumn>();
        }

        const fits::Table::Column &col = iter->second;

        const string val0  = what[3];
        const string delim = what[4];
        const string val1  = what[5];

        const uint32_t first = val0.empty() ? 0 : atoi(val0.c_str());
        const uint32_t last  = (val0.empty() && delim.empty()) ? col.num-1 : (val1.empty() ? first : atoi(val1.c_str()));

        if (first>=col.num)
        {
            cerr << "ERROR - First index " << first << " for column " << name << " exceeds number of elements " << col.num << endl;
            return vector<MyColumn>();
        }

        if (last>=col.num)
        {
            cerr << "ERROR - Last index " << last << " for column " << name << " exceeds number of elements " << col.num << endl;
            return vector<MyColumn>();
        }

        if (first>last)
        {
            cerr << "ERROR - Last index " << last << " for column " << name << " exceeds first index " << first << endl;
            return vector<MyColumn>();
        }

        MyColumn mycol;

        mycol.name  = name;
        mycol.col   = col;
        mycol.first = first;
        mycol.last  = last;

        vec.push_back(mycol);
    }

    for (auto it=vec.begin(); it!=vec.end(); it++)
        it->ptr = SetPtrAddress(it->name);

    return vec;
}

// --------------------------------------------------------------------------
//
//! Perform the actual dump, based on the current parameters
//
void FitsDumper::Dump(ofstream &fout, const vector<MyColumn> &cols, const string &filename)
{
    const fits::Table::Keys &fKeyMap = GetKeys();

    fout << "## --------------------------------------------------------------------------\n";
    fout << "## Fits file:\t" << fFilename << '\n';
    if (filename!="-")
        fout << "## File:    \t" << filename << '\n';
    fout << "## Table:   \t" << fKeyMap.find("EXTNAME")->second.value << '\n';
    fout << "## NumRows: \t" << GetInt("NAXIS2") << '\n';
    fout << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n';
    fout << "## --------------------------------------------------------------------------\n";
    ListKeywords(fout);
    fout << "## --------------------------------------------------------------------------\n";
    fout << "#\n";

    for (auto it=cols.begin(); it!=cols.end(); it++)
    {
        fout << "# " << it->name;

        if (it->first==it->last)
        {
            if (it->first!=0)
                fout << "[" << it->first << "]";
        }
        else
            fout << "[" << it->first << ":" << it->last << "]";

        fout << ": " << it->col.unit << '\n';
    }
    fout << "#" << endl;

    // -----------------------------------------------------------------

    while (GetNextRow())
    {
        const size_t row = GetRow();
        if (row==GetNumRows())
            break;

        for (auto it=cols.begin(); it!=cols.end(); it++)
        {
            string msg;
            for (uint32_t i=it->first; i<=it->last; i++)
            {
                switch (it->col.type)
                {
                case 'A':
                    msg += reinterpret_cast<const char*>(it->ptr)[i];
                    break;
                case 'B':
                    fout << (unsigned int)reinterpret_cast<const unsigned char*>(it->ptr)[i] << " ";
                    break;
                case 'L':
                    fout << reinterpret_cast<const bool*>(it->ptr)[i] << " ";
                    break;
                case 'I':
                    fout << reinterpret_cast<const int16_t*>(it->ptr)[i] << " ";
                    break;
                case 'J':
                    fout << reinterpret_cast<const int32_t*>(it->ptr)[i] << " ";
                    break;
                case 'K':
                    fout << reinterpret_cast<const int64_t*>(it->ptr)[i] << " ";
                    break;
                case 'E':
                    fout << reinterpret_cast<const float*>(it->ptr)[i] << " ";
                    break;
                case 'D':
                    fout << reinterpret_cast<const double*>(it->ptr)[i] << " ";
                    break;
                default:
                    ;
                }
            }

            if (it->col.type=='A')
                fout << "'" << msg << "' ";
        }
        fout << endl;
    }
}

void FitsDumper::DumpMinMax(ofstream &fout, const vector<MyColumn> &cols, bool fNoZeroPlease)
{
    vector<minMaxStruct> statData(cols.size());

    // Loop over all columns in our list of requested columns
    while (GetNextRow())
    {
        const size_t row = GetRow();
        if (row==GetNumRows())
            break;

        auto statsIt = statData.begin();

        for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
        {
            if ((it->name=="UnixTimeUTC" || it->name=="PCTime") && it->first==0 && it->last==1)
            {
                const uint32_t *val = reinterpret_cast<const uint32_t*>(it->ptr);
                if (fNoZeroPlease && val[0]==0 && val[1]==0)
                    continue;

                statsIt->add(Time(val[0], val[1]).Mjd());
                continue;
            }

            for (uint32_t i=it->first; i<=it->last; i++)
            {
                double cValue = 0;
                switch (it->col.type)
                {
                case 'L':
                        cValue = reinterpret_cast<const bool*>(it->ptr)[i];
                        break;
                case 'B':
                        cValue = reinterpret_cast<const int8_t*>(it->ptr)[i];
                        break;
                case 'I':
                        cValue = reinterpret_cast<const int16_t*>(it->ptr)[i];
                        break;
                case 'J':
                        cValue = reinterpret_cast<const int32_t*>(it->ptr)[i];
                        break;
                case 'K':
                        cValue = reinterpret_cast<const int64_t*>(it->ptr)[i];
                        break;
                case 'E':
                        cValue = reinterpret_cast<const float*>(it->ptr)[i];
                        break;
                case 'D':
                        cValue = reinterpret_cast<const double*>(it->ptr)[i];
                        break;
                default:
                    ;
                }

                if (fNoZeroPlease && cValue == 0)
                    continue;

                statsIt->add(cValue);
            }
        }
    }

    // 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)
    auto statsIt = statData.begin();
    for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
    {
        fout << "\n[" << it->name << ':' << it->first;
        if (it->first!=it->last)
            fout << ':' << it->last;
        fout << "]\n";

        if (statsIt->numValues == 0)
        {
            fout << "Min: -\nMax: -\nAvg: -\nRms: -" << endl;
            continue;
        }

        const long &num = statsIt->numValues;

        long double &avg = statsIt->average;
        long double &rms = statsIt->squared;

        avg /= num;
        rms  = sqrt(rms/num - avg*avg);

        fout << "Min: " << statsIt->min << '\n';
        fout << "Max: " << statsIt->max << '\n';
        fout << "Avg: " << avg << '\n';
        fout << "Rms: " << rms << endl;
    }
}

template<typename T>
void displayStats(vector<char> &array, ofstream& out)
{
    const size_t numElems = array.size()/sizeof(T);
    if (numElems == 0)
    {
        out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl;
        return;
    }

    T *val = reinterpret_cast<T*>(array.data());

    sort(val, val+numElems);

    out << "Min: " << double(val[0]) << '\n';
    out << "Max: " << double(val[numElems-1]) << '\n';

    if (numElems>2)
    {
        if (numElems%2 == 0)
            out << "Med: " << (double(val[numElems/2]) + double(val[numElems/2+1]))/2 << '\n';
        else
            out << "Med: " << double(val[numElems/2+1]) << '\n';
    }

    long double avg = 0;
    long double rms = 0;
    for (uint32_t i=0;i<numElems;i++)
    {
        avg += double(val[i]);
        rms += double(val[i])*double(val[i]);
    }

    avg /= numElems;
    rms  = sqrt(rms/numElems - avg*avg);

    out << "Avg: " << avg << '\n';
    out << "Rms: " << rms << endl;

}

void FitsDumper::DumpStats(ofstream &fout, const vector<MyColumn> &cols)
{
    // Loop over all columns in our list of requested columns
    vector<vector<char>> statData;

    for (auto it=cols.begin(); it!=cols.end(); it++)
        statData.push_back(vector<char>(it->col.size*GetNumRows()*(it->last-it->first+1)));

    while (GetNextRow())
    {
        const size_t row = GetRow();
        if (row==GetNumRows())
            break;

        auto statsIt = statData.begin();
        for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
        {
            const char *src = reinterpret_cast<const char*>(it->ptr);
            const size_t sz = (it->last-it->first+1)*it->col.size;
            memcpy(statsIt->data()+row*sz, src+it->first*it->col.size, sz);
        }
    }

    auto statsIt = statData.begin();
    for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++)
    {
        fout << "\n[" << it->name << ':' << it->first;
        if (it->last!=it->first)
            fout << ':' << it->last;
        fout << "]\n";

        switch (it->col.type)
        {
        case 'L':
            displayStats<bool>(*statsIt, fout);
            break;
        case 'B':
            displayStats<char>(*statsIt, fout);
            break;
        case 'I':
            displayStats<int16_t>(*statsIt, fout);
            break;
        case 'J':
            displayStats<int32_t>(*statsIt, fout);
            break;
        case 'K':
            displayStats<int64_t>(*statsIt, fout);
            break;
        case 'E':
            displayStats<float>(*statsIt, fout);
            break;
        case 'D':
            displayStats<double>(*statsIt, fout);
            break;
        default:
            ;
        }
    }
}

// --------------------------------------------------------------------------
//
//! Retrieves the configuration parameters
//! @param conf
//!             the configuration object
//
int FitsDumper::Exec(Configuration& conf)
{
    if (conf.Get<bool>("list"))
        List();

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


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

    // ------------------------------------------------------------

    if (conf.Get<bool>("minmax") && conf.Get<bool>("stat"))
    {
        cerr << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl;
        return -1;
    }
    if (conf.Get<bool>("stat") && conf.Get<bool>("nozero"))
    {
        cerr << "Invalid combination of options: nozero only works with minmax. Aborting" << endl;
        return -1;
    }

    // ------------------------------------------------------------

    const string filename = conf.Get<string>("outfile");

    ofstream fout(filename=="-"?"/dev/stdout":filename);
    if (!fout)
    {
        cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl;
        return false;
    }
    fout.precision(conf.Get<int>("precision"));

    const vector<MyColumn> cols = InitColumns(conf.Vec<string>("col"));
    if (cols.size()==0)
        return false;

    if (conf.Get<bool>("minmax"))
    {
        DumpMinMax(fout, cols, conf.Get<bool>("nozero"));
        return 0;
    }

    if (conf.Get<bool>("stat"))
    {
        DumpStats(fout, cols);
        return 0;
    }

    Dump(fout, cols, filename);

    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()
{
    // 
}

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")
        ("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")
        ("force",       po_switch(),                "Force reading the fits file even if END key is missing")
        ;

    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;

    if (!conf.Has("fitsfile"))
    {
        cerr << "Filename required." << endl;
        return -1;
    }

    FitsDumper loader(conf.Get<string>("fitsfile"));
    if (!loader)
    {
        cerr << "ERROR - Opening " << conf.Get<string>("fitsfile");
        cerr << " failed: " << strerror(errno) << endl;
        return -1;
    }

    return loader.Exec(conf);
}