Index: trunk/FACT++/src/fitsdump.cc =================================================================== --- trunk/FACT++/src/fitsdump.cc (revision 12883) +++ trunk/FACT++/src/fitsdump.cc (revision 12884) @@ -8,81 +8,60 @@ #include "Configuration.h" +#include + #include #include +#include + +#include "Time.h" #include "externals/fits.h" -//#define PLOTTING_PLEASE - -#ifdef PLOTTING_PLEASE -#include -#include -#include - -#include -#include -#include -#include -#include -#include - -#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()); - } +struct MyColumn +{ + string name; + + fits::Table::Column col; + + uint32_t first; + uint32_t last; + + void *ptr; }; -#endif - -class FitsDumper -{ -public: - FitsDumper(); - ~FitsDumper(); - + +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 (valmax) + max = val; + + numValues++; + } +}; + + +class FitsDumper : public fits +{ 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 CalculateOffsets() const; - - template - T PtrToValue(const unsigned char* &ptr) const; -// template -// 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 &, const vector &, unsigned char *, const vector >&) const; - - bool OpenFile(const string &, bool); /// Open a file - bool OpenTable(const string &); /// Open a table /// Lists all columns of an open file @@ -91,21 +70,16 @@ void ListKeywords(ostream &); - bool separateColumnsFromRanges(const vector& list, - vector >& ranges, - vector& listNamesOnly); - /// Perform the dumping, based on the current dump list - bool Dump(const string &, const vector &list, int); + vector InitColumns(const vector& list); + ///Display the selected columns values VS time -#ifdef PLOTTING_PLEASE - int doCurvesDisplay( const vector &list, const string& tableName); -#endif - int doMinMaxPlease(const string& filename, const vector& list, int precision); - int doStatsPlease(const string &filename, const vector& list, int precision); -// void doTBoundary(conf.Get("outfile"), conf.Get("precision"), true); - // bool Plot(const vector &list); + void Dump(ofstream &, const vector &, const string &); + void DumpMinMax(ofstream &, const vector &, bool); + void DumpStats(ofstream &, const vector &); public: + FitsDumper(const string &fname); + ///Configures the fitsLoader from the config file and/or command arguments. - int ExecConfig(Configuration& conf); + int Exec(Configuration& conf); }; @@ -116,18 +90,7 @@ //! the ostream where to redirect the outputs // -FitsDumper::FitsDumper() : fFile(0), fDotsPlease(false), fNoZeroPlease(false) -{ -} - -// -------------------------------------------------------------------------- -// -//! Destructor -// -FitsDumper::~FitsDumper() -{ - if (fFile) - delete fFile; -} - +FitsDumper::FitsDumper(const string &fname) : fits(fname), fFilename(fname) +{ +} string FitsDumper::ValueTypeToStr(char type) const @@ -148,86 +111,8 @@ } -int FitsDumper::ValueTypeToSize(char type) const -{ - switch (type) - { - case 'A': return sizeof(int8_t); - 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 -T FitsDumper::PtrToValue(const unsigned char* &ptr) const -{ - T t; - reverse_copy(ptr, ptr+sizeof(T), reinterpret_cast(&t)); - ptr += sizeof(T); - - return t; -} -// -------------------------------------------------------------------------- -// -//! Loads the fits file based on the current parameters -// -bool FitsDumper::OpenFile(const string &filename, bool force) -{ - if (fFile) - { - fFile->close(); - delete fFile; - } - - try { - fFile = new fits(filename, force); - } - 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; - } + const fits::Table::Keys &fKeyMap = GetKeys(); + const fits::Table::Columns &fColMap = GetColumns(); cout << "\nFile: " << fFilename << "\n"; @@ -245,9 +130,10 @@ cout << endl; - cout << flush; } void FitsDumper::ListKeywords(ostream &out) { + const fits::Table::Keys &fKeyMap = GetKeys(); + for (auto it=fKeyMap.begin(); it != fKeyMap.end(); it++) { out << "## " << setw(8) << it->first << " = " << setw(10); @@ -258,9 +144,4 @@ void FitsDumper::ListHeader(const string& filename) { - if (!fFile) - { - cerr << "No table open." << endl; - return; - } ofstream out(filename=="-"?"/dev/stdout":filename); if (!out) @@ -269,4 +150,6 @@ return; } + + const fits::Table::Keys &fKeyMap = GetKeys(); out << "\nTable: " << fKeyMap.find("EXTNAME")->second.value << " (rows=" << fKeyMap.find("NAXIS2")->second.value << ")\n"; @@ -279,96 +162,88 @@ } -bool FitsDumper::separateColumnsFromRanges(const vector& list, - vector >& ranges, - vector& listNamesOnly) -{ - for (vector::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; -} +vector FitsDumper::InitColumns(const vector &names) +{ + static const boost::regex expr("([[:word:].]+)(\\[([[:digit:]]+)?(:)?([[:digit:]]+)?\\])?"); + + const fits::Table::Columns &fColMap = GetColumns(); + + vector 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(); + } + + 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(); + } + + 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(); + } + + if (last>=col.num) + { + cerr << "ERROR - Last index " << last << " for column " << name << " exceeds number of elements " << col.num << endl; + return vector(); + } + + if (first>last) + { + cerr << "ERROR - Last index " << last << " for column " << name << " exceeds first index " << first << endl; + return vector(); + } + + 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 // -bool FitsDumper::Dump(const string &filename, const vector &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 > ranges; - vector 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); +void FitsDumper::Dump(ofstream &out, const vector &cols, const string &filename) +{ + const fits::Table::Keys &fKeyMap = GetKeys(); out << "## --------------------------------------------------------------------------\n"; + out << "## Fits file:\t" << fFilename << '\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 << "## NumRows: \t" << GetInt("NAXIS2") << '\n'; + out << "## Comment: \t" << ((fKeyMap.find("COMMENT") != fKeyMap.end()) ? fKeyMap.find("COMMENT")->second.value : "") << '\n'; out << "## --------------------------------------------------------------------------\n"; ListKeywords(out); @@ -376,76 +251,58 @@ out << "#\n"; - auto rangesIt = ranges.begin(); - for (vector::const_iterator it=listNamesOnly.begin(); it!=listNamesOnly.end(); it++, rangesIt++) - { - const fits::Table::Column& col = fColMap[*it]; -// const MyColumn *col = static_cast(fTable->column()[*it]); - if (rangesIt->first != 0 || rangesIt->second != (int)(col.num)) - { - out << "#"; - for (int i=rangesIt->first; isecond; i++) - out << " " << *it << "[" << i << "]"; - out << ": " << col.unit; + for (auto it=cols.begin(); it!=cols.end(); it++) + { + out << "# " << it->name; + + if (it->first==it->last) + { + if (it->first!=0) + out << "[" << it->first << "]"; } else - out << "# " << *it << "[" << col.num << "]: " << col.unit; - -// FIXME: retrive the column comment -// if (!col->comment().empty()) -// out << " (" <comment() << ")"; - out << '\n'; + out << "[" << it->first << ":" << it->last << "]"; + + out << ": " << it->col.unit << '\n'; } out << "#" << endl; - - vector ptrs; - for (auto namesIt=listNamesOnly.begin(); namesIt!=listNamesOnly.end(); namesIt++) - { - ptrs.push_back(fFile->SetPtrAddress(*namesIt)); - } - - while (fFile->GetNextRow()) - { - const size_t row = fFile->GetRow(); - if (row==fFile->GetNumRows()) - break; - - rangesIt = ranges.begin(); - auto ptrsIt=ptrs.begin(); - - for (vector::const_iterator in=listNamesOnly.begin(); in!=listNamesOnly.end(); in++, rangesIt++, ptrsIt++) - { - const void *ptr = *ptrsIt; - - const fits::Table::Column &cCol = fColMap.find(*in)->second; - + // ----------------------------------------------------------------- + + while (GetNextRow()) + { + const size_t row = GetRow(); + if (row==GetNumRows()) + break; + + for (auto it=cols.begin(); it!=cols.end(); it++) + { string msg; - for (int i=rangesIt->first; isecond; i++) + for (uint32_t i=it->first; i<=it->last; i++) { - switch (cCol.type) + switch (it->col.type) { case 'A': - msg += reinterpret_cast(ptr)[i]; + msg += reinterpret_cast(it->ptr)[i]; break; case 'B': - out << (unsigned int)reinterpret_cast(ptr)[i] << " "; + out << (unsigned int)reinterpret_cast(it->ptr)[i] << " "; break; case 'L': - out << reinterpret_cast(ptr)[i] << " "; + out << reinterpret_cast(it->ptr)[i] << " "; break; case 'I': - out << reinterpret_cast(ptr)[i] << " "; + out << reinterpret_cast(it->ptr)[i] << " "; break; case 'J': - out << reinterpret_cast(ptr)[i] << " "; + out << reinterpret_cast(it->ptr)[i] << " "; break; case 'K': - out << reinterpret_cast(ptr)[i] << " "; + out << reinterpret_cast(it->ptr)[i] << " "; break; case 'E': - out << reinterpret_cast(ptr)[i] << " "; + out << reinterpret_cast(it->ptr)[i] << " "; break; case 'D': - out << reinterpret_cast(ptr)[i] << " "; + out << reinterpret_cast(it->ptr)[i] << " "; break; default: @@ -454,10 +311,205 @@ } - if (cCol.type=='A') + if (it->col.type=='A') out << "'" << msg << "' "; } out << endl; } - return true; +} + +void FitsDumper::DumpMinMax(ofstream &out, const vector &cols, bool fNoZeroPlease) +{ + vector 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 in=cols.begin(); in!=cols.end(); in++, statsIt++) + { + if ((in->name=="UnixTimeUTC" || in->name=="PCTime") && in->first==0 && in->last==1) + { + const uint32_t *val = reinterpret_cast(in->ptr); + if (fNoZeroPlease && val[0]==0 && val[1]==0) + continue; + + statsIt->add(Time(val[0], val[1]).Mjd()); + continue; + } + + for (uint32_t i=in->first; i<=in->last; i++) + { + double cValue = 0; + switch (in->col.type) + { + case 'L': + cValue = reinterpret_cast(in->ptr)[i]; + break; + case 'B': + cValue = reinterpret_cast(in->ptr)[i]; + break; + case 'I': + cValue = reinterpret_cast(in->ptr)[i]; + break; + case 'J': + cValue = reinterpret_cast(in->ptr)[i]; + break; + case 'K': + cValue = reinterpret_cast(in->ptr)[i]; + break; + case 'E': + cValue = reinterpret_cast(in->ptr)[i]; + break; + case 'D': + cValue = reinterpret_cast(in->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++) + { + cout << "\n[" << it->name << ':' << it->first; + if (it->first!=it->last) + cout << ':' << it->last; + cout << "]\n"; + + if (statsIt->numValues == 0) + { + out << "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); + + out << "Min: " << statsIt->min << '\n'; + out << "Max: " << statsIt->max << '\n'; + out << "Avg: " << avg << '\n'; + out << "Rms: " << rms << endl; + } +} + +template +void displayStats(vector &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(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 &cols) +{ + // Loop over all columns in our list of requested columns + vector> statData; + + for (auto in=cols.begin(); in!=cols.end(); in++) + statData.push_back(vector(in->col.size*GetNumRows()*(in->last-in->first+1))); + + while (GetNextRow()) + { + const size_t row = GetRow(); + if (row==GetNumRows()) + break; + + auto statsIt = statData.begin(); + for (auto in=cols.begin(); in!=cols.end(); in++, statsIt++) + { + const char *src = reinterpret_cast(in->ptr); + const size_t sz = (in->last-in->first+1)*in->col.size; + memcpy(statsIt->data()+row*sz, src+in->first*in->col.size, sz); + } + } + + auto statsIt = statData.begin(); + for (auto in=cols.begin(); in!=cols.end(); in++, statsIt++) + { + out << "\n[" << in->name << ':' << in->first; + if (in->last!=in->first) + out << ':' << in->last; + out << "]\n"; + + switch (in->col.type) + { + case 'L': + displayStats(*statsIt, out); + break; + case 'B': + displayStats(*statsIt, out); + break; + case 'I': + displayStats(*statsIt, out); + break; + case 'J': + displayStats(*statsIt, out); + break; + case 'K': + displayStats(*statsIt, out); + break; + case 'E': + displayStats(*statsIt, out); + break; + case 'D': + displayStats(*statsIt, out); + break; + default: + ; + } + } } @@ -468,116 +520,64 @@ //! the configuration object // -int FitsDumper::ExecConfig(Configuration& conf) -{ - if (conf.Has("fitsfile")) - { - if (!OpenFile(conf.Get("fitsfile"), conf.Get("force"))) - return -1; - } -#ifdef PLOTTING_PLEASE - if (conf.Get("stat") && conf.Get("graph")) - { - cout << "Invalid conbination of options: cannot graph stats. Aborting" << endl; - return -1; - } -#endif - if (conf.Get("minmax") && conf.Get("stat")) - { - cout << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl; - return -1; - } -#ifdef PLOTTING_PLEASE - if (conf.Get("minmax") && conf.Get("graph")) - { - cout << "Invalid combination of options: cannot graph minmax. Aborting" << endl; - return -1; - } -#endif - if (conf.Get("stat") && conf.Get("nozero")) - { - cout << "Invalid combination of options: nozero only works with minmax. Aborting" << endl; - return -1; - } - - if (conf.Get("nozero")) - { - fNoZeroPlease = true; - } - - if (conf.Has("tablename")) - { - if (!OpenTable(conf.Get("tablename"))) - return -1; - } - +int FitsDumper::Exec(Configuration& conf) +{ if (conf.Get("list")) List(); - /* - if (conf.Get("tstart")) - { - doTBoundary(conf.Get("outfile"), conf.Get("precision"), true); - } - if (conf.Get("tstop")) - { - doTBoundary(conf.Get("outfile"), conf.Get("precision"), false); - } - */ - if (conf.Get("minmax")) - { - if (!conf.Has("col")) - { - cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl; - return 0; - } - doMinMaxPlease(conf.Get("outfile"), conf.Get>("col"), conf.Get("precision")); - return 0; - } - - if (conf.Get("stat")) - { - if (!conf.Has("col")) - { - cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl; - return 0; - } - doStatsPlease(conf.Get("outfile"), conf.Get>("col"), conf.Get("precision")); - return 0; - } - -#ifdef PLOTTING_PLEASE - if (conf.Get("graph")) - { - if (!conf.Has("col")) - { - cout << "Please specify the columns that should be dumped as arguments. Aborting" << endl; - return 0; - } - if (conf.Get("dots")) - fDotsPlease = true; - doCurvesDisplay(conf.Get>("col"), - conf.Get("tablename")); - return 1; - } -#endif if (conf.Get("header")) ListHeader(conf.Get("outfile")); + if (conf.Get("header") || conf.Get("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("outfile"), - conf.Get>("col"), - conf.Get("precision"))) - return -1; - } - + // ------------------------------------------------------------ + + if (conf.Get("minmax") && conf.Get("stat")) + { + cerr << "Invalid combination of options: cannot do stats and minmax. Aborting" << endl; + return -1; + } + if (conf.Get("stat") && conf.Get("nozero")) + { + cerr << "Invalid combination of options: nozero only works with minmax. Aborting" << endl; + return -1; + } + + // ------------------------------------------------------------ + + if (conf.Vec("col").size()==0) + { + cerr << "No columns specifiec." << endl; + return 0; + } + + const string filename = conf.Get("outfile"); + + ofstream out(filename=="-"?"/dev/stdout":filename); + if (!out) + { + cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl; + return false; + } + out.precision(conf.Get("precision")); + + const vector cols = InitColumns(conf.Vec("col")); + if (cols.size()==0) + return false; + + if (conf.Get("minmax")) + { + DumpMinMax(out, cols, conf.Get("nozero")); + return 0; + } + + if (conf.Get("stat")) + { + DumpStats(out, cols); + return 0; + } + + Dump(out, cols, filename); return 0; @@ -598,802 +598,4 @@ // } - -struct minMaxStruct { - double min; - double max; - double average; - double squared; - long numValues; - minMaxStruct() : min(1e10), max(-1e10), average(0), numValues(0) {} -}; - -int FitsDumper::doMinMaxPlease(const string& filename, const vector& 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 > ranges; - vector names; - - if (!separateColumnsFromRanges(list, ranges, names)) - { - 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); - - vector statData(names.size()); - - vector ptrs; - for (auto namesIt=names.begin(); namesIt!=names.end(); namesIt++) - { - ptrs.push_back(fFile->SetPtrAddress(*namesIt)); - } - - // Loop over all columns in our list of requested columns - while (fFile->GetNextRow()) - { - const size_t row = fFile->GetRow(); - if (row==fFile->GetNumRows()) - break; - - auto rangesIt = ranges.begin(); - auto statsIt = statData.begin(); - auto ptrsIt = ptrs.begin(); - - for (auto in=names.begin(); in!=names.end(); in++, rangesIt++, statsIt++, ptrsIt++) - { - const void *ptr = *ptrsIt; - - const fits::Table::Column &cCol = fColMap.find(*in)->second; - - if (*in=="UnixTimeUTC" || *in=="PCTime") - { - const uint32_t *val = reinterpret_cast(ptr); - if (fNoZeroPlease && val[0]==0 && val[1]==0) - continue; - - const double cValue = Time(val[0], val[1]).Mjd(); - - statsIt->average += cValue; - statsIt->squared += cValue*cValue; - - if (cValue < statsIt->min) - statsIt->min = cValue; - - if (cValue > statsIt->max) - statsIt->max = cValue; - - statsIt->numValues++; - - continue; - } - - for (int i=rangesIt->first; isecond; i++) - { - double cValue = 0; - switch (cCol.type) - { - case 'L': - cValue = reinterpret_cast(ptr)[i]; - break; - case 'B': - cValue = reinterpret_cast(ptr)[i]; - break; - case 'I': - cValue = reinterpret_cast(ptr)[i]; - break; - case 'J': - cValue = reinterpret_cast(ptr)[i]; - break; - case 'K': - cValue = reinterpret_cast(ptr)[i]; - break; - case 'E': - cValue = reinterpret_cast(ptr)[i]; - break; - case 'D': - cValue = reinterpret_cast(ptr)[i]; - break; - default: - ; - } - - if (fNoZeroPlease && cValue == 0) - continue; - - statsIt->average += cValue; - statsIt->squared += cValue*cValue; - - if (cValue < statsIt->min) - statsIt->min = cValue; - - if (cValue > statsIt->max) - statsIt->max = cValue; - - statsIt->numValues++; - } - } - } - - // 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=names.begin(); it!=names.end(); it++, statsIt++) - { - cout << "[" << *it << "]" << endl; - if (statsIt->numValues == 0) - { - out << "Min: -\nMax: -\nAvg: -\nRms: -" << endl; - continue; - } - - double &avg = statsIt->average; - double &rms = statsIt->squared; - long &num = statsIt->numValues; - - avg /= num; - rms = sqrt(rms/num - avg*avg); - - out << "Min: " << statsIt->min << '\n'; - out << "Max: " << statsIt->max << '\n'; - out << "Avg: " << avg << '\n'; - out << "Rms: " << rms << endl; - } - - /* - vector > columnsData; - vector statData; - int numRows = fFile->GetInt("NAXIS2"); - auto rangesIt = ranges.begin(); - for (vector::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; isecond; i++) - { - switch (it->first) { - case 'L': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'B': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'I': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'J': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'K': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'E': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'D': - cValue = reinterpret_cast(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; -*/ - 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 > ranges; - vector 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 > columnsData; - vector statData; - int numRows = fFile->GetInt("NAXIS2"); - auto rangesIt = ranges.begin(); - for (vector::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; isecond; i++) - { - switch (it->first) { - case 'L': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'B': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'I': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'J': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'K': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'E': - cValue = reinterpret_cast(it->second)[i]; - break; - case 'D': - cValue = reinterpret_cast(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 -void displayStats(char* array, int numElems, ofstream& out) -{ - if (numElems == 0) - { - out << "Min: -\nMax: -\nMed: -\nAvg: -\nRms: -" << endl; - return; - } - - T *val = reinterpret_cast(array); - - 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-1)/2] + val[(numElems-1)/2+1])/2 << '\n'; - else - out << "Med: " << double(val[numElems/2+1]) << '\n'; - } - - double avg = 0; - double rms = 0; - for (int i=0;i& 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 > ranges; - vector names; - - if (!separateColumnsFromRanges(list, ranges, names)) - { - 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 statData; - - auto rangesIt = ranges.begin(); - for (auto in=names.begin(); in!=names.end(); in++, rangesIt++) - { - cout << *in << endl; - const fits::Table::Column &cCol = fColMap.find(*in)->second; - statData.push_back(new char[(rangesIt->second - rangesIt->first)*cCol.size*fFile->GetNumRows()]); - } - vector ptrs; - for (auto namesIt=names.begin(); namesIt!=names.end(); namesIt++) - { - ptrs.push_back(fFile->SetPtrAddress(*namesIt)); - } - - while (fFile->GetNextRow()) - { - const size_t row = fFile->GetRow(); - if (row==fFile->GetNumRows()) - break; - - rangesIt = ranges.begin(); - auto statsIt = statData.begin(); - auto ptrsIt = ptrs.begin(); - - for (auto in=names.begin(); in!=names.end(); in++, rangesIt++, statsIt++, ptrsIt++) - { - const void *ptr = *ptrsIt; - - const fits::Table::Column &cCol = fColMap.find(*in)->second; - - const int span = rangesIt->second - rangesIt->first; - - for (int i=rangesIt->first; isecond; i++) - { - switch (cCol.type) - { - case 'L': - reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(ptr)[i]; - break; - case 'B': - reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(ptr)[i]; - break; - case 'I': - reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(ptr)[i]; - break; - case 'J': - reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(ptr)[i]; - break; - case 'K': - reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(ptr)[i]; - break; - case 'E': - reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(ptr)[i]; - break; - case 'D': - reinterpret_cast(*statsIt)[i - rangesIt->first + row*span] = reinterpret_cast(ptr)[i]; - break; - default: - ; - } - } - } - } - - //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(); - - for (auto it=list.begin(); it!=list.end(); it++, rangesIt++, statsIt++) - { - const fits::Table::Column &cCol = fColMap.find(*it)->second; - - const int span = rangesIt->second - rangesIt->first; - - out << "[" << *it << "]" << endl; - switch (cCol.type) - { - case 'L': - displayStats(*statsIt, fFile->GetNumRows()*span, out); - break; - case 'B': - displayStats(*statsIt, fFile->GetNumRows()*span, out); - break; - case 'I': - displayStats(*statsIt, fFile->GetNumRows()*span, out); - break; - case 'J': - displayStats(*statsIt, fFile->GetNumRows()*span, out); - break; - case 'K': - displayStats(*statsIt, fFile->GetNumRows()*span, out); - break; - case 'E': - displayStats(*statsIt, fFile->GetNumRows()*span, out); - break; - case 'D': - displayStats(*statsIt, fFile->GetNumRows()*span, out); - break; - default: - ; - } - } - - return true; -} -#ifdef PLOTTING_PLEASE -int FitsDumper::doCurvesDisplay( const vector &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 > ranges; - vector listNamesOnly; - if (!separateColumnsFromRanges(list, ranges, listNamesOnly)) - { - cerr << "Something went wrong while extracting the columns names from parameters. Aborting" << endl; - return false; - } - vector curvesNames; - stringstream str; - for (auto it=ranges.begin(), jt=listNamesOnly.begin(); it != ranges.end(); it++, jt++) - { - for (int i=it->first; isecond;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 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 > columnsData; - - for (vector::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 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;iGetNextRow(); - 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; jsecond; j++) - { - switch (it->first) { - case 'L': - str << reinterpret_cast(it->second)[j] << " "; - break; - case 'B': - str << reinterpret_cast(it->second)[j] << " "; - break; - case 'I': - str << reinterpret_cast(it->second)[j] << " "; - break; - case 'J': - str << reinterpret_cast(it->second)[j] << " "; - break; - case 'K': - str << reinterpret_cast(it->second)[j] << " "; - break; - case 'E': - str << reinterpret_cast(it->second)[j] << " "; - break; - case 'D': - str << reinterpret_cast(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 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) @@ -1406,9 +608,4 @@ #endif , "Name of FITS file") - ("tablename,t", var("DATA") -#if BOOST_VERSION >= 104200 - ->required() -#endif - , "Name of input table") ("col,c", vars(), "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("/dev/stdout"), "Name of output file (-:/dev/stdout)") @@ -1419,16 +616,10 @@ ("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") ("force", po_switch(), "Force reading the fits file even if END key is missing") -#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 + p.add("fitsfile", 1); // The first positional options + p.add("col", -1); // All others conf.AddOptions(configs); @@ -1445,5 +636,18 @@ return -1; - FitsDumper loader; - return loader.ExecConfig(conf); -} + if (!conf.Has("fitsfile")) + { + cerr << "Filename required." << endl; + return -1; + } + + FitsDumper loader(conf.Get("fitsfile")); + if (!loader) + { + cerr << "ERROR - Opening " << conf.Get("fitsfile"); + cerr << " failed: " << strerror(errno) << endl; + return -1; + } + + return loader.Exec(conf); +}