source: trunk/Mars/mcore/zfits.h@ 16815

/*
 * zfits.h
 *
 *  Created on: May 16, 2013
 *      Author: lyard
 */

#ifndef MARS_ZFITS
#define MARS_ZFITS

#include <stdexcept>

#include "fits.h"
#include "huffman.h"


#define FACT_RAW       0x0
#define FACT_SMOOTHING 0x1
#define FACT_HUFFMAN16 0x2

#define FACT_COL_MAJOR 'C'
#define FACT_ROW_MAJOR 'R'


#ifndef __MARS__
namespace std
{
#endif

class zfits : public fits
{
public:

    // Basic constructor
    zfits(const string& fname, const string& tableName="",
          bool force=false) : fits(fname, tableName, force),
        fBuffer(0),
        fTransposedBuffer(0),
        fCompressedBuffer(0),
        fNumTiles(0),
        fNumRowsPerTile(0),
        fCurrentRow(-1),
        fHeapOff(0),
        fTileSize(0)
    {
        InitCompressionReading();
    }

    // Alternative contstructor
    zfits(const string& fname, const string& fout, const string& tableName,
          bool force=false) : fits(fname, fout, tableName, force),
              fBuffer(0),
              fTransposedBuffer(0),
              fCompressedBuffer(0),
              fNumTiles(0),
              fNumRowsPerTile(0),
              fCurrentRow(-1),
              fHeapOff(0),
              fTileSize(0)
    {
        InitCompressionReading();
    }

    //  Skip the next row
    bool SkipNextRow()
    {
        if (!fTable.isCompressed)
            return fits::SkipNextRow();

        fRow++;
        return true;
    }
protected:

    //  Stage the requested row to internal buffer
    //  Does NOT return data to users
    virtual void StageRow(size_t row, char* dest)
    {
        if (!fTable.isCompressed)
        {
            fits::StageRow(row, dest);
            return;
        }

        ReadBinaryRow(row, dest);
    }

private:

