Index: ct/tools/pyscripts/pyfact/FITS.h =================================================================== --- /fact/tools/pyscripts/pyfact/FITS.h (revision 17730) +++ (revision ) @@ -1,112 +1,0 @@ -/* - * - * FITS.h - * - * Global fits header - * - * Author: lyard - * - */ - -#ifndef MARS_FITS -#define MARS_FITS - -#include - -namespace FITS -{ - - //Identifier of the compression schemes processes - enum CompressionProcess_t - { - kFactRaw = 0x0, - kFactSmoothing = 0x1, - kFactHuffman16 = 0x2 - }; - - //ordering of the columns / rows - enum RowOrdering_t - { - kOrderByCol = 'C', - kOrderByRow = 'R' - }; - -#ifdef __CINT__ - struct TileHeader; - struct BlockHeader; -#else - //Structure helper for tiles headers - struct TileHeader - { - char id[4]; - uint32_t numRows; - uint64_t size; - - TileHeader() {} - - TileHeader(uint32_t nRows, - uint64_t s) : id({'T', 'I', 'L', 'E'}), - numRows(nRows), - size(s) - { }; - } __attribute__((__packed__)); - - //Structure helper for blocks headers and compresion schemes - struct BlockHeader - { - uint64_t size; - char ordering; - unsigned char numProcs; - uint16_t processings[]; - - BlockHeader(uint64_t s=0, - char o=kOrderByRow, - unsigned char n=1) : size(s), - ordering(o), - numProcs(n) - {} - } __attribute__((__packed__)); -#endif - - //Helper structure to simplify the initialization and handling of compressed blocks headers - struct Compression - { - BlockHeader header; - std::vector sequence; - - Compression(const std::vector &seq, const RowOrdering_t &order=kOrderByCol) - : header(0, order, seq.size()), sequence(seq) - { - - } - - Compression(const CompressionProcess_t &compression=kFactRaw, const RowOrdering_t &order=kOrderByCol) - : header(0, order, 1), sequence(1) - { - sequence[0] = compression; - } - -#ifdef __MARS__ // needed for CINT - Compression(const int &compression) - : header(0, kOrderByCol, 1), sequence(1) - { - sequence[0] = compression; - } -#endif - - RowOrdering_t getOrdering() const { return RowOrdering_t(header.ordering); } - uint32_t getSizeOnDisk() const { return sizeof(BlockHeader) + sizeof(uint16_t)*header.numProcs; } - CompressionProcess_t getProc(uint32_t i) const { return CompressionProcess_t(sequence[i]); } - uint16_t getNumProcs() const { return header.numProcs; } - - void SetBlockSize(uint64_t size) { header.size = size; } - void Memcpy(char *dest) const - { - memcpy(dest, &header, sizeof(BlockHeader)); - memcpy(dest+sizeof(BlockHeader), sequence.data(), header.numProcs*sizeof(uint16_t)); - } - }; -}; - -#endif //_FITS_H_ - Index: ct/tools/pyscripts/pyfact/checksum.h =================================================================== --- /fact/tools/pyscripts/pyfact/checksum.h (revision 17730) +++ (revision ) @@ -1,242 +1,0 @@ -#ifndef MARS_checksum -#define MARS_checksum - -#ifndef __CINT__ -#include -#endif - -#include - -class Checksum -{ -public: - uint64_t buffer; - - void reset() { buffer = 0; } - Checksum() : buffer(0) { } - Checksum(const Checksum &sum) : buffer(sum.buffer) { } - Checksum(uint64_t v) : buffer(((v>>16)&0xffff) | ((v&0xffff)<<32)) { } - - uint32_t val() const { return (((buffer&0xffff)<<16) | ((buffer>>32)&0xffff)); } - - bool valid() const { return buffer==0xffff0000ffff; } - - void HandleCarryBits() - { - while (1) - { - const uint64_t carry = ((buffer>>48)&0xffff) | ((buffer&0xffff0000)<<16); - if (!carry) - break; - - buffer = (buffer&0xffff0000ffff) + carry; - } - } - - Checksum &operator+=(const Checksum &sum) - { - buffer += sum.buffer; - HandleCarryBits(); - return *this; - } - - Checksum operator+(Checksum sum) const - { - return (sum += *this); - } - - - bool add(const char *buf, size_t len, bool big_endian = true) - { - // Avoid overflows in carry bits - if (len>262140) // 2^18-4 - { - add(buf, 262140); - return add(buf+262140, len-262140); - } - - if (len%4>0) - { - std::ostringstream sout; - sout << "Length " << len << " not dividable by 4"; - -#ifdef __EXCEPTIONS - throw std::runtime_error(sout.str()); -#else - gLog << ___err___ << "ERROR - " << sout.str() << std::endl; - return false; -#endif - } - - const uint16_t *sbuf = reinterpret_cast(buf); - - uint32_t *hilo = reinterpret_cast(&buffer); - - - const uint16_t *end = sbuf + len/2; - - if (big_endian) - addLoopSwapping(sbuf, end, hilo); - else - addLoop(sbuf, end, hilo); - /*const uint16_t *end = sbuf + len/2; - while (1) - { - if (sbuf==end) - break; - - hilo[0] += ntohs(*sbuf++); - - if (sbuf==end) - break; - - hilo[1] += ntohs(*sbuf++); - }*/ - - HandleCarryBits(); - - return true; - } - - void addLoopSwapping(const uint16_t *sbuf, const uint16_t *end, uint32_t* hilo) - { - /* - for (size_t i = 0; i < len/2; i++) - { - //swap the bytes of the 32 bits value. but... - //the hi and lo values are stored in fits-like order. do not swap them - hilo[i%2] += ntohs(sbuf[i]); //(sbuf[i]&0xff00)>>8 | (sbuf[i]&0x00ff)<<8; - }*/ - - // This is about as twice as fast as the loop above - // ntohs is CPU optimized, i%2 doesn't need to be computed - while (1) - { - if (sbuf==end) - break; - - hilo[0] += ntohs(*sbuf++); - - if (sbuf==end) - break; - - hilo[1] += ntohs(*sbuf++); - } - } - - void addLoop(const uint16_t *sbuf, const uint16_t *end, uint32_t* hilo) - { - while (1) - { - if (sbuf==end) - break; - - hilo[0] += ntohs(*sbuf++); - - if (sbuf==end) - break; - - hilo[1] += ntohs(*sbuf++); - } - } - - bool add(const std::vector &v, bool big_endian = true) - { - return add(v.data(), v.size(), big_endian); - } - - std::string str(bool complm=true) const - { - std::string rc(16,0); - - const uint8_t exclude[13] = - { - 0x3a, 0x3b, 0x3c, 0x3d, 0x3e, 0x3f, 0x40, - 0x5b, 0x5c, 0x5d, 0x5e, 0x5f, 0x60 - }; - - const uint32_t value = complm ? ~val() : val(); // complement each bit of the value - - for (int ii = 0; ii < 4; ii++) - { - uint8_t byte = (value >> (24 - (8 * ii))); - - const uint8_t quotient = byte / 4 + '0'; - const uint8_t remainder = byte % 4; - - uint32_t ch[4] = { static_cast(quotient+remainder), quotient, quotient, quotient }; - - // avoid ASCII punctuation - while (1) - { - bool check = false; - for (int kk = 0; kk < 13; kk++) - { - for (int jj = 0; jj < 4; jj += 2) - { - if (ch[jj] != exclude[kk] && ch[jj+1] != exclude[kk]) - continue; - - ch[jj]++; - ch[jj+1]--; - check++; - } - } - - if (!check) - break; - } - - for (int jj = 0; jj < 4; jj++) // assign the bytes - rc[4*jj+ii] = ch[jj]; - } - - const char lastChar = rc[15]; - for (int i=15;i>0;i--) - rc[i] = rc[i-1]; - rc[0] = lastChar; - return rc; -/* - uint8_t *p = reinterpret_cast(&value); - //swap the bytes of the value - uint8_t temp; - temp = p[0]; - p[0] = p[3]; - p[3] = temp; - temp = p[1]; - p[1] = p[2]; - p[2] = temp; - - for (int i=0; i<4; i++) - { - rc[i+ 0] = '0' + p[i]/4 + p[i]%4; - rc[i+ 4] = '0' + p[i]/4; - rc[i+ 8] = '0' + p[i]/4; - rc[i+12] = '0' + p[i]/4; - } - - while(1) - { - bool ok = true; - for (int i=0; i<16-4; i++) - { - for (int j=0; j<13; j++) - if (rc[i]==exclude[j] || rc[(i+4)%16]==exclude[j]) - { - rc[i]++; - rc[(i+4)%16]--; - ok = false; - } - } - - if (ok) - break; - } - - - return rc; - */ - } -}; - -#endif Index: ct/tools/pyscripts/pyfact/factfits.h =================================================================== --- /fact/tools/pyscripts/pyfact/factfits.h (revision 17730) +++ (revision ) @@ -1,206 +1,0 @@ -/* - * factfits.h - * - * Created on: May 26, 2013 - * Author: lyard - */ - -#ifndef MARS_FACTFITS -#define MARS_FACTFITS - -#include "zfits.h" - -class factfits : public zfits -{ -public: - // Default constructor - factfits(const std::string& fname, const std::string& tableName="", bool force=false) : - zfits(fname, tableName, force), - fOffsetCalibration(0), - fOffsetStartCellData(0), - fOffsetData(0), - fNumRoi(0) - { - if (init()) - readDrsCalib(fname); - } - - // Alternative constructor - factfits(const std::string& fname, const std::string& fout, const std::string& tableName, bool force=false) : - zfits(fname, fout, tableName, force), - fOffsetCalibration(0), - fOffsetStartCellData(0), - fOffsetData(0), - fNumRoi(0) - { - if (init()) - readDrsCalib(fname); - } - -private: - - void StageRow(size_t row, char* dest) - { - zfits::StageRow(row, dest); - - // This file does not contain fact data or no calibration to be applied - if (fOffsetCalibration.empty()) - return; - - //re-get the pointer to the data to access the offsets - const uint8_t offset = (row*fTable.bytes_per_row)%4; - - int16_t *startCell = reinterpret_cast(fBufferRow.data() + offset + fOffsetStartCellData); - int16_t *data = reinterpret_cast(fBufferRow.data() + offset + fOffsetData); - - /* - for (uint32_t i=0; i<1440*1024; i+=1024, startCell++) - { - if (*startCell < 0) - { - data += fNumRoi; - continue; - } - for (uint32_t j=0; j1024) - { - while (cal("NPIX")!=1440) - return false; - - fNumRoi = Get("NROI"); - if (fNumRoi>1024) - return false; - - // check column details for Data - const Table::Columns::const_iterator it = fTable.cols.find("Data"); - if (it==fTable.cols.end() || it->second.num!=1440*fNumRoi || it->second.type!='I') - return false; - - // check column details for StartCellData - const Table::Columns::const_iterator is = fTable.cols.find("StartCellData"); - if (is==fTable.cols.end() || is->second.num!=1440 || is->second.type!='I') - return false; - - fOffsetStartCellData = is->second.offset; - fOffsetData = it->second.offset; - - return true; - } - - // Read the Drs calibration data - void readDrsCalib(const std::string& fileName) - { - //should not be mandatory, but improves the perfs a lot when reading not compressed, gzipped files - if (!IsCompressedFITS()) - return; - - zfits calib(fileName, "ZDrsCellOffsets"); - - if (calib.bad()) - { - clear(rdstate()|std::ios::badbit); - return; - } - - if (calib.eof()) - return; - - // Check correct size and format of table - if (calib.GetNumRows() != 1) - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("Table 'ZDrsCellOffsets' found, but not with one row as expected"); -#else - gLog << ___err___ << "ERROR - Table 'ZDrsCellOffsets' found, but not with one row as expected" << std::endl; - return; -#endif - } - if (calib.GetStr("TTYPE1") != "OffsetCalibration") - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("Table 'ZDrsCellOffsets' found, but first column is not the one expected"); -#else - gLog << ___err___ << "ERROR - Table 'ZDrsCellOffsets' found, but first column is not the one expected" << std::endl; - return; -#endif - } - bool isColumnPresent = false; - if (calib.HasKey("TFORM1") && calib.GetStr("TFORM1") == "1474560I") isColumnPresent = true; - if (calib.HasKey("ZFORM1") && calib.GetStr("ZFORM1") == "1474560I") isColumnPresent = true; - if (!isColumnPresent) // 1024*1440 - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("Table 'ZDrsCellOffsets' has wrong column format (TFROM1)"); -#else - gLog << ___err___ << "ERROR - Table 'ZDrsCellOffsets' has wrong column format (TFORM1)" << std::endl; - return; -#endif - } - - fOffsetCalibration.resize(1024*1440); - - calib.SetPtrAddress("OffsetCalibration", fOffsetCalibration.data()); - if (calib.GetNextRow()) - return; - - clear(rdstate()|std::ios::badbit); - -#ifdef __EXCEPTIONS - throw std::runtime_error("Reading column 'OffsetCalibration' failed."); -#else - gLog << ___err___ << "ERROR - Reading column 'OffsetCalibration' failed." << std::endl; -#endif - - } - - std::vector fOffsetCalibration; ///< integer values of the drs calibration used for compression - - size_t fOffsetStartCellData; - size_t fOffsetData; - - uint16_t fNumRoi; - -//#warning Time marker channels currently unhandled - -}; //class factfits - -#endif Index: ct/tools/pyscripts/pyfact/fits.h =================================================================== --- /fact/tools/pyscripts/pyfact/fits.h (revision 17730) +++ (revision ) @@ -1,1120 +1,0 @@ -#ifndef MARS_fits -#define MARS_fits - -#include - -#include -#include -#include -#include -#include -#include - -#define GCC_VERSION (__GNUC__ * 10000 + __GNUC_MINOR__ * 100 + __GNUC_PATCHLEVEL__) - -#ifndef __CINT__ -#include -#else -#define off_t size_t -namespace std -{ - template class unordered_map; -} -#endif - -#ifndef __MARS__ -#include -#include -#include -#ifndef gLog -#define gLog std::cerr -#define ___err___ "" -#define ___warn___ "" -#define ___all___ "" -#endif -#else -#include "MLog.h" -#include "MLogManip.h" -#define ___err___ err -#define ___warn___ warn -#define ___all___ all -#endif - -#define HAVE_ZLIB -#if defined(HAVE_ZLIB) || defined(__CINT__) -#include "izstream.h" -#else -#include -#define izstream ifstream -#warning Support for zipped FITS files disabled. -#endif - -#include "checksum.h" - -class fits : public izstream -{ -public: - //I know I know, you're going to yiell that this does not belong here. - //It will belong in the global scope eventually, and it makes the coding of zfits much simpler this way. - enum Compression_t - { - kCompUnknown, - kCompFACT - }; - - struct Entry - { - char type; - std::string value; - std::string comment; - std::string fitsString; - - template - T Get() const - { - T t; - - std::istringstream str(value); - str >> t; - - return t; - } - }; - - struct Table - { - off_t offset; - - bool is_compressed; - - std::string name; - size_t bytes_per_row; - size_t num_rows; - size_t num_cols; - size_t total_bytes; // NAXIS1*NAXIS2 - - struct Column - { - size_t offset; - size_t num; - size_t size; - size_t bytes; // num*size - char type; - std::string unit; - Compression_t comp; - }; - - typedef std::map Keys; - typedef std::map Columns; - typedef std::vector SortedColumns; - - Columns cols; - SortedColumns sorted_cols; - Keys keys; - - //------------ vectors containing the same info as cols and keys------ - // They are filled at the end of the contstuctor - // They are used py python, for retrieving the contents of the maps - // I did not find out how to return a map of via ROOT to python - // So I use this ugly hack. I guess its no problem, its just ugly - // - // btw: The way these variables are named, is imho as ugly as possible, - // but I will stick to it. - // usually a std::map contains key-value pairs - // that means the e.g. map contains keys of type string - // and values of type Entry. - // but now the map itself was called 'keys', which means one should say: - // the keys of 'keys' are of type string - // and - // the value of 'keys' are of type Entry. - // not to forget that: - // one field of a value of 'keys' is named 'value' - - // vectors for Keys keys; - vector Py_KeyKeys; - vector Py_KeyValues; - vector Py_KeyComments; - vector Py_KeyTypes; - - // vectors for Columns cols; - vector Py_ColumnKeys; - vector Py_ColumnOffsets; - vector Py_ColumnNums; - vector Py_ColumnSizes; - vector Py_ColumnTypes; - vector Py_ColumnUnits; - //-end of----- vectors containing the same info as cols and keys------ - - - int64_t datasum; - - std::string Trim(const std::string &str, char c=' ') const - { - // Trim Both leading and trailing spaces - const size_t pstart = str.find_first_not_of(c); // Find the first character position after excluding leading blank spaces - const size_t pend = str.find_last_not_of(c); // Find the first character position from reverse af - - // if all spaces or empty return an empty string - if (std::string::npos==pstart || std::string::npos==pend) - return std::string(); - - return str.substr(pstart, pend-pstart+1); - } - - bool Check(const std::string &key, char type, const std::string &value="") const - { - const Keys::const_iterator it = keys.find(key); - if (it==keys.end()) - { - std::ostringstream str; - str << "Key '" << key << "' not found."; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return false; -#endif - } - - if (it->second.type!=type) - { - std::ostringstream str; - str << "Wrong type for key '" << key << "': expected " << type << ", found " << it->second.type << "."; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return false; -#endif - } - - if (!value.empty() && it->second.value!=value) - { - std::ostringstream str; - str << "Wrong value for key '" << key << "': expected " << value << ", found " << it->second.value << "."; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return false; -#endif - } - - return true; - } - - Keys ParseBlock(const std::vector &vec) const - { - std::map rc; - - for (unsigned int i=0; i &vec, off_t off) : offset(off), - keys(ParseBlock(vec)) - { - is_compressed = HasKey("ZTABLE") && Check("ZTABLE", 'B', "T"); - - if (!Check("XTENSION", 'T', "BINTABLE") || - !Check("NAXIS", 'I', "2") || - !Check("BITPIX", 'I', "8") || - !Check("GCOUNT", 'I', "1") || - !Check("EXTNAME", 'T') || - !Check("NAXIS1", 'I') || - !Check("NAXIS2", 'I') || - !Check("TFIELDS", 'I')) - return; - - if (is_compressed) - { - if (!Check("ZNAXIS1", 'I') || - !Check("ZNAXIS2", 'I') || - !Check("ZPCOUNT", 'I', "0")) - return; - } - else - { - if (!Check("PCOUNT", 'I', "0")) - return; - } - - total_bytes = Get("NAXIS1")*Get("NAXIS2"); - bytes_per_row = is_compressed ? Get("ZNAXIS1") : Get("NAXIS1"); - num_rows = is_compressed ? Get("ZNAXIS2") : Get("NAXIS2"); - num_cols = Get("TFIELDS"); - datasum = is_compressed ? Get("DATASUM", -1) : Get("DATASUM", -1); -//cout << "IS COMPRESSED =-========= " << is_compressed << " " << Get("NAXIS1") << endl; - size_t bytes = 0; - - const std::string tFormName = is_compressed ? "ZFORM" : "TFORM"; - for (long long unsigned int i=1; i<=num_cols; i++) - { - const std::string num(std::to_string(i)); - - if (!Check("TTYPE"+num, 'T') || - !Check(tFormName+num, 'T')) - return; - - const std::string id = Get("TTYPE"+num); - const std::string fmt = Get(tFormName+num); - const std::string unit = Get("TUNIT"+num, ""); - const std::string comp = Get("ZCTYP"+num, ""); - - Compression_t compress = kCompUnknown; - if (comp == "FACT") - compress = kCompFACT; - - std::istringstream sin(fmt); - int n = 0; - sin >> n; - if (!sin) - n = 1; - - const char type = fmt[fmt.length()-1]; - - size_t size = 0; - switch (type) - { - // We could use negative values to mark floats - // otheriwse we could just cast them to int64_t? - case 'L': // logical - case 'A': // char - case 'B': size = 1; break; // byte - case 'I': size = 2; break; // short - case 'J': size = 4; break; // int - case 'K': size = 8; break; // long long - case 'E': size = 4; break; // float - case 'D': size = 8; break; // double - // case 'X': size = n; break; // bits (n=number of bytes needed to contain all bits) - // case 'C': size = 8; break; // complex float - // case 'M': size = 16; break; // complex double - // case 'P': size = 8; break; // array descriptor (32bit) - // case 'Q': size = 16; break; // array descriptor (64bit) - default: - { - std::ostringstream str; - str << "FITS format TFORM='" << fmt << "' not yet supported."; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return; -#endif - } - } - - const Table::Column col = { bytes, static_cast(n), size, n*size, type, unit, compress}; - - cols[id] = col; - sorted_cols.push_back(col); - bytes += n*size; - } - - if (bytes!=bytes_per_row) - { -#ifdef __EXCEPTIONS - throw std::runtime_error("Column size mismatch"); -#else - gLog << ___err___ << "ERROR - Column size mismatch" << std::endl; - return; -#endif - } - - name = Get("EXTNAME"); - - Fill_Py_Vectors(); - } - - void Fill_Py_Vectors(void) - { - // Fill the Py_ vectors see line: 90ff - // vectors for Keys keys; - for (Keys::const_iterator it=keys.begin(); it!=keys.end(); it++) - { - Py_KeyKeys.push_back(it->first); - Py_KeyValues.push_back(it->second.value); - Py_KeyComments.push_back(it->second.comment); - Py_KeyTypes.push_back(it->second.type); - } - - // vectors for Columns cols; - for (Columns::const_iterator it=cols.begin(); it!=cols.end(); it++) - { - Py_ColumnKeys.push_back(it->first); - Py_ColumnTypes.push_back(it->second.type); - Py_ColumnOffsets.push_back(it->second.offset); - Py_ColumnNums.push_back(it->second.num); - Py_ColumnSizes.push_back(it->second.size); - Py_ColumnUnits.push_back(it->second.unit); - } - } - - void PrintKeys(bool display_all=false) const - { - for (Keys::const_iterator it=keys.cbegin(); it!=keys.cend(); it++) - { - if (!display_all && - (it->first.substr(0, 6)=="TTYPE" || - it->first.substr(0, 6)=="TFORM" || - it->first.substr(0, 6)=="TUNIT" || - it->first=="TFIELDS" || - it->first=="XTENSION" || - it->first=="NAXIS" || - it->first=="BITPIX" || - it->first=="PCOUNT" || - it->first=="GCOUNT") - ) - continue; - - gLog << ___all___ << std::setw(2) << it->second.type << '|' << it->first << '=' << it->second.value << '/' << it->second.comment << '|' << std::endl; - }} - - void PrintColumns() const - { - typedef std::map, Column> Sorted; - - Sorted sorted; - - for (Columns::const_iterator it=cols.cbegin(); it!=cols.cend(); it++) - sorted[std::make_pair(it->second.offset, it->first)] = it->second; - - for (Sorted::const_iterator it=sorted.cbegin(); it!=sorted.cend(); it++) - { - gLog << ___all___ << std::setw(6) << it->second.offset << "| "; - gLog << it->second.num << 'x'; - switch (it->second.type) - { - case 'A': gLog << "char(8)"; break; - case 'L': gLog << "bool(8)"; break; - case 'B': gLog << "byte(8)"; break; - case 'I': gLog << "short(16)"; break; - case 'J': gLog << "int(32)"; break; - case 'K': gLog << "int(64)"; break; - case 'E': gLog << "float(32)"; break; - case 'D': gLog << "double(64)"; break; - } - gLog << ": " << it->first.second << " [" << it->second.unit << "]" << std::endl; - } - } - - operator bool() const { return !name.empty(); } - - bool HasKey(const std::string &key) const - { - return keys.find(key)!=keys.end(); - } - - bool HasColumn(const std::string& col) const - { - return cols.find(col)!=cols.end(); - } - - const Columns &GetColumns() const - { - return cols; - } - - const Keys &GetKeys() const - { - return keys; - } - - // Values of keys are always signed - template - T Get(const std::string &key) const - { - const std::map::const_iterator it = keys.find(key); - if (it==keys.end()) - { - std::ostringstream str; - str << "Key '" << key << "' not found."; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return T(); -#endif - } - return it->second.Get(); - } - - // Values of keys are always signed - template - T Get(const std::string &key, const T &deflt) const - { - const std::map::const_iterator it = keys.find(key); - return it==keys.end() ? deflt :it->second.Get(); - } - - size_t GetN(const std::string &key) const - { - const Columns::const_iterator it = cols.find(key); - return it==cols.end() ? 0 : it->second.num; - } - - - - // There may be a gap between the main table and the start of the heap: - // this computes the offset - streamoff GetHeapShift() const - { - if (!HasKey("ZHEAPPTR")) - return 0; - - const size_t shift = Get("ZHEAPPTR"); - return shift <= total_bytes ? 0 : shift - total_bytes; - } - - // return total number of bytes 'all inclusive' - streamoff GetTotalBytes() const - { - //get offset of special data area from start of main table - const streamoff shift = GetHeapShift(); - - //and special data area size - const streamoff size = HasKey("PCOUNT") ? Get("PCOUNT") : 0; - - // Get the total size - const streamoff total = total_bytes + size + shift; - - // check for padding - if (total%2880==0) - return total; - - // padding necessary - return total + (2880 - (total%2880)); - } - }; - -protected: - std::ofstream fCopy; - std::vector fListOfTables; // List of skipped tables. Last table is open table - - Table fTable; - - typedef std::pair Address; - typedef std::vector
Addresses; - //map fAddresses; - Addresses fAddresses; - - typedef std::unordered_map Pointers; - Pointers fPointers; - - std::vector> fGarbage; - - std::vector fBufferRow; - std::vector fBufferDat; - - size_t fRow; - - Checksum fChkHeader; - Checksum fChkData; - - bool ReadBlock(std::vector &vec) - { - int endtag = 0; - for (int i=0; i<36; i++) - { - char c[81]; - c[80] = 0; - read(c, 80); - if (!good()) - break; - - fChkHeader.add(c, 80); - -// if (c[0]==0) -// return vector(); - - std::string str(c); - -// if (!str.empty()) -// cout << setw(2) << i << "|" << str << "|" << (endtag?'-':'+') << endl; - - if (endtag==2 || str=="END ") - { - endtag = 2; // valid END tag found - continue; - } - - if (endtag==1 || str==" ") - { - endtag = 1; // end tag not found, but expected to be there - continue; - } - - vec.push_back(str); - } - - // Make sure that no empty vector is returned - if (endtag && vec.size()%36==0) - vec.emplace_back("END = '' / "); - - return endtag==2; - } - - std::string Compile(const std::string &key, int16_t i=-1) const - { -#if GCC_VERSION < 40603 - return i<0 ? key : key+std::to_string((long long int)(i)); -#else - return i<0 ? key : key+std::to_string(i); -#endif - } - - void Constructor(const std::string &fname, std::string fout="", const std::string& tableName="", bool force=false) - { - char simple[10]; - read(simple, 10); - if (!good()) - return; - - if (memcmp(simple, "SIMPLE = ", 10)) - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("File is not a FITS file."); -#else - gLog << ___err___ << "ERROR - File is not a FITS file." << std::endl; - return; -#endif - } - - seekg(0); - - while (good()) - { - std::vector block; - while (1) - { - // If we search for a table, we implicitly assume that - // not finding the table is not an error. The user - // can easily check that by eof() && !bad() - peek(); - if (eof() && !bad() && !tableName.empty()) - { - break; - } - // FIXME: Set limit on memory consumption - const int rc = ReadBlock(block); - if (!good()) - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("FITS file corrupted."); -#else - gLog << ___err___ << "ERROR - FITS file corrupted." << std::endl; - return; -#endif - } - - if (block.size()%36) - { - if (!rc && !force) - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("END keyword missing in FITS header."); -#else - gLog << ___err___ << "ERROR - END keyword missing in FITS file... file might be corrupted." << std::endl; - return; -#endif - } - break; - } - } - - if (block.empty()) - break; - - if (block[0].substr(0, 9)=="SIMPLE =") - { - fChkHeader.reset(); - continue; - } - - if (block[0].substr(0, 9)=="XTENSION=") - { - fTable = Table(block, tellg()); - fRow = (size_t)-1; - - if (!fTable) - { - clear(rdstate()|std::ios::badbit); - return; - } - - const std::string &tname = fTable.Get("EXTNAME"); - - fListOfTables.emplace_back(tname); - - // Check for table name. Skip until eof or requested table are found. - // skip the current table? - if ((!tableName.empty() && tableName!=tname) || (tableName.empty() && "ZDrsCellOffsets"==tname)) - { - const streamoff skip = fTable.GetTotalBytes(); - seekg(skip, std::ios_base::cur); - - fChkHeader.reset(); - - continue; - } - - fBufferRow.resize(fTable.bytes_per_row + 8-fTable.bytes_per_row%4); - fBufferDat.resize(fTable.bytes_per_row); - - break; - } - } - - if (fout.empty()) - return; - - if (*fout.rbegin()=='/') - { - const size_t p = fname.find_last_of('/'); - fout.append(fname.substr(p+1)); - } - - fCopy.open(fout); - if (!fCopy) - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("Could not open output file."); -#else - gLog << ___err___ << "ERROR - Failed to open output file." << std::endl; - return; -#endif - } - - const streampos p = tellg(); - seekg(0); - - std::vector buf(p); - read(buf.data(), p); - - fCopy.write(buf.data(), p); - if (!fCopy) - clear(rdstate()|std::ios::badbit); - } - -public: - fits(const std::string &fname, const std::string& tableName="", bool force=false) : izstream(fname.c_str()) - { - Constructor(fname, "", tableName, force); - if ((fTable.is_compressed ||fTable.name=="ZDrsCellOffsets") && !force) - { -#ifdef __EXCEPTIONS - throw std::runtime_error("You are trying to read a compressed fits with the base fits class. Please use factfits instead."); -#else - gLog << ___err___ << "ERROR - You are trying to read a compressed fits with the base fits class. Please use factfits instead." << std::endl; -#endif - clear(rdstate()|std::ios::badbit); - } - } - - fits(const std::string &fname, const std::string &fout, const std::string& tableName, bool force=false) : izstream(fname.c_str()) - { - Constructor(fname, fout, tableName, force); - if ((fTable.is_compressed || fTable.name=="ZDrsCellOffsets") && !force) - { -#ifdef __EXCEPTIONS - throw std::runtime_error("You are trying to read a compressed fits with the base fits class. Please use factfits instead."); -#else - gLog << ___err___ << "ERROR - You are trying to read a compressed fits with the base fits class. Please use factfits instead." << std::endl; -#endif - clear(rdstate()|std::ios::badbit); - } - } - - fits() : izstream() - { - - } - - ~fits() - { - std::copy(std::istreambuf_iterator(*this), - std::istreambuf_iterator(), - std::ostreambuf_iterator(fCopy)); - } - - virtual void StageRow(size_t row, char* dest) - { - // if (row!=fRow+1) // Fast seeking is ensured by izstream - seekg(fTable.offset+row*fTable.bytes_per_row); - read(dest, fTable.bytes_per_row); - //fin.clear(fin.rdstate()&~ios::eofbit); - } - - virtual void WriteRowToCopyFile(size_t row) - { - if (row==fRow+1) - { - const uint8_t offset = (row*fTable.bytes_per_row)%4; - - fChkData.add(fBufferRow); - if (fCopy.is_open() && fCopy.good()) - fCopy.write(fBufferRow.data()+offset, fTable.bytes_per_row); - if (!fCopy) - clear(rdstate()|std::ios::badbit); - } - else - if (fCopy.is_open()) - clear(rdstate()|std::ios::badbit); - } - - void ZeroBufferForChecksum(std::vector& vec, const uint64_t extraZeros=0) - { - auto ib = vec.begin(); - auto ie = vec.end(); - - *ib++ = 0; - *ib++ = 0; - *ib++ = 0; - *ib = 0; - - for (uint64_t i=0;i - void revcpy(char *dest, const char *src, const int &num) - { - const char *pend = src + num*N; - for (const char *ptr = src; ptr(dest, src, c.num); break; - case 4: revcpy<4>(dest, src, c.num); break; - case 8: revcpy<8>(dest, src, c.num); break; - } - } - - virtual bool GetRow(size_t row, bool check=true) - { - if (check && row>=fTable.num_rows) - return false; - - const uint8_t offset = ReadRow(row); - if (!good()) - return good(); - - const char *ptr = fBufferRow.data() + offset; - - for (Addresses::const_iterator it=fAddresses.cbegin(); it!=fAddresses.cend(); it++) - { - const Table::Column &c = it->second; - - const char *src = ptr + c.offset; - char *dest = reinterpret_cast(it->first); - - MoveColumnDataToUserSpace(dest, src, c); - } - - return good(); - } - - bool GetNextRow(bool check=true) - { - return GetRow(fRow+1, check); - } - - virtual bool SkipNextRow() - { - seekg(fTable.offset+(++fRow)*fTable.bytes_per_row); - return good(); - } - - static bool Compare(const Address &p1, const Address &p2) - { - return p1.first>p2.first; - } - - template - const T &GetAs(const std::string &name) - { - return *reinterpret_cast(fPointers[name]); - } - - void *SetPtrAddress(const std::string &name) - { - if (fTable.cols.count(name)==0) - { - std::ostringstream str; - str << "SetPtrAddress('" << name << "') - Column not found."; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return NULL; -#endif - } - - Pointers::const_iterator it = fPointers.find(name); - if (it!=fPointers.end()) - return it->second; - - fGarbage.emplace_back(fTable.cols[name].bytes); - - void *ptr = fGarbage.back().data(); - - fPointers[name] = ptr; - fAddresses.emplace_back(ptr, fTable.cols[name]); - sort(fAddresses.begin(), fAddresses.end(), Compare); - return ptr; - } - - template - bool SetPtrAddress(const std::string &name, T *ptr, size_t cnt) - { - if (fTable.cols.count(name)==0) - { - std::ostringstream str; - str << "SetPtrAddress('" << name << "') - Column not found."; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return false; -#endif - } - - if (sizeof(T)!=fTable.cols[name].size) - { - std::ostringstream str; - str << "SetPtrAddress('" << name << "') - Element size mismatch: expected " - << fTable.cols[name].size << " from header, got " << sizeof(T); -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return false; -#endif - } - - if (cnt!=fTable.cols[name].num) - { - std::ostringstream str; - str << "SetPtrAddress('" << name << "') - Element count mismatch: expected " - << fTable.cols[name].num << " from header, got " << cnt; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return false; -#endif - } - - // if (fAddresses.count(ptr)>0) - // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl; - - //fAddresses[ptr] = fTable.cols[name]; - fPointers[name] = ptr; - fAddresses.emplace_back(ptr, fTable.cols[name]); - sort(fAddresses.begin(), fAddresses.end(), Compare); - return true; - } - - template - bool SetRefAddress(const std::string &name, T &ptr) - { - return SetPtrAddress(name, &ptr, sizeof(ptr)/sizeof(T)); - } - - template - bool SetVecAddress(const std::string &name, std::vector &vec) - { - return SetPtrAddress(name, vec.data(), vec.size()); - } - - template - T Get(const std::string &key) const - { - return fTable.Get(key); - } - - template - T Get(const std::string &key, const std::string &deflt) const - { - return fTable.Get(key, deflt); - } - - bool SetPtrAddress(const std::string &name, void *ptr, size_t cnt=0) - { - if (fTable.cols.count(name)==0) - { - std::ostringstream str; - str <<"SetPtrAddress('" << name << "') - Column not found."; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return false; -#endif - } - - if (cnt && cnt!=fTable.cols[name].num) - { - std::ostringstream str; - str << "SetPtrAddress('" << name << "') - Element count mismatch: expected " - << fTable.cols[name].num << " from header, got " << cnt; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - " << str.str() << std::endl; - return false; -#endif - } - - // if (fAddresses.count(ptr)>0) - // gLog << warn << "SetPtrAddress('" << name << "') - Pointer " << ptr << " already assigned." << endl; - - //fAddresses[ptr] = fTable.cols[name]; - fPointers[name] = ptr; - fAddresses.emplace_back(ptr, fTable.cols[name]); - sort(fAddresses.begin(), fAddresses.end(), Compare); - return true; - } - - bool HasKey(const std::string &key) const { return fTable.HasKey(key); } - bool HasColumn(const std::string& col) const { return fTable.HasColumn(col);} - const Table::Columns &GetColumns() const { return fTable.GetColumns();} - const Table::SortedColumns& GetSortedColumns() const { return fTable.sorted_cols;} - const Table::Keys &GetKeys() const { return fTable.GetKeys();} - - int64_t GetInt(const std::string &key) const { return fTable.Get(key); } - uint64_t GetUInt(const std::string &key) const { return fTable.Get(key); } - double GetFloat(const std::string &key) const { return fTable.Get(key); } - std::string GetStr(const std::string &key) const { return fTable.Get(key); } - - size_t GetN(const std::string &key) const - { - return fTable.GetN(key); - } - -// size_t GetNumRows() const { return fTable.num_rows; } - size_t GetRow() const { return fRow==(size_t)-1 ? 0 : fRow; } - - operator bool() const { return fTable && fTable.offset!=0; } - - void PrintKeys() const { fTable.PrintKeys(); } - void PrintColumns() const { fTable.PrintColumns(); } - - // 'Wrappers' for the Table's Getters of the Py Vectors - vector GetPy_KeyKeys() {return fTable.Py_KeyKeys; } - vector GetPy_KeyValues() {return fTable.Py_KeyValues; } - vector GetPy_KeyComments() {return fTable.Py_KeyComments; } - vector GetPy_KeyTypes() {return fTable.Py_KeyTypes; } - - vector GetPy_ColumnKeys() {return fTable.Py_ColumnKeys; } - vector GetPy_ColumnOffsets() {return fTable.Py_ColumnOffsets; } - vector GetPy_ColumnNums() {return fTable.Py_ColumnNums; } - vector GetPy_ColumnSizes() {return fTable.Py_ColumnSizes;} - vector GetPy_ColumnTypes() {return fTable.Py_ColumnTypes;} - vector GetPy_ColumnUnits() {return fTable.Py_ColumnUnits;} - - bool IsHeaderOk() const { return fTable.datasum<0?false:(fChkHeader+Checksum(fTable.datasum)).valid(); } - virtual bool IsFileOk() const { return (fChkHeader+fChkData).valid(); } - - bool IsCompressedFITS() const { return fTable.is_compressed;} - - virtual size_t GetNumRows() const - { - return fTable.Get("NAXIS2"); - } - - virtual size_t GetBytesPerRow() const - { - return fTable.Get("NAXIS1"); - } - - const std::vector &GetTables() const - { - return fListOfTables; - } -}; - -#endif Index: ct/tools/pyscripts/pyfact/huffman.h =================================================================== --- /fact/tools/pyscripts/pyfact/huffman.h (revision 17730) +++ (revision ) @@ -1,435 +1,0 @@ -#ifndef FACT_huffman -#define FACT_huffman - -#include -#include - -#include -#include -#include - -#define MAX_SYMBOLS (1<<16) - -// ================================================================ - -namespace Huffman -{ - static unsigned long numbytes_from_numbits(unsigned long numbits) - { - return numbits / 8 + (numbits % 8 ? 1 : 0); - } - - struct TreeNode - { - TreeNode *parent; - union - { - struct - { - TreeNode *zero, *one; - }; - uint16_t symbol; - }; - - size_t count; - bool isLeaf; - - TreeNode(uint16_t sym, size_t cnt=0) : parent(0), isLeaf(true) - { - symbol = sym; - count = cnt; - } - - TreeNode(TreeNode *n0=0, TreeNode *n1=0) : parent(0), isLeaf(false) - { - count = n0 && n1 ? n0->count + n1->count : 0; - zero = n0 && n1 ? (n0->count > n1->count ? n0 : n1) : NULL; - one = n0 && n1 ? (n0->count > n1->count ? n1 : n0) : NULL; - - if (n0) - n0->parent = this; - - if (n1) - n1->parent = this; - } - - ~TreeNode() - { - if (isLeaf) - return; - - if (zero) - delete zero; - if (one) - delete one; - } - - bool operator() (const TreeNode *hn1, const TreeNode *hn2) const - { - return hn1->count < hn2->count; - } - }; - - - struct Encoder - { - struct Code - { - size_t bits; - uint8_t numbits; - - Code() : numbits(0) { } - }; - - size_t count; - Code lut[1<<16]; - - bool CreateEncoder(const TreeNode *n, size_t bits=0, uint8_t nbits=0) - { - if (n->isLeaf) - { -#ifdef __EXCEPTIONS - if (nbits>sizeof(size_t)*8) - throw std::runtime_error("Too many different symbols - this should not happen!"); -#else - if (nbits>sizeof(size_t)*8) - { - count = 0; - return false; - } -#endif - lut[n->symbol].bits = bits; - lut[n->symbol].numbits = nbits==0 ? 1 : nbits; - count++; - return true; - } - - return - CreateEncoder(n->zero, bits, nbits+1) && - CreateEncoder(n->one, bits | (1<numbits; - const uint8_t *bits = (uint8_t*)&code->bits; - - while (nbits>0) - { - // Number of bits available in the current byte - const uint8_t free_bits = 8 - curbit; - - // Write bits to current byte - curbyte |= *bits<=free_bits) - { - out += curbyte; - curbyte = *bits>>free_bits; - - bits++; - } - - // Adapt the number of available bits, the number of consumed bits - // and the bit-pointer accordingly - const uint8_t consumed = nbits>8 ? 8 : nbits; - nbits -= consumed; - curbit += consumed; - curbit %= 8; - } - } - - // If the buffer-byte is half-full, also add it to the output buffer - if (curbit>0) - out += curbyte; - } - - Encoder(const uint16_t *bufin, size_t bufinlen) : count(0) - { - uint64_t counts[MAX_SYMBOLS]; - memset(counts, 0, sizeof(uint64_t)*MAX_SYMBOLS); - - // Count occurances - for (const uint16_t *p=bufin; p set; - for (int i=0; i1) - { - auto it = set.begin(); - - auto it1 = it++; - auto it2 = it; - - TreeNode *nn = new TreeNode(*it1, *it2); - - set.erase(it1, ++it2); - - set.insert(nn); - } - - // get the root of the tree - const TreeNode *root = *set.begin(); - - CreateEncoder(root); - - // This will delete the whole tree - delete root; - } - - }; - - - - struct Decoder - { - uint16_t symbol; - uint8_t nbits; - bool isLeaf; - - Decoder *lut; - - Decoder() : isLeaf(false), lut(NULL) - { - } - - ~Decoder() - { - if (lut) - delete [] lut; - } - - void Set(uint16_t sym, uint8_t n=0, size_t bits=0) - { - if (!lut) - lut = new Decoder[256]; - - if (n>8) - { - lut[bits&0xff].Set(sym, n-8, bits>>8); - return; - } - - const int nn = 1<<(8-n); - - for (int i=0; ilut->symbol; - return in_ptr; - } - - uint8_t curbit = 0; - while (in_ptr> curbit); - -#ifdef __EXCEPTIONS - if (!p->lut) - throw std::runtime_error("Unknown bitcode in stream!"); -#else - if (!p->lut) - return NULL; -#endif - - p = p->lut + curbyte; - if (!p->isLeaf) - { - in_ptr++; - continue; - } - - *out_ptr++ = p->symbol; - curbit += p->nbits; - - p = this; - - if (curbit>=8) - { - curbit %= 8; - in_ptr++; - } - - } - - return curbit ? in_ptr+1 : in_ptr; - } - - Decoder(const uint8_t* bufin, int64_t &pindex) : isLeaf(false), lut(NULL) - { - // FIXME: Sanity check for size missing.... - - // Read the number of entries. - size_t count=0; - memcpy(&count, bufin + pindex, sizeof(count)); - pindex += sizeof(count); - - // Read the entries. - for (size_t i=0; isizeof(size_t)) - throw std::runtime_error("Number of bytes for a single symbol exceeds maximum."); -#else - if (numbytes>sizeof(size_t)) - { - pindex = -1; - return; - } -#endif - size_t bits=0; - memcpy(&bits, bufin+pindex, numbytes); - pindex += numbytes; - - Set(sym, numbits, bits); - } - } - }; - - inline bool Encode(std::string &bufout, const uint16_t *bufin, size_t bufinlen) - { - const Encoder encoder(bufin, bufinlen); - -#ifndef __EXCEPTIONS - if (encoder.count==0) - return false; -#endif - - bufout.append((char*)&bufinlen, sizeof(size_t)); - encoder.WriteCodeTable(bufout); - encoder.Encode(bufout, bufin, bufinlen); - - return true; - } - - inline int64_t Decode(const uint8_t *bufin, - size_t bufinlen, - std::vector &pbufout) - { - int64_t i = 0; - - // Read the number of data bytes this encoding represents. - size_t data_count = 0; - memcpy(&data_count, bufin, sizeof(size_t)); - i += sizeof(size_t); - - - pbufout.resize(data_count); - - const Decoder decoder(bufin, i); - -#ifndef __EXCEPTIONS - if (i==-1) - return -1; -#endif - - const uint8_t *in_ptr = - decoder.Decode(bufin+i, bufin+bufinlen, - pbufout.data(), pbufout.data()+data_count); - -#ifndef __EXCEPTIONS - if (!in_ptr) - return -1; -#endif - - return in_ptr-bufin; - } -}; - -#endif Index: ct/tools/pyscripts/pyfact/izstream.h =================================================================== --- /fact/tools/pyscripts/pyfact/izstream.h (revision 17730) +++ (revision ) @@ -1,170 +1,0 @@ -#ifndef MARS_izstream -#define MARS_izstream - -#include - -#include -#include - -#ifdef __CINT__ -typedef void *gzFile; -#else -#include -#endif - -class izstream : public std::streambuf, public std::istream -{ -private: - static const int fgBufferSize = 2048*1024*2; - - gzFile fFile; // file handle for compressed file - char *fBuffer; // data buffer - - int underflow() - { - if (gptr() && gptr()(gptr()); - - if (!is_open()) - return EOF; - - // gptr()-eback(): if more than four bytes are already flushed - const int iputback = gptr()-eback()>4 ? 4 : gptr()-eback(); - - // Copy the last four bytes flushed into the putback area - memcpy(fBuffer+(4-iputback), gptr()-iputback, iputback); - - // Fill the buffer starting at the current file position and reset buffer - // pointers by calling setg - const int num = gzread(fFile, fBuffer+4, fgBufferSize-4); - if (num <= 0) // ERROR or EOF - return EOF; - - // reset buffer pointers - setg(fBuffer+(4-iputback), fBuffer+4, fBuffer+4+num); - - // return next character - return *reinterpret_cast(gptr()); - } - - -public: - izstream() : std::istream(this), fFile(0) - { - fBuffer = new char[fgBufferSize]; - setg(fBuffer+4, fBuffer+4, fBuffer+4); - } - izstream(const char *name) : std::istream(this), fFile(0) - { - fBuffer = new char[fgBufferSize]; - setg(fBuffer+4, fBuffer+4, fBuffer+4); - open(name); - } - ~izstream() { izstream::close(); delete [] fBuffer; } - - int is_open() { return fFile!=0; } - - // -------------------------------------------------------------------------- - // - // Open a file by name. Test if it is open like for an ifstream - // It doesn't matter whether the file is gzip compressed or not. - // - void open(const char* name) - { - if (is_open()) - { - clear(rdstate()|std::ios::failbit); - return; - } - - fFile = gzopen(name, "rb"); - if (fFile == 0) - { - clear(rdstate()|std::ios::failbit); - return; - } - } - // -------------------------------------------------------------------------- - // - // Close an open file. - // - void close() - { - if (!is_open()) - return; - - if (gzclose(fFile) != Z_OK) - clear(rdstate()|std::ios::failbit); - - fFile = 0; - } - - std::streambuf::pos_type seekoff(std::streambuf::off_type offset, std::ios_base::seekdir dir, - std::ios_base::openmode = std::ios_base::in) - { - // Using a switch instead results in: - // In member function `virtual std::streampos izstream::seekoff(long int, std::_Ios_Seekdir, std::_Ios_Openmode)': - // warning: enumeration value `_M_ios_seekdir_end' not handled in switch - // warning: case value `0' not in enumerated type `_Ios_Seekdir' - // warning: case value `1' not in enumerated type `_Ios_Seekdir' - // warning: case value `2' not in enumerated type `_Ios_Seekdir' - - if (dir==std::ios::end) - { - clear(rdstate()|std::ios::failbit); - return EOF; - } - - // We only do relative seeking to avoid unnecessary decompression - // of the whole file - if (dir==std::ios::beg) - offset -= tellg(); - - // Calculate future position in streambuffer - const char *ptr = gptr()+offset; - - // This is the number of bytes still available in the buffer - const size_t sbuf = egptr()-gptr(); - - // Check if the new position will still be in the buffer - // In this case the target data was already decompressed. - if (ptr>=eback() && ptr("ZNAXIS2"); - else - return fTable.Get("NAXIS2"); - } - size_t GetBytesPerRow() const - { - if (fTable.is_compressed) - return fTable.Get("ZNAXIS1"); - else - return fTable.Get("NAXIS1"); - } - -protected: - - // Stage the requested row to internal buffer - // Does NOT return data to users - virtual void StageRow(size_t row, char* dest) - { - if (!fTable.is_compressed) - { - fits::StageRow(row, dest); - return; - } - ReadBinaryRow(row, dest); - } - -private: - - // Do what it takes to initialize the compressed structured - void InitCompressionReading() - { - fCatalogInitialized = true; - - if (!fTable.is_compressed) - return; - - //The constructor may have failed - if (!good()) - return; - - if (fTable.is_compressed) - for (auto it=fTable.sorted_cols.cbegin(); it!= fTable.sorted_cols.cend(); it++) - { - if (it->comp == kCompFACT) - continue; - - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("Only the FACT compression scheme is handled by this reader."); -#else - gLog << ___err___ << "ERROR - Only the FACT compression scheme is handled by this reader." << std::endl; - return; -#endif - } - - fColumnOrdering.resize(fTable.sorted_cols.size(), FITS::kOrderByRow); - - //Get compressed specific keywords - fNumTiles = fTable.is_compressed ? GetInt("NAXIS2") : 0; - fNumRowsPerTile = fTable.is_compressed ? GetInt("ZTILELEN") : 0; - - //read the file's catalog - ReadCatalog(); - - //give it some space for uncompressing - AllocateBuffers(); - } - - // Copy decompressed data to location requested by user - void MoveColumnDataToUserSpace(char* dest, const char* src, const Table::Column& c) - { - if (!fTable.is_compressed) - { - fits::MoveColumnDataToUserSpace(dest, src, c); - return; - } - - memcpy(dest, src, c.num*c.size); - } - - bool fCatalogInitialized; - - std::vector fBuffer; /// fTransposedBuffer; /// fCompressedBuffer; /// fColumnOrdering; ///< ordering of the column's rows. Can change from tile to tile. - - size_t fNumTiles; ///< Total number of tiles - size_t fNumRowsPerTile; ///< Number of rows per compressed tile - int64_t fCurrentRow; ///< current row in memory signed because we need -1 - size_t fShrinkFactor; ///< shrink factor - - streamoff fHeapOff; ///< offset from the beginning of the file of the binary data - streamoff fHeapFromDataStart; ///< offset from the beginning of the data table - - std::vector>> fCatalog; ///< Catalog, i.e. the main table that points to the compressed data. - std::vector fTileSize; ///< size in bytes of each compressed tile - std::vector> fTileOffsets; ///< offset from start of tile of a given compressed column - - Checksum fRawsum; ///< Checksum of the uncompressed, raw data - - // Get buffer space - void AllocateBuffers() - { - uint32_t buffer_size = fTable.bytes_per_row*fNumRowsPerTile; - uint32_t compressed_buffer_size = fTable.bytes_per_row*fNumRowsPerTile + - //use a bit more memory for block headers. 256 char coding the compression sequence max. - fTable.num_cols*(sizeof(FITS::BlockHeader)+256) + - //a bit more for the tile headers - sizeof(FITS::TileHeader) + - //and a bit more for checksuming - 8; - - if (buffer_size % 4 != 0) - buffer_size += 4 - (buffer_size%4); - - if (compressed_buffer_size % 4 != 0) - compressed_buffer_size += 4 - (compressed_buffer_size%4); - - fBuffer.resize(buffer_size); - - fTransposedBuffer.resize(buffer_size); - fCompressedBuffer.resize(compressed_buffer_size); - } - - // Read catalog data. I.e. the address of the compressed data inside the heap - void ReadCatalog() - { - std::vector readBuf(16); - fCatalog.resize(fNumTiles); - - const streampos catalogStart = tellg(); - - fChkData.reset(); - - //do the actual reading - for (uint32_t i=0;i(reinterpret_cast(tempValues), readBuf.data(), 2); - if (tempValues[0] < 0 || tempValues[1] < 0) - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("Negative value in the catalog"); -#else - gLog << ___err___ << "ERROR - negative value in the catalog" << std::endl; - return; -#endif - } - //add catalog entry - fCatalog[i].emplace_back(tempValues[0], tempValues[1]); - } - - //see if there is a gap before heap data - fHeapOff = tellg()+fTable.GetHeapShift(); - fHeapFromDataStart = fNumTiles*fTable.num_cols*2*sizeof(int64_t) + fTable.GetHeapShift(); - - //check if the catalog has been shrinked - fShrinkFactor = HasKey("ZSHRINK") ? GetUInt("ZSHRINK") : 1; - - if (fNumRowsPerTile%fShrinkFactor) - { - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("Rows per tile and shrink factor do not match"); -#else - gLog << ___err___ << "ERROR - Rows per tile and shrink factor do not match" << std::endl; - return; -#endif - } - - if (fShrinkFactor>0) - fNumRowsPerTile /= fShrinkFactor; - - //compute the total size of each compressed tile - fTileSize.resize(fNumTiles); - fTileOffsets.resize(fNumTiles); - for (uint32_t i=0;i buf(catSize); - read(buf.data(), catSize); - - fCopy.write(buf.data(), catSize); - if (!fCopy) - clear(rdstate()|std::ios::badbit); - } - - //overrides fits.h method with empty one - //work is done in ReadBinaryRow because it requires volatile data from ReadBinaryRow - virtual void WriteRowToCopyFile(size_t row) - { - if (row == fRow+1) - fRawsum.add(fBufferRow); - } - - // Compressed version of the read row, even files with shrunk catalogs - // can be read fully sequentially so that streaming, e.g. through - // stdout/stdin, is possible. - bool ReadBinaryRow(const size_t &rowNum, char *bufferToRead) - { - if (rowNum >= GetNumRows()) - return false; - - if (!fCatalogInitialized) - InitCompressionReading(); - - // Book keeping, where are we? - const int64_t requestedTile = rowNum / fNumRowsPerTile; - const int64_t currentTile = fCurrentRow / fNumRowsPerTile; - - const int64_t requestedSuperTile = requestedTile / fShrinkFactor; - //const int64_t currentSuperTile = currentTile / fShrinkFactor; - - const int64_t requestedSubTile = requestedTile % fShrinkFactor; - //const int64_t currentSubTile = currentTile % fShrinkFactor; - - // Is this the first tile we read at all? - const bool isFirstTile = fCurrentRow<0; - - // Is this just the next tile in the sequence? - const bool isNextTile = requestedTile==currentTile+1 || isFirstTile; - - fCurrentRow = rowNum; - - // Do we have to read a new tile from disk? - if (requestedTile!=currentTile || isFirstTile) - { - //skip to the beginning of the tile - const int64_t superTileStart = fCatalog[requestedSuperTile][0].second - sizeof(FITS::TileHeader); - - std::vector offsets = fTileOffsets[requestedSuperTile]; - - // If this is a sub tile we might have to step forward a bit and - // seek for the sub tile. If we were just reading the previous one - // we can skip that. - if (!isNextTile || isFirstTile) - { - // step to the beginnig of the super tile - seekg(fHeapOff+superTileStart); - - // If there are sub tiles we might have to seek through the super tile - for (uint32_t k=0; k0) - { - // Read header - read(destBuffer, sizeof(FITS::TileHeader)); - - // Get size of tile - currentTileSize = reinterpret_cast(destBuffer)->size; - - // now read the remaining bytes of this tile - read(destBuffer+sizeof(FITS::TileHeader), currentTileSize-sizeof(FITS::TileHeader)); - - // Calculate the offsets recursively - offsets[0] = 0; - - //skip through the columns - for (size_t i=0; i(pos)->size; - } - } - else - { - // If we are reading the first tile of a super tile, all information - // is already available. - currentTileSize = fTileSize[requestedSuperTile] + sizeof(FITS::TileHeader); - read(destBuffer, currentTileSize); - } - - - // If we are reading sequentially, calcualte checksum - if (isNextTile) - { - // Padding for checksum calculation - memset(fCompressedBuffer.data(), 0, offset); - memset(destBuffer+currentTileSize, 0, fCompressedBuffer.size()-currentTileSize-offset); - fChkData.add(fCompressedBuffer); - } - - // Check if we are writing a copy of the file - if (isNextTile && fCopy.is_open() && fCopy.good()) - { - fCopy.write(fCompressedBuffer.data()+offset, currentTileSize); - if (!fCopy) - clear(rdstate()|std::ios::badbit); - } - else - if (fCopy.is_open()) - clear(rdstate()|std::ios::badbit); - - - // uncompress the buffer - const uint32_t thisRoundNumRows = (GetNumRows()offset; // pointer to column (destination buffer) - - switch (fColumnOrdering[i]) - { - case FITS::kOrderByRow: - // regular, "semi-transposed" copy - for (char *dest=buffer; destbytes); - src += it->bytes; // next column - } - break; - - case FITS::kOrderByCol: - // transposed copy - for (char *elem=buffer; elembytes; elem+=it->size) // element-by-element (arrays) - { - for (char *dest=elem; destsize); - src += it->size; // next element - } - } - break; - - default: - clear(rdstate()|std::ios::badbit); -#ifdef __EXCEPTIONS - throw std::runtime_error("Unkown column ordering scheme found"); -#else - gLog << ___err___ << "ERROR - unkown column ordering scheme" << std::endl; - return false; -#endif - break; - }; - } - } - - //Data loaded and uncompressed. Copy it to destination - memcpy(bufferToRead, fBuffer.data()+fTable.bytes_per_row*(fCurrentRow%fNumRowsPerTile), fTable.bytes_per_row); - return good(); - } - - // Read a bunch of uncompressed data - uint32_t UncompressUNCOMPRESSED(char* dest, - const char* src, - uint32_t numElems, - uint32_t sizeOfElems) - { - memcpy(dest, src, numElems*sizeOfElems); - return numElems*sizeOfElems; - } - - // Read a bunch of data compressed with the Huffman algorithm - uint32_t UncompressHUFFMAN16(char* dest, - const char* src, - uint32_t numChunks) - { - std::vector uncompressed; - - //read compressed sizes (one per row) - const uint32_t* compressedSizes = reinterpret_cast(src); - src += sizeof(uint32_t)*numChunks; - - //uncompress the rows, one by one - uint32_t sizeWritten = 0; - for (uint32_t j=0;j(src), compressedSizes[j], uncompressed); - - memcpy(dest, uncompressed.data(), uncompressed.size()*sizeof(uint16_t)); - - sizeWritten += uncompressed.size()*sizeof(uint16_t); - dest += uncompressed.size()*sizeof(uint16_t); - src += compressedSizes[j]; - } - return sizeWritten; - } - - // Apply the inverse transform of the integer smoothing - uint32_t UnApplySMOOTHING(int16_t* data, - uint32_t numElems) - { - //un-do the integer smoothing - for (uint32_t j=2;j &offsets, - const uint32_t &thisRoundNumRows, - const uint32_t offset) - { - char *dest = fTransposedBuffer.data(); - - //uncompress column by column - for (uint32_t i=0; i(&fCompressedBuffer[compressedOffset]); - - fColumnOrdering[i] = head->ordering; - - const uint32_t numRows = (head->ordering==FITS::kOrderByRow) ? thisRoundNumRows : col.num; - const uint32_t numCols = (head->ordering==FITS::kOrderByCol) ? thisRoundNumRows : col.num; - - const char *src = fCompressedBuffer.data()+compressedOffset+sizeof(FITS::BlockHeader)+sizeof(uint16_t)*head->numProcs; - - for (int32_t j=head->numProcs-1;j >= 0; j--) - { - uint32_t sizeWritten=0; - - switch (head->processings[j]) - { - case FITS::kFactRaw: - sizeWritten = UncompressUNCOMPRESSED(dest, src, numRows*numCols, col.size); - break; - - case FITS::kFactSmoothing: - sizeWritten = UnApplySMOOTHING(reinterpret_cast(dest), numRows*numCols); - break; - - case FITS::kFactHuffman16: - sizeWritten = UncompressHUFFMAN16(dest, src, numRows); - break; - - default: - clear(rdstate()|std::ios::badbit); - - std::ostringstream str; - str << "Unkown processing applied to data. Col " << i << " proc " << j << " out of " << (int)head->numProcs; -#ifdef __EXCEPTIONS - throw std::runtime_error(str.str()); -#else - gLog << ___err___ << "ERROR - Unknown processing applied to data. Aborting" << std::endl; - return; -#endif - } - //increment destination counter only when processing done. - if (j==0) - dest+= sizeWritten; - } - } - } - - void CheckIfFileIsConsistent(bool update_catalog=false) - { - //goto start of heap - const streamoff whereAreWe = tellg(); - seekg(fHeapOff); - - //init number of rows to zero - uint64_t numRows = 0; - - //get number of columns from header - const size_t numCols = fTable.num_cols; - - std::vector > > catalog; - - FITS::TileHeader tileHead; - FITS::BlockHeader columnHead; - - streamoff offsetInHeap = 0; - //skip through the heap - while (true) - { - read((char*)(&tileHead), sizeof(FITS::TileHeader)); - //end of file - if (!good()) - break; - - //padding or corrupt data - if (memcmp(tileHead.id, "TILE", 4)) - { - clear(rdstate()|std::ios::badbit); - break; - } - - //a new tile begins here - catalog.emplace_back(); - offsetInHeap += sizeof(FITS::TileHeader); - - //skip through the columns - for (size_t i=0;i