//**************************************************************** /** @class FitsDumper @brief Dumps contents of fits tables to stdout or a file */ //**************************************************************** #include "Configuration.h" #include #include #include #include #include "Time.h" #include "externals/fits.h" #ifdef HAVE_ROOT #include "TFormula.h" #endif 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 (valmax) 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 InitColumns(vector list); vector InitColumnsRoot(vector &list); double GetDouble(const MyColumn &, size_t) const; ///Display the selected columns values VS time void Dump(ofstream &, const vector &, const string &, size_t, size_t, const string &); void DumpRoot(ofstream &, const vector &, const string &, size_t, size_t, const string &); void DumpMinMax(ofstream &, const vector &, size_t, size_t, bool); void DumpStats(ofstream &, const vector &, size_t, size_t); 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 FitsDumper::InitColumns(vector 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 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 = atol(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; mycol.ptr = SetPtrAddress(name); vec.push_back(mycol); } return vec; } double FitsDumper::GetDouble(const MyColumn &it, size_t i) const { switch (it.col.type) { case 'A': return reinterpret_cast(it.ptr)[i]; case 'L': return reinterpret_cast(it.ptr)[i]; case 'B': return (unsigned int)reinterpret_cast(it.ptr)[i]; case 'I': return reinterpret_cast(it.ptr)[i]; case 'J': return reinterpret_cast(it.ptr)[i]; case 'K': return reinterpret_cast(it.ptr)[i]; case 'E': return reinterpret_cast(it.ptr)[i]; case 'D': return reinterpret_cast(it.ptr)[i]; } return 0; } // -------------------------------------------------------------------------- // //! Perform the actual dump, based on the current parameters // void FitsDumper::Dump(ofstream &fout, const vector &cols, const string &filter, size_t first, size_t limit, const string &filename) { const fits::Table::Keys &fKeyMap = GetKeys(); #ifdef HAVE_ROOT TFormula select; if (!filter.empty() && select.Compile(filter.c_str())) throw runtime_error("Syntax Error: TFormula::Compile failed for '"+filter+"'"); #endif 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'; #ifdef HAVE_ROOT if (!filter.empty()) fout << "## Selection: \t" << select.GetExpFormula() << '\n'; #endif fout << "## --------------------------------------------------------------------------\n"; ListKeywords(fout); fout << "## --------------------------------------------------------------------------\n"; fout << "#\n"; size_t num = 0; 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 << "]"; if (!it->col.unit.empty()) fout << ": " << it->col.unit; fout << '\n'; num += it->last-it->first+1; } fout << "#" << endl; // ----------------------------------------------------------------- #ifdef HAVE_ROOT vector data(num); #endif const size_t last = limit ? first + limit : size_t(-1); while (GetRow(first++)) { const size_t row = GetRow(); if (row==GetNumRows() || row==last) break; size_t p = 0; #ifdef HAVE_ROOT data[p++] = first-1; #endif ostringstream sout; sout.precision(fout.precision()); for (auto it=cols.begin(); it!=cols.end(); it++) { string msg; for (uint32_t i=it->first; i<=it->last; i++, p++) { switch (it->col.type) { case 'A': msg += reinterpret_cast(it->ptr)[i]; break; case 'B': sout << (unsigned int)reinterpret_cast(it->ptr)[i] << " "; break; case 'L': sout << reinterpret_cast(it->ptr)[i] << " "; break; case 'I': sout << reinterpret_cast(it->ptr)[i] << " "; break; case 'J': sout << reinterpret_cast(it->ptr)[i] << " "; break; case 'K': sout << reinterpret_cast(it->ptr)[i] << " "; break; case 'E': sout << reinterpret_cast(it->ptr)[i] << " "; break; case 'D': sout << reinterpret_cast(it->ptr)[i] << " "; break; default: ; } #ifdef HAVE_ROOT if (!filter.empty()) data[p] = GetDouble(*it, i); #endif } if (it->col.type=='A') sout << "'" << msg << "' "; } #ifdef HAVE_ROOT if (!filter.empty() && select.EvalPar(0, data.data())<0.5) continue; #endif fout << sout.str() << endl; } } vector FitsDumper::InitColumnsRoot(vector &names) { static const boost::regex expr("[^\\[]([[:word:].]+)(\\[([[:digit:]]+)\\])?"); const fits::Table::Columns &cols = GetColumns(); vector vec; for (auto it=names.begin(); it!=names.end(); it++) { if (it->empty()) continue; *it = ' '+*it; string::const_iterator ibeg = it->begin(); string::const_iterator iend = it->end(); boost::smatch what; while (boost::regex_search(ibeg, iend, what, expr, boost::match_extra)) { const string all = what[0]; const string name = what[1]; const size_t idx = atol(string(what[3]).c_str()); // Check if found colum is valid const auto ic = cols.find(name); if (ic==cols.end()) { ibeg++; //cout << "Column '" << name << "' does not exist." << endl; //return vector(); continue; } if (idx>=ic->second.num) { cout << "Column '" << name << "' has no index " << idx << "." << endl; return vector(); } // find index if column already exists size_t p = 0; for (; preplace(ibeg-it->begin()+what.position(1), what.length()-1, id.str()); ibeg = what[0].first+3; iend = it->end(); if (psecond; mycol.first = idx; mycol.last = idx; mycol.ptr = SetPtrAddress(name); vec.push_back(mycol); } } ostringstream id; id << '[' << vec.size() << ']'; for (auto it=names.begin(); it!=names.end(); it++) { while (1) { auto p = it->find_first_of('#'); if (p==string::npos) break; it->replace(p, 1, id.str()); } } //cout << endl; //for (size_t i=0; isecond.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"; if (!filter.empty()) fout << "## Selection: " << form[0].GetExpFormula() << "\n##\n"; size_t num = 0; for (auto it=vec.begin(); it!=vec.end(); it++, num++) { fout << "## [" << num << "] = " << it->name; if (it->first==it->last) { if (it->first!=0) fout << "[" << it->first << "]"; } else fout << "[" << it->first << ":" << it->last << "]"; if (!it->col.unit.empty()) fout << ": " << it->col.unit; fout << '\n'; } fout << "##\n"; fout << "## --------------------------------------------------------------------------\n"; fout << "#\n"; fout << "# "; for (auto it=form.begin()+1; it!=form.end(); it++) fout << " \"" << it->GetExpFormula() << "\""; fout << "\n#" << endl; // ----------------------------------------------------------------- vector data(vec.size()+1); const size_t last = limit ? first + limit : size_t(-1); while (GetRow(first++)) { const size_t row = GetRow(); if (row==GetNumRows() || row==last) break; size_t p = 0; for (auto it=vec.begin(); it!=vec.end(); it++, p++) data[p] = GetDouble(*it, it->first); data[p] = first; if (!filter.empty() && form[0].EvalPar(0, data.data())<0.5) continue; for (auto iform=form.begin()+1; iform!=form.end(); iform++) fout << iform->EvalPar(0, data.data()) << " "; fout << endl; } #endif } void FitsDumper::DumpMinMax(ofstream &fout, const vector &cols, size_t first, size_t limit, bool fNoZeroPlease) { vector statData(cols.size()); // Loop over all columns in our list of requested columns const size_t last = limit ? first + limit : size_t(-1); while (GetRow(first++)) { const size_t row = GetRow(); if (row==GetNumRows() || row==last) 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(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++) { const double cValue = GetDouble(*it, i); 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 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, size_t first, size_t limit) { // Loop over all columns in our list of requested columns vector> statData; const size_t num = limit==0 || GetNumRows()(it->col.size*num*(it->last-it->first+1))); // Loop over all columns in our list of requested columns const size_t last = limit ? first + limit : size_t(-1); while (GetRow(first++)) { const size_t row = GetRow(); if (row==GetNumRows() || row==last) break; auto statsIt = statData.begin(); for (auto it=cols.begin(); it!=cols.end(); it++, statsIt++) { const char *src = reinterpret_cast(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(*statsIt, fout); break; case 'B': displayStats(*statsIt, fout); break; case 'I': displayStats(*statsIt, fout); break; case 'J': displayStats(*statsIt, fout); break; case 'K': displayStats(*statsIt, fout); break; case 'E': displayStats(*statsIt, fout); break; case 'D': displayStats(*statsIt, fout); break; default: ; } } } // -------------------------------------------------------------------------- // //! Retrieves the configuration parameters //! @param conf //! the configuration object // int FitsDumper::Exec(Configuration& conf) { if (conf.Get("list")) List(); if (conf.Get("header")) ListHeader(conf.Get("outfile")); if (conf.Get("header") || conf.Get("list")) 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; } // ------------------------------------------------------------ const string filename = conf.Get("outfile"); ofstream fout(filename=="-"?"/dev/stdout":filename); if (!fout) { cerr << "Cannot open file " << filename << ": " << strerror(errno) << endl; return false; } fout.precision(conf.Get("precision")); const string filter = conf.Has("filter") ? conf.Get("filter") : ""; const size_t first = conf.Get("first"); const size_t limit = conf.Get("limit"); if (conf.Get("root")) { DumpRoot(fout, conf.Vec("col"), filter, first, limit, filename); return 0; } const vector cols = InitColumns(conf.Vec("col")); if (cols.size()==0) return false; if (conf.Get("minmax")) { DumpMinMax(fout, cols, first, limit, conf.Get("nozero")); return 0; } if (conf.Get("stat")) { DumpStats(fout, cols, first, limit); return 0; } Dump(fout, cols, filter, first, limit, 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" "\n" "Addressing a column:\n" " ColumnName: Will address all fields of a column\n" " ColumnName[n]: Will address the n-th field of a column (starts with 0)\n" " ColumnName[n1:n2]: Will address all fields between n1 and including n2\n" #ifdef HAVE_ROOT "\n" "Selecting a column:\n" " Commandline option: --filter\n" " Explanation: Such a selection is evaluated using TFormula, hence, every " "mathematical operation allowed in TFormula is allowed there, too. " "The reference is the column index as printed in the output stream, " "starting with 1. The index 0 is reserved for the row number.\n" #endif ; cout << endl; } void PrintHelp() { #ifdef HAVE_ROOT cout << "\n\n" "Examples:\n" "In --root mode, fitsdump support TFormula's syntax for all columns and the filter " "You can then refer to a column or a (single) index of the column just by its name " "If the index is omitted, 0 is assumed. Note that the [x:y] syntax in this mode is " "not supported\n" "\n" " fitsdump Zd --filter=\"[0]>20 && cos([1])*TMath::RadToDeg()<45\"\n" "\n" "The columns can also be addressed with their names\n" "\n" " fitsdump -r \"(Zd+Err)*TMath::DegToRad()\" --filter=\"[0]<25 && [1]<0.05\"\n" "\n" "is identical to\n" "\n" " fitsdump -r \"(Zd[0]+Err[0])*TMath::DegToRad()\" --filter=\"[0]<25 && [1]<0.05\"\n" "\n" "A special placeholder exists for the row number\n" "\n" " fitsdump -r \"#\" --filter=\"#>10 && #<100\"\n" "\n"; cout << endl; #endif } void SetupConfiguration(Configuration& conf) { po::options_description configs("Fitsdump options"); configs.add_options() ("fitsfile,f", var() #if BOOST_VERSION >= 104200 ->required() #endif , "Name of FITS file") ("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("-"), "Name of output file (-:/dev/stdout)") ("precision,p", var(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") ("first", var(size_t(0)), "First number of row to read") ("limit", var(size_t(0)), "Limit for the maximum number of rows to read (0=unlimited)") #ifdef HAVE_ROOT ("root,r", po_switch(), "Enable root mode") ("filter,f", var(""), "Filter to restrict the selection of events (e.g. '[0]>10 && [0]<20'; does not work with stat and minmax yet)") #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; 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); }