    //Structure helper for reading tiles headers
    typedef 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__)) TileHeader;

    //Structure helper for reading blocks headers and compresion schemes
    typedef struct BlockHeader
    {
        uint64_t      size;
        char          ordering;
        unsigned char numProcs;
        uint16_t      processings[];

        BlockHeader(uint64_t      s=0,
                    char          o=FACT_ROW_MAJOR,
                    unsigned char n=1) : size(s),
                                         ordering(o),
                                         numProcs(n)
        {}
    } __attribute__((__packed__)) BlockHeader;

    // Do what it takes to initialize the compressed structured
    void InitCompressionReading()
    {
        //The constructor may have failed
        if (!good())
            return;

        if (fTable.isCompressed)
        for (auto it=fTable.sortedCols.begin(); it!= fTable.sortedCols.end(); it++)
            if (it->comp != FACT)
            {
#ifdef __EXCEPTIONS
                    throw runtime_error("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." << endl;
                    return;
#endif
            }

        fColumnOrdering.resize(fTable.sortedCols.size());
        for (auto it=fColumnOrdering.begin(); it != fColumnOrdering.end(); it++)
            (*it) = FACT_ROW_MAJOR;
        //Get compressed specific keywords
        fNumTiles       = fTable.isCompressed ? GetInt("NAXIS2") : 0;
        fNumRowsPerTile = fTable.isCompressed ? GetInt("ZTILELEN") : 0;

        //give it some space for uncompressing
        AllocateBuffers();

        //read the file's catalog
        ReadCatalog();
    }

    // Copy decompressed data to location requested by user
    void MoveColumnDataToUserSpace(char* dest, const char* src, const Table::Column& c)
    {
        if (!fTable.isCompressed)
        {
            fits::MoveColumnDataToUserSpace(dest, src, c);
            return;
        }

        memcpy(dest, src, c.num*c.size);
    }

    vector<char> fBuffer;           ///<store the uncompressed rows
    vector<char> fTransposedBuffer; ///<intermediate buffer to transpose the rows
    vector<char> fCompressedBuffer; ///<compressed rows
    vector<char> fColumnOrdering; ///< ordering of the column's rows

    size_t fNumTiles;       ///< Total number of tiles
    size_t fNumRowsPerTile; ///< Number of rows per compressed tile
    size_t fCurrentRow;     ///< current row in memory.

    streamoff fHeapOff;     ///< offset from the beginning of the file of the binary data

    vector<vector<pair<int64_t, int64_t> > > fCatalog;///< Catalog, i.e. the main table that points to the compressed data.
    vector<size_t>                           fTileSize; ///< size in bytes of each compressed tile
    vector<vector<size_t> >                  fTileOffsets; ///< offset from start of tile of a given compressed column

    void AllocateBuffers()
    {
        if (!fTable.isCompressed)
            return;

        fBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile);

        fTransposedBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile);
        fCompressedBuffer.resize(fTable.bytes_per_row*fNumRowsPerTile + fTable.num_cols*(sizeof(BlockHeader)+256)); //use a bit more memory for block headers
    }

    // Read catalog data. I.e. the address of the compressed data inside the heap
    void ReadCatalog()
    {
        if (!fTable.isCompressed)
            return;

        char readBuf[16];
        fCatalog.resize(fNumTiles);

        const streampos catalogStart = tellg();

        //do the actual reading
        for (uint32_t i=0;i<fNumTiles;i++)
            for (uint32_t j=0;j<fTable.num_cols;j++)
            {
                read(readBuf, 2*sizeof(int64_t));

                //swap the bytes
                int64_t tempValues[2] = {0,0};
                revcpy<8>(reinterpret_cast<char*>(tempValues), readBuf, 2);
                if (tempValues[0] < 0 || tempValues[1] < 0)
                {
                    clear(rdstate()|ios::badbit);
#ifdef __EXCEPTIONS
                    throw runtime_error("Negative value in the catalog");
#else
                    gLog << ___err___ << "ERROR - negative value in the catalog" << endl;
                    return;
#endif
                }
                //add catalog entry
                fCatalog[i].emplace_back(tempValues[0], tempValues[1]);
            }

        //compute the total size of each compressed tile
        fTileSize.resize(fNumTiles);
        fTileOffsets.resize(fNumTiles);
        for (uint32_t i=0;i<fNumTiles;i++)
        {
            fTileSize[i] = 0;
            for (uint32_t j=0;j<fTable.num_cols;j++)
            {
                fTileSize[i] += fCatalog[i][j].first;
                fTileOffsets[i].emplace_back(fCatalog[i][j].second - fCatalog[i][0].second);
            }
        }
        //see if there is a gap before heap data
        fHeapOff = tellg()+fTable.GetHeapShift();

        if (!fCopy.is_open())
            return;

        //write catalog and heap gap to target file
        seekg(catalogStart);

        const size_t catSize = fTable.GetHeapShift() + fTable.total_bytes;

        vector<char> buf(catSize);
        read(buf.data(), catSize);

        fCopy.write(buf.data(), catSize);
        if (!fCopy)
            clear(rdstate()|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 )
    {

    }
    // Compressed versin of the read row
    bool ReadBinaryRow(const size_t &rowNum, char *bufferToRead)
    {
        if (rowNum >= GetNumRows())
            return false;

        const uint32_t requestedTile = rowNum/fNumRowsPerTile;
        const uint32_t currentTile   = fCurrentRow/fNumRowsPerTile;
        const size_t   previousRow   = fCurrentRow;

        fCurrentRow = rowNum;

        //should we read yet another chunk of data ?
        if (requestedTile != currentTile)
        {
            //read yet another chunk from the file
            const int64_t sizeToRead = fTileSize[requestedTile];

            //skip to the beginning of the tile
            seekg(fHeapOff+fCatalog[requestedTile][0].second - sizeof(TileHeader));
            TileHeader tHead;
            read((char*)(&tHead), sizeof(TileHeader));
            read(fCompressedBuffer.data(), sizeToRead);

            if (fCurrentRow == previousRow+1 &&
                fCopy.is_open() &&
                fCopy.good())
            {
                fCopy.write((char*)(&tHead), sizeof(TileHeader));
                fCopy.write(fCompressedBuffer.data(), sizeToRead);
                if (!fCopy)
                    clear(rdstate()|ios::badbit);
            }
            else
                if (fCopy.is_open())
                    clear(rdstate()|ios::badbit);

            const uint32_t thisRoundNumRows = (GetNumRows()<fCurrentRow + fNumRowsPerTile) ? GetNumRows()%fNumRowsPerTile : fNumRowsPerTile;

            //uncompress it
            UncompressBuffer(requestedTile, thisRoundNumRows);

            // pointer to column (source buffer)
            const char *src = fTransposedBuffer.data();

            uint32_t i=0;
            for (auto it=fTable.sortedCols.begin(); it!=fTable.sortedCols.end(); it++, i++)
            {
                char *buffer = fBuffer.data() + it->offset; // pointer to column (destination buffer)

                switch (fColumnOrdering[i])
                {
                    case FACT_ROW_MAJOR:
                        // regular, "semi-transposed" copy
                        for (char *dest=buffer; dest<buffer+thisRoundNumRows*fTable.bytes_per_row; dest+=fTable.bytes_per_row) // row-by-row
                        {
                            memcpy(dest, src, it->bytes);
                            src += it->bytes;  // next column
                        }
                    break;

                    case FACT_COL_MAJOR:
                        // transposed copy
                        for (char *elem=buffer; elem<buffer+it->bytes; elem+=it->size) // element-by-element (arrays)
                        {
                            for (char *dest=elem; dest<elem+thisRoundNumRows*fTable.bytes_per_row; dest+=fTable.bytes_per_row) // row-by-row
                            {
                                memcpy(dest, src, it->size);
                                src += it->size; // next element
                            }
                        }
                    break;
                    default:
    #ifdef __EXCEPTIONS
                        throw runtime_error("Unkown column ordering scheme");
    #else
                        gLog << ___err___ << "ERROR - unkown column ordering scheme" << endl;
                        return;
    #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)
    {
        vector<uint16_t> uncompressed;

        //read compressed sizes (one per row)
        const uint32_t* compressedSizes = reinterpret_cast<const uint32_t*>(src);
        src += sizeof(uint32_t)*numChunks;

        //uncompress the rows, one by one
        uint32_t sizeWritten = 0;
        for (uint32_t j=0;j<numChunks;j++)
        {
            Huffman::Decode(reinterpret_cast<const unsigned char*>(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;
    }

    uint32_t UnApplySMOOTHING(int16_t*   data,
                              uint32_t numElems)
    {
        //un-do the integer smoothing
        for (uint32_t j=2;j<numElems;j++)
            data[j] = data[j] + (data[j-1]+data[j-2])/2;

        return numElems*sizeof(uint16_t);
    }
    // Data has been read from disk. Uncompress it !
    void UncompressBuffer(const uint32_t &catalogCurrentRow, const uint32_t &thisRoundNumRows)
    {
        char *dest = fTransposedBuffer.data();

        //uncompress column by column
        for (uint32_t i=0; i<fTable.sortedCols.size(); i++)
        {
            const fits::Table::Column &col = fTable.sortedCols[i];
            if (col.num == 0)
                continue;

            //get the compression flag
            const int64_t compressedOffset = fTileOffsets[catalogCurrentRow][i];

            BlockHeader* head = reinterpret_cast<BlockHeader*>(&fCompressedBuffer[compressedOffset]);

            fColumnOrdering[i] = head->ordering;

            uint32_t numRows = (head->ordering==FACT_ROW_MAJOR) ? thisRoundNumRows : col.num;
            uint32_t numCols = (head->ordering==FACT_COL_MAJOR) ? thisRoundNumRows : col.num;

            const char *src = fCompressedBuffer.data()+compressedOffset+sizeof(BlockHeader)+sizeof(uint16_t)*head->numProcs;

            for (uint32_t j=head->numProcs;j != 0; j--)
            {
                uint32_t sizeWritten=0;

                switch (head->processings[j-1])
                {
                    case FACT_RAW:
                            if (head->numProcs == 1)
                                sizeWritten = UncompressUNCOMPRESSED(dest, src, numRows*numCols, col.size);
                    break;
                    case FACT_SMOOTHING:
                            sizeWritten = UnApplySMOOTHING(reinterpret_cast<int16_t*>(dest), numRows*numCols);
                    break;
                    case FACT_HUFFMAN16:
                            sizeWritten = UncompressHUFFMAN16(dest, src, numRows);
                    break;
                    default:
#ifdef __EXCEPTIONS
                    throw runtime_error("Unknown processing applied to data. Aborting");
#else
                    gLog << ___err___ << "ERROR - Unknown processing applied to data. Aborting" << endl;
                    return;
#endif
                    break;
                }
                //increment destination counter only when processing done.
                if (j==1) dest+= sizeWritten;
            }
        }
    }

};//class zfits

#ifndef __MARS__
}; //namespace std
#endif

#endif