GHdf5.cpp

Go to the documentation of this file.

/***************************************************************************
 *                   GHdf5.cpp - HDF5 file handling class                  *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2026 by Juergen Knoedlseder                              *
 * ----------------------------------------------------------------------- *
 *                                                                         *
 *  This program is free software: you can redistribute it and/or modify   *
 *  it under the terms of the GNU General Public License as published by   *
 *  the Free Software Foundation, either version 3 of the License, or      *
 *  (at your option) any later version.                                    *
 *                                                                         *
 *  This program is distributed in the hope that it will be useful,        *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 *  GNU General Public License for more details.                           *
 *                                                                         *
 *  You should have received a copy of the GNU General Public License      *
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
 *                                                                         *
 ***************************************************************************/
/**
 * @file GHdf5.hpp
 * @brief HDF5 file handling class implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <string>
#include <cstring>          // for std::memcpy, std::memmove, std::strncmp
#include <typeinfo>
#include "GHdf5.hpp"
#include "GException.hpp"
#include "GTools.hpp"
 
/* __ Method name definitions ____________________________________________ */
#define G_LOAD                                      "GHdf5::load(GFilename&)"
#define G_XML_HDF5_ENTRY                "GHdf5::xml_hdf5_entry(std::string&)"
#define G_READ_SUPERBLOCK                     "GHdf5::read_superblock(FILE*)"
#define G_READ_SYMBOL_TABLE_NODE      "GHdf5::read_symbol_table_node(FILE*, "\
                    "uint64_t&, uint64_t&, std::string&, GXmlElement*, int&)"
#define G_READ_SYMBOL_TABLE_ENTRY    "GHdf5::read_symbol_table_entry(FILE*, "\
                                             "uint64_t&, GXmlElement*, int&)"
#define G_READ_GROUP        "GHdf5::read_group(FILE*, uint64_t&, uint64_t&, "\
                                                        "GXmlElement*, int&)"
#define G_READ_BTREE_CHUNKED                      "GHdf5::read_btree(FILE*, "\
                                      "uint64_t&, GXmlElement*, int& indent)"
#define G_READ_OBJECT_HEADER              "GHdf5::read_object_header(FILE*, "\
                                            "uint64_t&,  GXmlElement*, int&)"
#define G_READ_OBJECT_HEADER_V1   "GHdf5::read_object_header_v1(FILE*, int&)"
#define G_READ_OBJECT_HEADER_V2   "GHdf5::read_object_header_v2(FILE*, int&)"
#define G_READ_OBJECT_HEADER_MESSAGE     "GHdf5::read_object_header_message("\
                                                         "FILE*, int&, int&)"
#define G_READ_MESSAGE_DATASPACE      "GHdf5::read_message_dataspace(FILE*, "\
                                                "GXmlElement*, std::string&)"
#define G_READ_MESSAGE_DATATYPE        "GHdf5::read_message_datatype(FILE*, "\
                                                "GXmlElement*, std::string&)"
#define G_READ_MESSAGE_LAYOUT            "GHdf5::read_message_layout(FILE*, "\
                                                "GXmlElement*, std::string&)"
#define G_READ_MESSAGE_FILTER            "GHdf5::read_message_filter(FILE*, "\
                                                "GXmlElement*, std::string&)"
#define G_READ_MESSAGE_ATTRIBUTE      "GHdf5::read_message_attribute(FILE*, "\
                                                "GXmlElement*, std::string&)"
#define G_GLOBAL_HEAP_STRING   "GHdf5::global_heap_string(FILE*, uint64_t&, "\
                                                           "int&, uint32_t&)"
#define G_FREAD_DATA1                        "GHdf5::fread_data(FILE*, int&)"
#define G_FREAD_DATA2      "gammalib::hdf5::fread_data(FILE*, GXmlElement*, "\
                                  "GXmlElement*, GXmlElement*, GXmlElement*)"
#define G_FREAD_DATA_CHUNK         "gammalib::hdf5::fread_data_chunk(FILE*, "\
      "GXmlElement*, GXmlElement*, GXmlElement*, GXmlElement*, GXmlElement*)"
#define G_FREAD_INT                  "gammalib::hdf5::fread_int(FILE*, int&)"
#define G_FREAD_UINT32            "gammalib::hdf5::fread_uint32(FILE*, int&)"
#define G_FREAD_UINT64            "gammalib::hdf5::fread_uint64(FILE*, int&)"
#define G_FREAD_STRING            "gammalib::hdf5::fread_string(FILE*, int&)"
#define G_FREAD_DATA_AS_STRING "gammalib::hdf5::fread_data_as_string(FILE*, "\
                                                              "GXmlElement*)"
#define G_FREAD_ZERO                "gammalib::hdf5::fread_zero(FILE*, int&)"
#define G_DATA_TO_INT      "gammalib::hdf5::data_to_int(char*, GXmlElement*)"
#define G_DATA_TO_DOUBLE             "gammalib::hdf5::data_to_double(char*, "\
                                                              "GXmlElement*)"
#define G_DATA_TO_STRING             "gammalib::hdf5::data_to_string(char*, "\
                                                              "GXmlElement*)"
#define G_DATA_FILTER            "gammalib::hdf5::data_filter(std::string&, "\
                                                              "GXmlElement*)"
#define G_DATA_FILTER_BITSHUFFLE \
         "gammalib::hdf5::data_filter_bitshuffle(std::string&, GXmlElement*)"
#define G_XML_MSG_TYPE          "gammalib::hdf5::xml_msg_type(GXmlElement*, "\
                                                   "std::string& type, int&)"
#define G_XML_MSG_ATTRIBUTE \
              "gammalib::hdf5::xml_msg_attribute(GXmlElement*, std::string&)"
#define G_DECOMPRESS_LZ4     "decompress_lz4(char*, char*, size_t&, size_t&)"
#define G_READ_VARIABLE_LENGTH   "read_variable_length(uint8_t**, uint8_t*, "\
                                                                     "bool&)"
#define G_MEMCPY_USING_OFFSET      "memcpy_using_offset(uint8_t*, uint8_t*, "\
                                                         "uint8_t*, size_t&)"
 
 
/* __ Macros _____________________________________________________________ */
 
/* __ Coding definitions _________________________________________________ */
 
/* __ Debug definitions __________________________________________________ */
//#define G_HDF5_DEBUG                          //!< Debug HDF5 file decoding
//#define G_BITSHUFFLE_DEBUG                    //!< Debug Bitshuffe filter
 
/* __ Constants __________________________________________________________ */
 
/* __ Prototypes for bitshuffle filter ___________________________________ */
void          bitshuffle_decompress(char*         input,
                                    char*         output,
                                    const size_t& elements,
                                    const size_t& element_size,
                                    size_t&       block_size);
uint64_t      bitshuffle_read_uint64(const char* buffer);
uint32_t      bitshuffle_read_uint32(const char* buffer);
void          bitshuffle_elements(const char*   input,
                                  char*         output,
                                  const size_t& elements,
                                  const size_t& element_size);
int           decompress_lz4(const char*   src,
                             char*         dst,
                             const size_t& srcSize,
                             const size_t& outputSize);
size_t        read_variable_length(const uint8_t** ip,
                                   const uint8_t*  ilimit,
                                   const bool&     initial_check);
uint16_t      read_uint16(const void* memPtr);
void          memcpy_using_offset(const uint8_t* srcPtr,
                                  uint8_t*       dstPtr,
                                  uint8_t*       dstEnd,
                                  const size_t&  offset);
void          memcpy_beyond8(const uint8_t* srcPtr,
                             uint8_t*       dstPtr,
                             uint8_t*       dstEnd);
void          memory_beyond32(const uint8_t* srcPtr,
                              uint8_t*       dstPtr,
                              uint8_t*       dstEnd);
 
/* __ Using of namespaces ________________________________________________ */
using namespace gammalib::hdf5;
 
 
/*==========================================================================
 =                                                                         =
 =                       Constructors/destructors                          =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Void constructor
 *
 * Creates an empty instance of the HDF5 class.
 ***************************************************************************/
GHdf5::GHdf5(void)
{
    // Initialise members
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Load constructor
 *
 * @param[in] filename File name of HDF5 file.
 *
 * Create HDF5 instance by loading data from a HDF5 file.
 ***************************************************************************/
GHdf5::GHdf5(const GFilename& filename)
{
    // Initialise members
    init_members();
 
    // Load file
    load(filename);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy constructor
 *
 * @param[in] file HDF5 file instance.
 **************************************************************************/
GHdf5::GHdf5(const GHdf5& file)
{
    // Initialise members
    init_members();
 
    // Copy members
    copy_members(file);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Destructor
 ***************************************************************************/
GHdf5::~GHdf5(void)
{
    // Free members
    free_members();
 
    // Return
    return;
}
 
 
/*==========================================================================
 =                                                                         =
 =                                Operators                                =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Assignment operator
 *
 * @param[in] file HDF5 file instance.
 * @return HDF5 file instance.
 *
 * Assign HDF5 file instance to this object. The assignment performs a deep
 * copy of all information, hence the original object from which the
 * assignment was made can be destroyed after this operation without any loss
 * of information.
 ***************************************************************************/
GHdf5& GHdf5::operator=(const GHdf5& file)
{
    // Execute only if object is not identical
    if (this != &file) {
 
        // Free members
        free_members();
 
        // Initialise members
        init_members();
 
        // Copy members
        copy_members(file);
 
    } // endif: object was not identical
 
    // Return this object
    return *this;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Public methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Clear instance
 *
 * Clears HDF5 file instance by resetting all class members to an initial
 * state. Any information that was present before will be lost.
 ***************************************************************************/
void GHdf5::clear(void)
{
    // Free class members
    free_members();
 
    // Initialise members
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Clone instance
 *
 * @return Pointer to deep copy of HDF5 file instance.
 ***************************************************************************/
GHdf5* GHdf5::clone(void) const
{
    return new GHdf5(*this);
}
 
 
/***********************************************************************//**
 * @brief Load data from HDF5 file into instance
 *
 * @param[in] filename File name of HDF5 file.
 *
 * @exception GException::file_error
 *            Unable to open specified file.
 *
 * Load data from HDF5 file into a HDF5 file.
 ***************************************************************************/
void GHdf5::load(const GFilename& filename)
{
    // Clear instance
    clear();
 
    // Store filename
    m_filename = filename;
 
    // Expand environment variables in filename
    std::string fname = gammalib::expand_env(filename.url());
 
    // Open HDF5 file (read-only)
    FILE* fptr = std::fopen(fname.c_str(), "rb");
    if (fptr == NULL) {
        std::string msg = "Unable to open file \""+fname+"\" for read access. "
                          "Please specify a readable file.";
        throw GException::file_error(G_LOAD, msg);
    }
 
    // Read file
    read(fptr);
 
    // Close file
    std::fclose(fptr);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read data from HDF5 file into instance
 *
 * @param[in] fptr File pointer to HDF5 file.
 *
 * Read data from HDF5 file into a HDF5 file.
 ***************************************************************************/
void GHdf5::read(FILE* fptr)
{
    // Clear instance
    clear();
 
    // Move file pointer to start of file
    std::fseek(fptr, 0, SEEK_SET);
 
    // Parse HDF5 file by reading its superblock
    read_superblock(fptr);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Return XML element to entry tag
 *
 * @param[in] name Entry name.
 * @return Pointer to XML element of given name.
 *
 * @exception GException::invalid_value
 *            No HDF5 metadata were loaded.
 *            No entry element with specified @p name found.
 *
 * Returns a XML pointer to an entry element with a given @p name.
 *
 * The method works recursevely down the entire XML structure of the HDF5
 * file and stops once the first entry element with a "name" attribute equal
 * to @p name was found. If no such element could be found, the method throws
 * an exception.
 ***************************************************************************/
const GXmlElement* GHdf5::xml_hdf5_entry(const std::string& name) const
{
    // Throw an exception if XML instance is empty
    if (m_xml.is_empty()) {
        std::string msg = "No HDF5 metadata were loaded. Please first "
                          "load metadata from a HDF5 file.";
        throw GException::invalid_value(G_XML_HDF5_ENTRY, msg);
    }
 
    // Get ROOT element
    const GXmlElement* root = m_xml.element(0);
 
    // Search element
    const GXmlElement* element = gammalib::xml_get_element_by_attribute(root,
                                                                        "entry",
                                                                        "name",
                                                                        name);
 
    // Throw an exception if no entry was found
    if (element == NULL) {
        std::string msg = "No entry with attribute "+name+"=\""+name+
                          "\" found in HDF5 metadata.";
        throw GException::invalid_value(G_XML_HDF5_ENTRY, msg);
    }
 
    // Return element
    return element;
}
 
 
/***********************************************************************//**
 * @brief Print HDF5 file
 *
 * @param[in] chatter Chattiness.
 * @return String containing HDF5 information.
 ***************************************************************************/
std::string GHdf5::print(const GChatter& chatter) const
{
    // Initialise result string
    std::string result;
 
    // Continue only if chatter is not silent
    if (chatter != SILENT) {
 
        // Append header
        result.append("=== GHdf5 ===");
 
        // Append information
        result.append("\n"+gammalib::parformat("Filename")+m_filename);
        result.append("\n"+gammalib::parformat("Size of Offsets")+gammalib::str(m_offset));
        result.append("\n"+gammalib::parformat("Size of Lengths")+gammalib::str(m_length));
 
        // If file name is not empty then append HDF5 metadata
        if (m_filename.length() > 0) {
 
            //TODO: Implement generic printing
 
        } // endif: file name was not empty
 
    } // endif: chatter was not silent
 
    // Return result
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Private methods                             =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Initialise class members
 ***************************************************************************/
void GHdf5::init_members(void)
{
    // Initialise members
    m_filename.clear();
    m_offset = 0;
    m_length = 0;
    m_xml.clear();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy class members
 *
 * @param[in] file HDF5 file instance.
 ***************************************************************************/
void GHdf5::copy_members(const GHdf5& file)
{
    // Copy members
    m_filename = file.m_filename;
    m_offset   = file.m_offset;
    m_length   = file.m_length;
    m_xml      = file.m_xml;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Delete class members
 ***************************************************************************/
void GHdf5::free_members(void)
{
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read HDF5 file superblock
 *
 * @param[in] fptr File pointer.
 *
 * @exception GException::invalid_value
 *            File is not a HDF5 file
 *            Superblock version out of valid range
 *            Invalid offset and length values
 *
 * Reads the superblock of the HDF5 file and walks done the hierarchical HDF5
 * file structure to extract all metadata. The information is stored in the
 * XML file member m_xml that will echo the hierarchical information of
 * the HDF5 file.
 *
 * The superblock may begin at certain predefined offsets within the HDF5
 * file, allowing a block of unspecified content for users to place
 * additional information at the beginning (and end) of the HDF5 file
 * without limiting the HDF5 library's ability to manage the objects within
 * the file itself.
 *
 * This method assumes that the @p fptr file pointer points to the
 * superblock. If an offset is needed, the client needs to position the
 * file pointer accordingly.
 *
 * This method fully supports superblock versions 0 and 1. For superblock
 * versions 2 and 3, entering of the group graph still needs to be
 * implemented.
 ***************************************************************************/
void GHdf5::read_superblock(FILE* fptr)
{
    // Set HDF5 signature
    const char hdf5_signature[] = {static_cast<char>(0x89), 0x48, 0x44, 0x46, 0x0d, 0x0a, 0x1a, 0x0a};
 
    // Read signature
    std::string signature = fread_string(fptr, 8);
 
    // Check signature to verify that the file is indeed a HDF5 file
    if (std::strncmp(signature.c_str(), &hdf5_signature[0], 8) != 0) {
        std::string msg = "File does not contain a valid HDF5 signature. "
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
    }
 
    // Read version
    uint32_t version = fread_uint32(fptr, 1);
 
    // Check version
    if (version > 3) {
        std::string msg = "Superblock version "+gammalib::str(version)+" is "
                          "not comprised between 0 and 3. Please specify a "
                          "valid HDF5 file.";
        throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
    }
 
    // Parse version 0 or 1
    if (version < 2) {
 
        // Read version fields
        int version_free = fread_int(fptr, 1);
        int version_root = fread_int(fptr, 1);
        fread_zero(fptr, 1);
        int version_shared = fread_int(fptr, 1);
 
        // Read offset and length
        m_offset = fread_int(fptr, 1);
        m_length = fread_int(fptr, 1);
 
        // Check offset and length
        if (m_offset < 1) {
            std::string msg = "Offset "+gammalib::str(m_offset)+" is smaller "
                              "than 1 Byte. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
        }
        if (m_offset > 8) {
            std::string msg = "Offset "+gammalib::str(m_offset)+" is larger "
                              "than 8 Bytes. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
        }
        if (m_length < 1) {
            std::string msg = "Length "+gammalib::str(m_offset)+" is smaller "
                              "than 1 Byte. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
        }
        if (m_length > 8) {
            std::string msg = "Length "+gammalib::str(m_offset)+" is larger "
                              "than 8 Bytes. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
        }
 
        // Read further fields
        fread_zero(fptr, 1);                          // Zero padding
        int group_leaf_node     = fread_int(fptr, 2); // Group Leaf Node K
        int group_internal_node = fread_int(fptr, 2); // Group Internal Node K
        int flags               = fread_int(fptr, 4); // File Consistency Flags
        int indexed_storage     = 0;
 
        // Read optional field
        if (version == 1) {
            indexed_storage = fread_int(fptr, 2);
            fread_zero(fptr, 2);
        }
 
        // Read addresses
        uint64_t address_base = fread_uint64(fptr, m_offset); // Base Address
        uint64_t address_free = fread_uint64(fptr, m_offset); // Address of File Free Space Info
        uint64_t address_eof  = fread_uint64(fptr, m_offset); // End of File Address
        uint64_t address_drv  = fread_uint64(fptr, m_offset); // Driver Information Block Address
 
        // Debugging: print superblock information
        #if defined(G_HDF5_DEBUG)
        std::string msg = "Superblock(version="+gammalib::str(version)+
                          ", offset="+gammalib::str(m_offset)+
                          ", length="+gammalib::str(m_length)+")";
        debug_msg(msg, 0);
        #endif
 
        // Get current file pointer position as address of Root Group
        // Symbol Table Entry
        uint64_t pos = std::ftell(fptr);
 
        // Create symbol table entry XML element and append it to XML file
        GXmlElement* entry = m_xml.append("entry");
 
        // Add attributes to "ROOT" entry
        entry->attribute("version", gammalib::str(version));
        entry->attribute("offset", gammalib::str(m_offset));
        entry->attribute("length", gammalib::str(m_length));
        entry->attribute("version_free_space_storage", gammalib::str(version_free));
        entry->attribute("version_root_group_symbol_table", gammalib::str(version_root));
        entry->attribute("version_shared_header_message", gammalib::str(version_shared));
        entry->attribute("group_leaf_node", gammalib::str(group_leaf_node));
        entry->attribute("group_internal_node", gammalib::str(group_internal_node));
        entry->attribute("flags", gammalib::str(flags));
        if (version == 1) {
            entry->attribute("indexed_storage_internal_node", gammalib::str(indexed_storage));
        }
        entry->attribute("address_base", gammalib::str(address_base));
        entry->attribute("address_free_space_storage", gammalib::str(address_free));
        entry->attribute("address_eof", gammalib::str(address_eof));
        entry->attribute("address_driver_info", gammalib::str(address_drv));
 
        // Root Group Symbol Table Entry
        read_symbol_table_entry(fptr, pos, 0, "ROOT", entry, 1);
 
    } // endif: version 0 or 1
 
    // Parse version 2 or 3
    else {
 
        // Read offset and length
        m_offset = fread_int(fptr, 1);
        m_length = fread_int(fptr, 1);
 
        // Check offset and length
        if (m_offset < 1) {
            std::string msg = "Offset "+gammalib::str(m_offset)+" is smaller "
                              "than 1 Byte. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
        }
        if (m_offset > 8) {
            std::string msg = "Offset "+gammalib::str(m_offset)+" is larger "
                              "than 8 Bytes. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
        }
        if (m_length < 1) {
            std::string msg = "Length "+gammalib::str(m_offset)+" is smaller "
                              "than 1 Byte. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
        }
        if (m_length > 8) {
            std::string msg = "Length "+gammalib::str(m_offset)+" is larger "
                              "than 8 Bytes. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_SUPERBLOCK, msg);
        }
 
        // Read further fields
        int flags = fread_int(fptr, 1); // File Consistency Flags
 
        // Read addresses
        uint64_t address_base = fread_uint64(fptr, m_offset); // Base Address
        uint64_t address_ext  = fread_uint64(fptr, m_offset); // Superblock Extension Address
        uint64_t address_eof  = fread_uint64(fptr, m_offset); // End of File Address
        uint64_t address_root = fread_uint64(fptr, m_offset); // Root Group Object Header Address
 
        // Read further fields
        uint32_t chksum = fread_uint32(fptr, 4); // Superblock Checksum
 
        // Create symbol table entry XML element and append it to XML file
        GXmlElement* entry = m_xml.append("entry");
 
        // Add attributes to "ROOT" entry
        entry->attribute("version", gammalib::str(version));
        entry->attribute("offset", gammalib::str(m_offset));
        entry->attribute("length", gammalib::str(m_length));
        entry->attribute("flags", gammalib::str(flags));
        entry->attribute("address_base", gammalib::str(address_base));
        entry->attribute("address_superblock_extension", gammalib::str(address_ext));
        entry->attribute("address_eof", gammalib::str(address_eof));
        entry->attribute("address_root_group_object_header", gammalib::str(address_root));
        entry->attribute("checksum", gammalib::str(chksum));
 
        // Enter group graph
        //TODO
 
    } // endelse: version 2 or 3
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read HDF5 symbol table node
 *
 * @param[in] fptr File pointer.
 * @param[in] address Address of symbol table node.
 * @param[in] offset Byte offset in heap.
 * @param[in] heap Heap.
 * @param[in] table Pointer to XML element for table.
 * @param[in] indent Indent for debugging messages.
 *
 * @exception GException::invalid_value
 *            Symbol table node signature not found
 *            Invalid symbol table node version
 *
 * This method reads a symbol table node from the file. The symbol table
 * node is a leaf of a B-tree.
 *
 * Each symbol table node starts with the signature "SNOD" and contains an
 * arbitrary number of symbols (group entries). All symbol table nodes
 * have the same length of 2*m_offset+24.
 ***************************************************************************/
void GHdf5::read_symbol_table_node(FILE*              fptr,
                                   const uint64_t&    address,
                                   const uint64_t&    offset,
                                   const std::string& heap,
                                   GXmlElement*       table,
                                   const int&         indent)
{
    // Set symbol table node signature
    const char table_signature[] = {'S', 'N', 'O', 'D'};
 
    // Move file pointer to symbol table node address
    std::fseek(fptr, address, SEEK_SET);
 
    // Read symbol table node signature
    std::string signature = fread_string(fptr, 4);
 
    // Check signature to verify that we point indeed to a symbol table node
    if (std::strncmp(signature.c_str(), &table_signature[0], 4) != 0) {
        std::string msg = "No \"SNOD\" signature found at address "+
                          hexaddress(address)+" in HDF5 file. "
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_SYMBOL_TABLE_NODE, msg);
    }
 
    // Read version
    uint32_t version = fread_uint32(fptr, 1);
 
    // Check version
    if (version != 1) {
        std::string msg = "Symbol table node version "+gammalib::str(version)+
                          " is not 1. Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_SYMBOL_TABLE_NODE, msg);
    }
 
    // Read fields
    fread_zero(fptr, 1);
    int symbols = fread_int(fptr, 2);
 
    // Debugging: print symbol table node information
    #if defined(G_HDF5_DEBUG)
    std::string msg = "Symbol_table_node(address=" + hexaddress(address) +
                      ", version=" + gammalib::str(version) +
                      ", symbols=" + gammalib::str(symbols)+")";
    debug_msg(msg, indent);
    #endif
 
    // Get current file pointer position as address of first symbol table entry
    uint64_t pos = std::ftell(fptr);
 
    // Loop over symbol table node entries
    for (int i = 0; i < symbols; ++i) {
 
        // Create entry XML element and append it to node
        GXmlElement* entry = table->append("entry");
 
        // Read symbol table entry
        read_symbol_table_entry(fptr, pos, offset, heap, entry, indent+1);
 
        // Increment pointer
        pos += 2 * m_offset + 24;
 
    } // endfor: looped over symbol table node entries
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read and handle HDF5 symbol table entry
 *
 * @param[in] fptr File pointer.
 * @param[in] address Address of symbol table entry.
 * @param[in] offset Byte offset in heap.
 * @param[in] heap Heap.
 * @param[in] entry Pointer to XML element for entry.
 * @param[in] indent Indent for debugging messages.
 *
 * @exception GException::invalid_value
 *            Invalid cache type encountered
 *
 * Read and handle a symbol table entry. If the entry has a cache type of 1
 * the method will read the corresponding group.
 *
 * Each symbol table entry in a symbol table node is designed to allow for
 * very fast browsing of stored objects. Toward that design goal, the symbol
 * table entries include space for caching certain constant metadata from the
 * object header.
 ***************************************************************************/
void GHdf5::read_symbol_table_entry(FILE*              fptr,
                                    const uint64_t&    address,
                                    const uint64_t&    offset,
                                    const std::string& heap,
                                    GXmlElement*       entry,
                                    const int&         indent)
{
    // Initialise cache-type dependent information
    uint64_t address_btree     = 0;
    uint64_t address_heap      = 0;
    uint32_t offset_link_value = 0;
 
    // Move pointer to address
    std::fseek(fptr, address, SEEK_SET);
 
    // Read table entry
    uint64_t link_name_offset      = fread_uint64(fptr, m_offset);
    uint64_t object_header_address = fread_uint64(fptr, m_offset);
    uint32_t cache_type            = fread_uint32(fptr, 4);
    fread_zero(fptr, 4);
 
    // Get link name
    std::string link_name = extract(heap, offset+link_name_offset);
 
    // Set link name as attribute of XML element
    entry->attribute("name", link_name);
 
    // Decode scratch pad according to cache type
    switch (cache_type) {
 
    // No data is cached in the scratch space
    case 0:
        break;
 
    // Group object header metadata is cached in the scratch-pad space.
    // This implies that the symbol table entry refers to another group.
    case 1:
        address_btree = fread_uint64(fptr, m_offset);
        address_heap  = fread_uint64(fptr, m_offset);
        break;
 
    // The entry is a symbolic link. The first four bytes of the scratch-pad
    // space are the offset into the local heap for the link value. The object
    // header address will be undefined.
    case 2:
        offset_link_value = fread_uint32(fptr, 4);
        break;
 
    // Invalid cache type
    default:
        std::string msg = "Invalid cache type "+gammalib::str(cache_type)+
                          " encountered in symbol table entry located at "
                          "address "+hexaddress(address)+" in HDF5 file. "
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_SYMBOL_TABLE_ENTRY, msg);
        break;
    } // endswitch: decoded scratch pad
 
    // Debugging: print symbol table entry information
    #if defined(G_HDF5_DEBUG)
    std::string msg = "Symbol_table_entry(address="+hexaddress(address)+
                      ", name="+link_name+
                      ", type="+gammalib::str(cache_type)+
                      ", header_address="+hexaddress(object_header_address);
    if (cache_type == 1) {
        msg += ", btree_address="+hexaddress(address_btree)+
               ", heap_address="+hexaddress(address_heap);
    }
    else if (cache_type == 2) {
        msg += "offset_link_value="+hexaddress(offset_link_value);
    }
    msg += ")";
    debug_msg(msg, indent);
    #endif
 
    // Create symbol table header XML element and append it to entry
    GXmlElement* header = entry->append("header");
 
    // Read header
    read_object_header(fptr, object_header_address, header, indent+1);
 
    // If group metadata is cached then read group
    if (cache_type == 1) {
 
        // Create group XML element and append it to entry
        GXmlElement* group = entry->append("group");
 
        // Read group
        read_group(fptr, address_btree, address_heap, group, indent+1);
 
    } // endif: read group
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read and handle HDF5 group
 *
 * @param[in] fptr File pointer.
 * @param[in] btree_address Address of B-tree.
 * @param[in] heap_address Address of heap.
 * @param[in] group Pointer to XML element for group.
 * @param[in] indent Indent for debugging messages.
 *
 * @exception GException::invalid_value
 *            TREE or HEAP signatures not found
 *            Heap version is not zero
 *            B-tree node type is not 0 or 1
 *
 * This method only supports version 1 B-trees. Version 2 B-trees have the
 * signature "BTHD" and need to be implemented if required.
 ***************************************************************************/
void GHdf5::read_group(FILE*           fptr,
                       const uint64_t& address_btree,
                       const uint64_t& address_heap,
                       GXmlElement*    group,
                       const int&      indent)
{
    // Set TREE and HEAP signatures
    const char tree_signature[] = {'T', 'R', 'E', 'E'};
    const char heap_signature[] = {'H', 'E', 'A', 'P'};
 
    // Debugging: on entry, print group information
    #if defined(G_HDF5_DEBUG)
    std::string msg = "Group(btree_address="+hexaddress(address_btree)+
                      ", heap_address="+hexaddress(address_heap)+")";
    debug_msg(msg, indent);
    #endif
 
    // Move file pointer to heap address
    std::fseek(fptr, address_heap, SEEK_SET);
 
    // Read heap signature
    std::string signature = fread_string(fptr, 4);
 
    // Check signature to verify that we point indeed to a B-tree
    if (std::strncmp(signature.c_str(), &heap_signature[0], 4) != 0) {
        std::string msg = "No \"HEAP\" signature found at address "+
                          hexaddress(address_heap)+" in HDF5 file. "
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_GROUP, msg);
    }
 
    // Read heap version
    int version = fread_int(fptr, 1);
 
    // Check heap version
    if (version != 0) {
        std::string msg = "Heap version "+gammalib::str(version)+" is not 0."
                          " Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_GROUP, msg);
    }
 
    // Read heap fields
    fread_zero(fptr, 3);
    uint64_t heap_size = fread_uint64(fptr, m_length);
    fread_uint64(fptr, m_length);
    uint64_t heap_address = fread_uint64(fptr, m_offset);
 
    // Move file pointer to heap data
    std::fseek(fptr, heap_address, SEEK_SET);
 
    // Read heap data
    std::string heap = fread_data(fptr, heap_size);
 
    // Debugging: print heap information
    #if defined(G_HDF5_DEBUG)
    msg = "Heap(address="+hexaddress(address_heap)+
          ", data_address="+hexaddress(heap_address)+
          ", size="+gammalib::str(heap_size)+")";
    debug_msg(msg, indent+1);
    #endif
 
    // Move pointer to B-tree address
    std::fseek(fptr, address_btree, SEEK_SET);
 
    // Read B-tree signature
    signature = fread_string(fptr, 4);
 
    // Check signature to verify that we point indeed to a B-tree
    if (std::strncmp(signature.c_str(), &tree_signature[0], 4) != 0) {
        std::string msg = "No \"TREE\" signature found at address "+
                          hexaddress(address_btree)+" in HDF5 file. "
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_GROUP, msg);
    }
 
    // Read node type
    int node_type = fread_int(fptr, 1);
 
    // Check node type
    if (node_type < 0 || node_type > 1) {
        std::string msg = "B-tree node type "+gammalib::str(node_type)+" is not 0 "
                          "or 1. Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_GROUP, msg);
    }
 
    // Read fields
    int      node_level    = fread_int(fptr, 1);
    int      entries       = fread_int(fptr, 2);
    uint64_t address_left  = fread_uint64(fptr, m_offset);
    uint64_t address_right = fread_uint64(fptr, m_offset);
 
    // Debugging: print tree information
    #if defined(G_HDF5_DEBUG)
    msg = "B-tree(address="+hexaddress(address_btree)+
          ", target=";
    if (node_type == 0) {
        msg += "\"group nodes\"";
    }
    else {
        msg += "\"raw data chunk\"";
    }
    msg += ", entries="+gammalib::str(entries)+")";
    debug_msg(msg, indent+1);
    #endif
 
    // Loop over entries
    for (int i = 0; i < entries; ++i) {
 
        // Handle node type 0
        if (node_type == 0) {
 
            // Get byte offset in heap and child address
            uint64_t offset        = fread_uint64(fptr, m_length);
            uint64_t address_child = fread_uint64(fptr, m_offset);
 
            // Create table XML element and append it to group
            GXmlElement* table = group->append("table");
 
            // Read symbol table node
            read_symbol_table_node(fptr, address_child, offset, heap, table, indent+2);
 
        } // endif: node type was 0
 
        // Handle node type 1
        else {
 
        } // endelse: node type was 1
 
    } // endfor: looped over entries
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read HDF5 B-tree for chunked data
 *
 * @param[in] fptr File pointer.
 * @param[in] address Address of B-tree.
 * @param[in] dimensions Dataset dimensions + 1.
 * @param[in] btree Pointer to XML element for btree.
 *
 * @exception GException::invalid_value
 *            TREE signature not found
 *            B-tree node type is not 1
 *
 * This method only supports version 1 B-trees. Version 2 B-trees have the
 * signature "BTHD" and need to be implemented if required.
 ***************************************************************************/
void GHdf5::read_btree_chunked(FILE*           fptr,
                               const uint64_t& address,
                               const int&      dimensions,
                               GXmlElement*    btree)
{
    // Set TREE signature
    const char tree_signature[] = {'T', 'R', 'E', 'E'};
 
    // Move pointer to B-tree address
    std::fseek(fptr, address, SEEK_SET);
 
    // Read B-tree signature
    std::string signature = fread_string(fptr, 4);
 
    // Check signature to verify that we point indeed to a B-tree
    if (std::strncmp(signature.c_str(), &tree_signature[0], 4) != 0) {
        std::string msg = "No \"TREE\" signature found at address "+
                          hexaddress(address)+" in HDF5 file. "
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_BTREE_CHUNKED, msg);
    }
 
    // Read node type
    int node_type = fread_int(fptr, 1);
 
    // Check node type
    if (node_type != 1) {
        std::string msg = "B-tree node type "+gammalib::str(node_type)+
                          " for chunked data is not 1."
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_BTREE_CHUNKED, msg);
    }
 
    // Read fields
    int      node_level    = fread_int(fptr, 1);
    int      entries       = fread_int(fptr, 2);
    uint64_t address_left  = fread_uint64(fptr, m_offset);
    uint64_t address_right = fread_uint64(fptr, m_offset);
 
    // Add attributes to XML element
    btree->attribute("type", gammalib::str(node_type));
    btree->attribute("level", gammalib::str(node_level));
    btree->attribute("entries", gammalib::str(entries));
    btree->attribute("dimensions", gammalib::str(dimensions));
 
    // Loop over entries
    for (int i = 0; i < entries; ++i) {
 
        // Read fields
        int chunk_size  = fread_int(fptr, 4);
        int filter_mask = fread_int(fptr, 4);
 
        // Create key XML element and append it to B-tree
        GXmlElement* key = btree->append("key");
 
        // Add attributes to XML element
        key->attribute("size", gammalib::str(chunk_size));
        key->attribute("filter_mask", gammalib::str(filter_mask));
 
        // Read keys
        for (int k = 0; k < dimensions; ++k) {
 
            // Read fields
            uint64_t index = fread_uint64(fptr, 8);
 
            // Add attribute to XML element
            if (k < dimensions-1) {
                key->attribute("index"+gammalib::str(k), gammalib::str(index));
            }
            else {
                key->attribute("offset", gammalib::str(index));
            }
 
        } // endfor: looped over dimensions
 
        // Read child address
        uint64_t address_child = fread_uint64(fptr, m_offset);
 
        // Add attributes to XML element
        key->attribute("address", gammalib::str(address_child));
 
        // If node level is not zero this nodes points to another B-tree, hence
        // recursively analyse this subtree
        if (node_level != 0) {
 
            // Get current address
            uint64_t address = std::ftell(fptr);
 
            // Create B-tree XML element and append it to key
            GXmlElement* subbtree = key->append("btree");
 
            // Read B-tree
            read_btree_chunked(fptr, address_child, dimensions, subbtree);
 
            // Restore current address for reading of next key
            std::fseek(fptr, address, SEEK_SET);
 
        } // endif: node was not a leaf node
 
    } // endfor: looped over entries
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read object header
 *
 * @param[in] fptr File pointer.
 * @param[in] address Address of header.
 * @param[in] header Pointer to XML element for header.
 * @param[in] indent Indent for debugging messages.
 *
 * @exception GException::invalid_value
 *            Invalid object header version
 *
 * The header information of an object is designed to encompass all the
 * information about an object, except for the data itself. This information
 * includes the dataspace, datatype, information about how the data is stored
 * on disk (in external files, compressed, broken up in blocks, and so on),
 * as well as other information used by the library to speed up access to the
 * data objects or maintain a file’s integrity. Information stored by user
 * applications as attributes is also stored in the object’s header. The
 * header of each object is not necessarily located immediately prior to the
 * object’s data in the file and in fact may be located in any position in
 * the file. The order of the messages in an object header is not
 * significant.
 *
 * Object headers are composed of a prefix and a set of messages. The prefix
 * contains the information needed to interpret the messages and a small
 * amount of metadata about the object, and the messages contain the majority
 * of the metadata about the object.
 ***************************************************************************/
void GHdf5::read_object_header(FILE*           fptr,
                               const uint64_t& address,
                               GXmlElement*    header,
                               const int&      indent)
{
    // Move pointer to address
    std::fseek(fptr, address, SEEK_SET);
 
    // Read version
    int version = fread_int(fptr, 1);
 
    // Version dependent reading
    switch (version) {
    case 1:
        read_object_header_v1(fptr, header, indent);
        break;
    case 2:
        read_object_header_v2(fptr, header, indent);
        break;
    default:
        std::string msg = "Object header version "+gammalib::str(version)+" is "
                          "not comprised between 1 and 2. Please specify a "
                          "valid HDF5 file.";
        throw GException::invalid_value(G_READ_OBJECT_HEADER, msg);
        break;
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read object header version 1
 *
 * @param[in] fptr File pointer.
 * @param[in] header Pointer to XML element for header.
 * @param[in] indent Indent for debugging messages.
 *
 * @exception GException::invalid_value
 *            Invalid message type
 ***************************************************************************/
void GHdf5::read_object_header_v1(FILE*        fptr,
                                  GXmlElement* header,
                                  const int&   indent)
{
    // Read header prefix fields
    fread_zero(fptr, 1);
    int      messages  = fread_int(fptr, 2);
    int      ref_count = fread_int(fptr, 4);
    uint64_t size      = fread_int(fptr, 4);
    fread_zero(fptr, 4);
 
    // Initialise address of first message
    uint64_t address = std::ftell(fptr);
 
    // Debugging: print header information
    #if defined(G_HDF5_DEBUG)
    std::string msg = "Object_header(address="+hexaddress(address)+
                      ", version=1"
                      ", messages="+gammalib::str(messages)+
                      ", ref_count="+gammalib::str(ref_count)+
                      ", size="+gammalib::str(size)+")";
    debug_msg(msg, indent);
    #endif
 
    // Loop over all header messages
    for (int i = 0; i < messages; ++i) {
 
        // Move pointer to message address
        std::fseek(fptr, address, SEEK_SET);
 
        // Read message fields
        int msg_type  = fread_int(fptr, 2);
        int msg_size  = fread_int(fptr, 2);
        int msg_flags = fread_int(fptr, 1);
        fread_zero(fptr, 3);
 
        // Check message type
        if (msg_type < 0 || msg_type > 23) {
            std::string msg = "Invalid message type "+hexaddress(msg_type)+
                              " encountered at address "+hexaddress(address)+
                              " in HDF5 file. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_OBJECT_HEADER_V1, msg);
        }
 
        // Kludge: If message size is not a multiple of 8 this is an invalid
        // message. Set the message size to 8 Bytes and skip.
        if (msg_size % 8 != 0) {
            #if defined(G_HDF5_DEBUG)
            std::string msg = "*** invalid message size "+gammalib::str(msg_size)+
                              " for message type "+gammalib::str(msg_type)+
                              " with flags "+gammalib::str(msg_flags)+
                              " at address "+hexaddress(address)+
                              ". Skip message.";
            debug_msg(msg, indent+1);
            #endif
            msg_size = 8;
            address += msg_size + 8;
            continue;
        }
 
        // Set address to next message
        address += msg_size + 8;
 
        // Handle continuation message
        if (msg_type == 16) {
 
            // Update address and size
            address = fread_uint64(fptr, m_offset);
            size    = fread_uint64(fptr, m_length);
 
            // Debugging: print continuation message
            #if defined(G_HDF5_DEBUG)
            std::string msg = "- message(address=" + hexaddress(address) +
                              ", Continuation, size= " + gammalib::str(size) + ")";
            debug_msg(msg, indent+1);
            #endif
 
            // Continue
            continue;
 
        } // endif: handled continuation message
 
        // Create message XML element and append it to header
        GXmlElement* message = header->append("message");
 
        // Handle remaining messages
        read_object_header_message(fptr, msg_type, message, indent+1);
 
    } // endfor: looped over messages
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read object header version 2
 *
 * @param[in] fptr File pointer.
 * @param[in] header Pointer to XML element for header.
 * @param[in] indent Indent for debugging messages.
 *
 * @exception GException::feature_not_implemented
 *            Method not yet implemented.
 *
 * Method not yet implemented.
 ***************************************************************************/
void GHdf5::read_object_header_v2(FILE*        fptr,
                                  GXmlElement* header,
                                  const int&   indent)
{
    // Throw exception
    std::string msg = "Decoding of header version 2 not yet implemented.";
    throw GException::feature_not_implemented(G_READ_OBJECT_HEADER_V2, msg);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read one object header message
 *
 * @param[in] fptr File pointer.
 * @param[in] type Message type.
 * @param[in] message Pointer to XML element for message.
 * @param[in] indent Indent for debugging messages.
 *
 * @exception GException::invalid_value
 *            Invalid message
 ***************************************************************************/
void GHdf5::read_object_header_message(FILE*        fptr,
                                       const int&   type,
                                       GXmlElement* message,
                                       const int&   indent)
{
    // Debugging: initialise debug message
    #if defined(G_HDF5_DEBUG)
    std::string msg = "- message(address=" + hexaddress(std::ftell(fptr)) + ", ";
    #else
    std::string msg;
    #endif
 
    // Handle messages according to type
    switch (type) {
 
    // NIL message
    case 0: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "NIL";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "NIL");
        break;
    }
 
    // Dataspace message
    case 1: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Dataspace";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "Dataspace");
 
        // Handle dataspace message
        read_message_dataspace(fptr, message, msg);
        break;
    }
 
    // Link Info message
    case 2: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Link Info - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "LinkInfo");
        break;
    }
 
    // Datatype message
    case 3: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Datatype";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "Datatype");
 
        // Handle datatype message
        read_message_datatype(fptr, message, msg);
        break;
    }
 
    // Fill Value (old) message
    case 4: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Fill Value (old) - NOT FULLY DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "FillValueOld");
 
        // Read fields
        int size = fread_int(fptr, 4);
 
        // Read fill value
        //TODO
 
        // Add attributes to XML element
        message->attribute("size", gammalib::str(size));
 
        // Debugging: add dimensions and flags
        #if defined(G_HDF5_DEBUG)
        msg += ", size=" + gammalib::str(size);
        #endif
        break;
    }
 
    // Fill Value message
    case 5: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Fill Value - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "FillValue");
        break;
    }
 
    // Link message
    case 6: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Link - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "Link");
        break;
    }
 
    // External Data Files message
    case 7: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "External Data Files - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "ExternalDataFile");
        break;
    }
 
    // Data Storage - Layout message
    case 8: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Data Storage - Layout";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "DataStorageLayout");
 
        // Handle message
        read_message_layout(fptr, message, msg);
        break;
    }
 
    // Bogus message
    case 9: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Bogus - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "Bogus");
        break;
    }
 
    // Group Info message
    case 10: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Group Info - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "GroupInfo");
        break;
    }
 
    // Data Storage - Filter Pipeline message
    case 11: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Data Storage - Filter Pipeline";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "DataStorageFilterPipeline");
 
        // Handle message
        read_message_filter(fptr, message, msg);
        break;
    }
 
    // Attribute message
    case 12: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Attribute";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "Attribute");
 
        // Handle message
        read_message_attribute(fptr, message, msg);
        break;
    }
 
    // Object Comment message
    case 13: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Object Comment - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "ObjectComment");
        break;
    }
 
    // Object Modification Time (Old) message
    case 14: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Object Modification Time (Old) - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "ObjectModificationTimeOld");
        break;
    }
 
    // Shared Message Table message
    case 15: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Shared Message Table - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "SharedMessageTable");
        break;
    }
 
    // Object Header Continuation message
    case 16: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Continuation";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "Continuation");
        break;
    }
 
    // Symbol Table Message message
    case 17: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Symbol Table Message";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "SymbolTableMessage");
 
        // Read fields
        uint64_t address_tree = fread_uint64(fptr, m_offset);
        uint64_t address_heap = fread_uint64(fptr, m_offset);
 
        // Add attributes to XML element
        message->attribute("address_tree", gammalib::str(address_tree));
        message->attribute("address_heap", gammalib::str(address_heap));
 
        // Debugging: add dimensions and flags
        #if defined(G_HDF5_DEBUG)
        msg += ", address_tree=" + hexaddress(address_tree) +
               ", address_heap=" + hexaddress(address_heap);
        #endif
        break;
    }
 
    // Object Modification Time message
    case 18: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Object Modification Time - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "ObjectModificationTime");
        break;
    }
 
    // B-tree K Values message
    case 19: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "B-tree K Values - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "BtreeKValue");
        break;
    }
 
    // Driver Info message
    case 20: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Driver Info - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "DriverInfo");
        break;
    }
 
    // Attribute Info message
    case 21: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Attribute Info - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "AttributeInfo");
        break;
    }
 
    // Object Reference Count message
    case 22: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "Object Reference Count - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "ObjectReferenceCount");
        break;
    }
 
    // File Space Info message
    case 23: {
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += "File Space Info - NOT DECODED!";
        #endif
 
        // Add message type to XML element
        message->attribute("type", "FileSpaceInfo");
        break;
    }
 
    // Code should never arrive here
    default:
        break;
 
    } // endwsitch
 
    // Debugging: finalise debug message
    #if defined(G_HDF5_DEBUG)
    msg += ")";
    debug_msg(msg, indent);
    #endif
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read dataspace message
 *
 * @param[in] fptr File pointer.
 * @param[in] type Message type.
 * @param[in] msg String for debugging messages (not used if empty length).
 * @return Number of elements in dataspace.
 *
 * The dataspace message describes the number of dimensions (in other words,
 * "rank") and size of each dimension that the data object has. This message
 * is only used for datasets which have a simple, rectilinear, array-like
 * layout; datasets requiring a more complex layout are not yet supported.
 *
 * This method supports versions 1 and 2 of the dataspace message. If other
 * version numbers are encountered the validity flag is set to 0 and the
 * decoding is stopped.
 ***************************************************************************/
void GHdf5::read_message_dataspace(FILE*        fptr,
                                   GXmlElement* message,
                                   std::string& msg)
{
    // Initialise number of elements in dataspace
    uint64_t elements = 0;
 
    // Read version
    int version = fread_int(fptr, 1);
 
    // If version is not valid then just add version number to field and
    // set number of elements to 0. This is a kludge since the test file
    // contained an invalud version number of 24, and this version number
    // is hence filtered.
    if (version < 1 || version > 2) {
        message->attribute("version", gammalib::str(version));
        message->attribute("elements", "0");
        #if defined(G_HDF5_DEBUG)
        if (msg.length() > 0) {
            msg += ", version=" + gammalib::str(version) +
                   ", elements=0";
        }
        #endif
    }
 
    // ... otherwise handle dataspace
    else {
 
        // Read fields
        int dimensions = fread_int(fptr, 1);
        int flags      = fread_int(fptr, 1);
        int type       = fread_int(fptr, 1); // Only exists in version 2
        fread_zero(fptr, 4);
 
        // Add attributes to XML element
        message->attribute("version", gammalib::str(version));
        message->attribute("dimensions", gammalib::str(dimensions));
        message->attribute("flags",      gammalib::str(flags));
        if (version == 2) {
            message->attribute("type", gammalib::str(type));
        }
 
        // Debugging: add dimensions and flags
        #if defined(G_HDF5_DEBUG)
        if (msg.length() > 0) {
            msg += ", version=" + gammalib::str(version) +
                   ", dimensions=" + gammalib::str(dimensions) +
                   ", flags=" + gammalib::str(flags);
            if (version == 2) {
                msg += ", type=" + gammalib::str(type);
            }
        }
        #endif
 
        // If there are no dimensions then set number of elements to
        // one
        if (dimensions == 0) {
            elements = 1;
        }
 
        // ... otherwise loop over dimensions to extract sizes and compute
        // the number of elements
        else {
 
            // Loop over dimensions
            for (int i = 0; i < dimensions; ++i) {
 
                // Read dimension size
                uint64_t size = fread_uint64(fptr, m_length);
 
                // Compute number of elements
                if (i == 0) {
                    elements = int(size);
                }
                else {
                    elements *= int(size);
                }
 
                // Add size to XML element
                message->attribute("size"+gammalib::str(i), gammalib::str(size));
 
                // Debugging: add size
                #if defined(G_HDF5_DEBUG)
                if (msg.length() > 0) {
                    msg += ", size"+gammalib::str(i)+"=" + gammalib::str(size);
                }
                #endif
 
            } // endfor: looped over dimensions
 
        } // endelse: there were dimensions
 
        // Optionally extract maximum sizes
        if ((flags & 1) == 1) {
 
            // Loop over dimensions to extract maximum sizes
            for (int i = 0; i < dimensions; ++i) {
 
                // Read dimension size
                uint64_t size = fread_uint64(fptr, m_length);
 
                // Add maximum size to XML element
                message->attribute("max_size"+gammalib::str(i), gammalib::str(size));
 
                // Debugging: add maximum size
                #if defined(G_HDF5_DEBUG)
                if (msg.length() > 0) {
                    msg += ", max_size"+gammalib::str(i)+"=" + gammalib::str(size);
                }
                #endif
 
            } // endfor: looped over dimensions
 
        } // endif: maximum sizes were present
 
        // Optionally extract permutation index
        if ((version == 1) && ((flags & 2) == 2)) {
 
            // Loop over dimensions to extract permutation index
            for (int i = 0; i < dimensions; ++i) {
 
                // Read dimension size
                uint64_t size = fread_uint64(fptr, m_length);
 
                // Add permutation to XML element
                message->attribute("permutation"+gammalib::str(i), gammalib::str(size));
 
                // Debugging: add permutation
                #if defined(G_HDF5_DEBUG)
                if (msg.length() > 0) {
                    msg += ", permutation"+gammalib::str(i)+"=" + gammalib::str(size);
                }
                #endif
 
            } // endfor: looped over dimensions
 
        } // endif: permutation indices were present
 
        // Add elements to XML element
        message->attribute("elements", gammalib::str(elements));
 
        // Debugging: add elements
        #if defined(G_HDF5_DEBUG)
        if (msg.length() > 0) {
            msg += ", elements=" + gammalib::str(elements);
        }
        #endif
 
    } // endelse: version was valid
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read datatype message
 *
 * @param[in] fptr File pointer.
 * @param[in] message Pointer to XML element.
 * @param[in] msg String for debugging messages (not used if empty length).
 * @return Size of datatype elements in Bytes.
 *
 * @exception GException::invalid_value
 *            Invalid datatype version or class
 *
 * The datatype message defines the datatype for each element of a dataset
 * or a common datatype for sharing between multiple datasets. A datatype
 * can describe an atomic type like a fixed- or floating-point type or more
 * complex types like a C struct (compound datatype), array (array datatype)
 * or C++ vector (variable-length datatype).
 *
 * Datatype messages that are part of a dataset object do not describe how
 * elements are related to one another; the dataspace message is used for
 * that purpose. Datatype messages that are part of a committed datatype
 * (formerly named datatype) message describe a common datatype that can be
 * shared by multiple datasets in the file.
 ***************************************************************************/
void GHdf5::read_message_datatype(FILE*        fptr,
                                  GXmlElement* message,
                                  std::string& msg)
{
    // Read fields
    int class_version = fread_int(fptr, 1);
    int bits1         = fread_int(fptr, 1);
    int bits2         = fread_int(fptr, 1);
    int bits3         = fread_int(fptr, 1);
    int size          = fread_int(fptr, 4); // Size of datatype elements in Bytes
 
    // Extract version and class
    int version = class_version / 16;
    int cls     = class_version & 15;
 
    // Check version
    if (version < 1 || version > 5) {
        std::string msg = "Datatype version "+gammalib::str(version)+
                          " is not between 1 and 5."
                          " Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_MESSAGE_DATATYPE, msg);
    }
 
    // Check class
    if (cls < 0 || cls > 11) {
        std::string msg = "Datatype class "+gammalib::str(cls)+
                          " is not between 0 and 11."
                          " Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_MESSAGE_DATATYPE, msg);
    }
 
    // Get classname
    std::string clsname = gammalib::hdf5::classname(cls);
 
    // Add version, class and size to XML element
    message->attribute("version",   gammalib::str(version));
    message->attribute("class",     gammalib::str(cls));
    message->attribute("classname", clsname);
    message->attribute("size",      gammalib::str(size));
 
    // Debugging: add version, class and size
    #if defined(G_HDF5_DEBUG)
    if (msg.length() > 0) {
        msg += ", version=" + gammalib::str(version) +
               ", class=" + gammalib::str(cls) +
               ", classname=\"" + clsname + "\"" +
               ", size=" + gammalib::str(size);
    }
    #endif
 
    // Handle class specific information
    switch (cls) {
 
    // Fixed point
    case 0: {
        // Decode Bit fields
        int order = bits1 & 1;       // 0=Little-endian, 1=Big-endian
        int lopad = (bits1 & 2) / 2; // Lo-pad
        int hipad = (bits1 & 4) / 4; // High-pad
        int sign  = (bits1 & 8) / 8; // 0=Unsigned, 1=Signed
 
        // Read properties
        int offset    = fread_int(fptr, 2);
        int precision = fread_int(fptr, 2);
 
        // Add XML attributes
        message->attribute("order",     gammalib::str(order));
        if (order == 0) {
            message->attribute("ordername", "little-endian");
        }
        else {
            message->attribute("ordername", "big-endian");
        }
        message->attribute("signed",    gammalib::str(sign));
        message->attribute("offset",    gammalib::str(offset));
        message->attribute("precision", gammalib::str(precision));
        message->attribute("lopad",     gammalib::str(lopad));
        message->attribute("hipad",     gammalib::str(hipad));
        break;
    }
 
    // Floating point
    case 1: {
        break;
    }
 
    // Time
    case 2: {
        break;
    }
 
    // String
    case 3: {
        break;
    }
 
    // Bit field
    case 4: {
        break;
    }
 
    // Opaque
    case 5: {
        break;
    }
 
    // Compound
    case 6: {
        break;
    }
 
    // Reference
    case 7: {
        break;
    }
 
    // Enumerated
    case 8: {
        // Decode Bit fields
        int number = bits1 + bits2 * 256;
 
        // Add XML attributes
        message->attribute("number", gammalib::str(number));
 
        // Debugging: add attributes
        #if defined(G_HDF5_DEBUG)
        if (msg.length() > 0) {
            msg += ", number=" + gammalib::str(number);
        }
        #endif
 
        // Version 1 or 2
        if ((version == 1) || (version == 2)) {
 
            // Read base type datatype
            GXmlElement* basetype = message->append("basetype");
            std::string  empty;
            read_message_datatype(fptr, basetype, empty);
 
            // Allocate string vectors for name and values
            std::vector<std::string> names;
            std::vector<std::string> values;
 
            // Read names
            for (int i = 0; i < number; ++i) {
 
                // Read multiples of 8 Bytes until a NULL character is
                // encountered
                std::string name;
                bool        looping = true;
                while (looping) {
 
                    // Read next 8 Bytes
                    std::string segment = fread_data(fptr, 8);
 
                    // Search for termination NULL
                    for (int k = 0; k < 8; ++k) {
                        if (segment[k] == '\0') {            // found
                            if (k > 0) {                     // string is not empty
                                name += segment.substr(0,k); // then add substring
                            }
                            looping = false;                 // stop looping
                            break;                           // exit loop
                        }
                    }
 
                    // If we are still looping no termination NULL was found.
                    // We hence add the entire segment
                    if (looping) {
                        name += segment;
                    }
 
                } // endwhile: looping over 8 Byte junks
 
                // Append name
                names.push_back(name);
 
            } // endfor: looped over names
 
            // Read values
            for (int i = 0; i < number; ++i) {
 
                // Read value
                std::string value = fread_data_as_string(fptr, basetype);
 
                // Append value
                values.push_back(value);
 
            } // endfor: looped over values
 
            // Create XML elements for all name/value pairs
            for (int i = 0; i < number; ++i) {
 
                // Create name XML element
                GXmlElement* enumeration = message->append("enumeration");
 
                // Set name and value
                enumeration->attribute("name",  names[i]);
                enumeration->attribute("value", values[i]);
 
            } // endfor: looped over name/value pairs
 
        } // endif: version 1 or 2
        break;
    }
 
    // Variable-Length
    case 9: {
        // Decode Bit fields
        int type    = bits1 & 15; // 0=Sequence, 1=String
        int padding = bits1 / 16; // 0=Null terminated, 1=Null pad, 2=Space pad
        int charset = bits2 & 15; // 0=ASCII, 1=UTF-8
 
        // Read parent type datatype (not used for string)
        GXmlElement* parenttype = message->append("parenttype");
        std::string  empty;
        read_message_datatype(fptr, parenttype, empty);
 
        // Add XML attributes
        message->attribute("type",    gammalib::str(type));
        message->attribute("padding", gammalib::str(padding));
        message->attribute("charset", gammalib::str(charset));
 
        // Debugging: add version, class and size
        #if defined(G_HDF5_DEBUG)
        if (msg.length() > 0) {
            msg += ", type=" + gammalib::str(type) +
                   ", padding=" + gammalib::str(padding) +
                   ", charset=" + gammalib::str(charset);
        }
        #endif
        break;
    }
 
    // Array
    case 10: {
        break;
    }
 
    // Complex
    case 11: {
        break;
    }
 
    // This should never be reacked
    default:
        break;
 
    } // endwsitch: Handled class specific information
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read Data Storage - Layout message
 *
 * @param[in] fptr File pointer.
 * @param[in] message Pointer to XML element.
 * @param[in] msg String for debugging messages.
 *
 * @exception GException::invalid_value
 *            Invalid data storage layout version or class
 ***************************************************************************/
void GHdf5::read_message_layout(FILE*        fptr,
                                GXmlElement* message,
                                std::string& msg)
{
    // Read version
    int version = fread_int(fptr, 1);
 
    // Version dependent treatment
    switch (version) {
 
    // Version 1 or 2
    case 1:
    case 2: {
        break;
    }
 
    case 3: {
        // Read layout class
        int cls = fread_int(fptr, 1);
 
        // Add attribute to XML element
        message->attribute("class", gammalib::str(cls));
 
        // Debugging: add class
        #if defined(G_HDF5_DEBUG)
        msg += ", class=" + gammalib::str(cls);
        #endif
 
        // Class dependent treatment
        switch (cls) {
 
        // Compact layout
        case 0: {
            // Add attribute to XML element
            message->attribute("classname", "compact");
 
            // Debugging: add layout class
            #if defined(G_HDF5_DEBUG)
            msg += ", compact";
            #endif
            break;
        }
 
        // Contiguous layout
        case 1: {
            // Add attribute to XML element
            message->attribute("classname", "contiguous");
 
            // Debugging: add layout class
            #if defined(G_HDF5_DEBUG)
            msg += ", contiguous";
            #endif
 
            // Read fields
            uint64_t address = fread_uint64(fptr, m_offset);
            uint64_t size    = fread_uint64(fptr, m_length);
 
            // Add attributes to XML element
            message->attribute("address", gammalib::str(address));
            message->attribute("size", gammalib::str(size));
 
            // Debugging: add address and size
            #if defined(G_HDF5_DEBUG)
            msg += ", address=" + hexaddress(address) +
                   ", size=" + gammalib::str(size);
            #endif
            break;
        }
 
        // Chunked layout
        case 2: {
            // Add attribute to XML element
            message->attribute("classname", "chunked");
 
            // Debugging: add layout class
            #if defined(G_HDF5_DEBUG)
            msg += ", chunked";
            #endif
 
            // Read fields
            int           dimensions = fread_int(fptr, 1);
            uint64_t address    = fread_uint64(fptr, m_offset);
 
            // Add attributes to XML element
            message->attribute("address", gammalib::str(address));
            message->attribute("dimensions", gammalib::str(dimensions));
 
            // Debugging: add address and dimension
            #if defined(G_HDF5_DEBUG)
            msg += ", address=" + hexaddress(address) +
                   ", dimensions=" + gammalib::str(dimensions);
            #endif
 
            // Initialise number of elements
            int elements = 0;
 
            // Loop over dimensions
            for (int i = 0; i < dimensions; ++i) {
 
                // Read size
                int size = fread_int(fptr, 4);
 
                // Update number of elements
                if (i == 0) {
                    elements = size;
                }
                else {
                    elements *= size;
                }
 
                // Add attribute to XML element
                message->attribute("size"+gammalib::str(i), gammalib::str(size));
 
                // Debugging: add size
                #if defined(G_HDF5_DEBUG)
                msg += ", size" + gammalib::str(i) +
                       "=" + gammalib::str(size);
                #endif
 
            } // endfor: looped over dimensions
 
            // Add attributes to XML element
            message->attribute("elements", gammalib::str(elements));
 
            // Debugging: add elements
            #if defined(G_HDF5_DEBUG)
            msg += ", elements=" + gammalib::str(elements);
            #endif
 
            // Create XML element for B-tree that holds addresses of chunked
            // data
            GXmlElement* btree = message->append("btree");
 
            // Read B-tree
            read_btree_chunked(fptr, address, dimensions, btree);
            break;
        }
 
        // Invalid layout class
        default: {
            std::string msg = "Data storage layout class "+gammalib::str(cls)+
                              " is not between 0 and 2."
                              " Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_READ_MESSAGE_LAYOUT, msg);
            break;
        }
        } // endswitch: class
        break;
    }
 
    // Invalid version
    default: {
        std::string msg = "Data storage layout version "+gammalib::str(version)+
                          " is not between 1 and 3."
                          " Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_MESSAGE_LAYOUT, msg);
        break;
    }
    } // endswitch: version
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read Data Storage - Filter pipeline message
 *
 * @param[in] fptr File pointer.
 * @param[in] type Message type.
 * @param[in] msg String for debugging messages.
 *
 * @exception GException::invalid_value
 *            Invalid data filter pipeline message version
 ***************************************************************************/
void GHdf5::read_message_filter(FILE*        fptr,
                                GXmlElement* message,
                                std::string& msg)
{
    // Read version
    int version = fread_int(fptr, 1);
 
    // Version dependent treatment
    switch (version) {
 
    // Version 1
    case 1: {
 
        // Read fields
        int filters = fread_int(fptr, 1);
        fread_zero(fptr, 6);
 
        // Add attribute to XML element
        message->attribute("filters", gammalib::str(filters));
 
        // Debugging: add filters
        #if defined(G_HDF5_DEBUG)
        msg += ", filters=" + gammalib::str(filters);
        #endif
 
        // Loop over filter descriptions
        for (int i = 0; i < filters; ++i) {
 
            // Add filter XML element
            GXmlElement* filter = message->append("filter");
 
            // Read fields
            int id               = fread_int(fptr, 2);
            int name_length      = fread_int(fptr, 2);
            int flags            = fread_int(fptr, 2);
            int number_of_values = fread_int(fptr, 2);
 
            // Read filter name
            std::string name;
            if (name_length > 0) {
                name = fread_string(fptr, name_length);
            }
 
            // Add attributes to XML element
            filter->attribute("name", name);
            filter->attribute("id", gammalib::str(id));
            filter->attribute("flags", gammalib::str(flags));
 
            // Loop over client data values
            for (int k = 0; k < number_of_values; ++k) {
 
                // Read client data
                int client_data = fread_int(fptr, 4);
 
                // Create XML element to homd client data
                GXmlElement* data = filter->append("data");
 
                // Create XML text node
                GXmlText text(gammalib::str(client_data));
 
                // Append text node
                data->append(text);
 
            } // endfor: looped over client data
 
        } // endfor: looped over filters
        break;
    }
 
    // Version 2
    case 2: {
        break;
    }
 
    // Invalid version
    default: {
        std::string msg = "Data storage filter pipeline version " +
                          gammalib::str(version) + " is not one of 1 or 2."
                          " Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_MESSAGE_FILTER, msg);
        break;
    }
    } // endswitch: version
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read Attribute message
 *
 * @param[in] fptr File pointer.
 * @param[in] type Message type.
 * @param[in] msg String for debugging messages.
 *
 * @exception GException::invalid_value
 *            Invalid attribute version
 *
 * The Attribute message is used to store objects in the HDF5 file which are
 * used as attributes, or "metadata" about the current object. An attribute
 * is a small dataset; it has a name, a datatype, a dataspace, and raw data.
 *
 * Since attributes are stored in the object header, they should be
 * relatively small (in other words, less than 64KB). They can be associated
 * with any type of object which has an object header (groups, datasets, or
 * committed (named) datatypes).
 *
 * Note: Attributes on an object must have unique names: the HDF5 Library
 * currently enforces this by causing the creation of an attribute with a
 * duplicate name to fail. Attributes on different objects may have the same
 * name, however.
 ***************************************************************************/
void GHdf5::read_message_attribute(FILE*        fptr,
                                   GXmlElement* message,
                                   std::string& msg)
{
    // Read version
    int version = fread_int(fptr, 1);
 
    // Version dependent treatment
    switch (version) {
 
    // Version 1
    case 1: {
        // Read fields
        fread_zero(fptr, 1);
        int name_size      = pad2eight(fread_int(fptr, 2));
        int datatype_size  = pad2eight(fread_int(fptr, 2));
        int dataspace_size = pad2eight(fread_int(fptr, 2));
 
        // Read attribute name
        std::string name = fread_string(fptr, name_size);
 
        // Add attributes to XML element
        message->attribute("name", name);
 
        // Debugging: add message name
        #if defined(G_HDF5_DEBUG)
        msg += ", name=" + name;
        #endif
 
        // Read datatype
        uint64_t     address  = std::ftell(fptr);
        GXmlElement* datatype = message->append("datatype");
        read_message_datatype(fptr, datatype, msg);
        std::fseek(fptr, address+datatype_size, SEEK_SET);
 
        // Read dataspace
        address                = std::ftell(fptr);
        GXmlElement* dataspace = message->append("dataspace");
        read_message_dataspace(fptr, dataspace, msg);
        std::fseek(fptr, address+dataspace_size, SEEK_SET);
 
        // Optionally read message data
        if (datatype->has_attribute("size") && dataspace->has_attribute("elements")) {
 
            // Extract data element size and number of elements
            int size     = gammalib::toint(datatype->attribute("size"));
            int cls      = gammalib::toint(datatype->attribute("class"));
            int sign     = gammalib::toint(datatype->attribute("signed"));
            int elements = gammalib::toint(dataspace->attribute("elements"));
 
            // Special handling for variable-length data
            if (cls == 9) {
 
                // Get variable-length type
                int type = gammalib::toint(datatype->attribute("type"));
 
                // Read fields (sequence length, collection address, object index)
                uint32_t sequence_length    = fread_uint32(fptr, 4);
                uint64_t collection_address = fread_uint64(fptr, m_offset);
                int      object_index       = fread_int(fptr, 4);
 
                // Add attributes to XML element
                message->attribute("sequence_length", gammalib::str(sequence_length));
                message->attribute("collection_address", gammalib::str(collection_address));
                message->attribute("object_index", gammalib::str(object_index));
 
                // Debugging: add attributes
                #if defined(G_HDF5_DEBUG)
                msg += ", sequence_length="+gammalib::str(sequence_length) +
                       ", collection_address="+hexaddress(collection_address) +
                       ", object_index="+gammalib::str(object_index);
                #endif
 
                // If type is string then get attribute string from global heap
                if (type == 1) {
 
                    // Get attribute string
                    std::string value = global_heap_string(fptr,
                                                           collection_address,
                                                           object_index,
                                                           sequence_length);
 
                    // Add value to XML element
                    message->attribute("value", value);
 
                    // Debugging: add value
                    #if defined(G_HDF5_DEBUG)
                    msg += ", value=\""+value+"\"";
                    #endif
 
                } // endif: type was string
                break;
 
            } // endif: handled variable-length data
 
            // ... otherwise handle remaining classes
            else {
 
                // Compute data size
                int data_size = size * elements;
 
                // If we have data that decode it
                if (data_size > 0) {
 
                    // If there is a single element then set it as value attribute
                    if (elements == 1) {
 
                        // Read value as string
                        std::string value = fread_data_as_string(fptr, datatype);
 
                        // Set attribute
                        message->attribute("value", value);
 
                    } // endif: there was a single attribute
 
                    // ... otherwise add value tags
                    else {
 
                        // Loop over elements
                        for (int i = 0; i < elements; ++i) {
 
                            // Create value XML element
                            GXmlElement* value = message->append("value");
 
                            // Read value as string
                            std::string string = fread_data_as_string(fptr, datatype);
 
                            // Create XML text node
                            GXmlText text(string);
 
                            // Append text node
                            value->append(text);
 
                        } // endfor: looped over elements
 
                    } // endelse: added value tags
 
                } // endif: there were data
 
            } // endelse: handle remaining classes
 
        } // endif: relevant attributes present
        break;
    }
 
    // Version 2
    case 2: {
        break;
    }
 
    // Version 3
    case 3: {
        break;
    }
 
    // Invalid version
    default: {
        std::string msg = "Attribute version "+gammalib::str(version)+
                          " is not between 1 and 3."
                          " Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_READ_MESSAGE_ATTRIBUTE, msg);
        break;
    }
 
    } // endswitch: version
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Get string from global heap
 *
 * @param[in] fptr File pointer.
 * @param[in] address Collection address (address of global heap).
 * @param[in] id Heap object index.
 * @param[in] length Length of string.
 * @return String from global heap.
 *
 * @exception GException::invalid_value
 *            Invalid group heap signature
 *            Invalid group heap version
 *            Object index not found in global heap collection
 *            Length of string does not correspond to length of object
 *
 * Get string from the global heap using it's collection address @p address
 * and the heap object index @p id.
 *
 * If the requested global heap does not yet exist in the XML object it is
 * read from the file and stored into a "gobal_heap" element. Multiple global
 * heap elements may exist and are identified by their file address.
 ***************************************************************************/
std::string GHdf5::global_heap_string(FILE*           fptr,
                                      const uint64_t& address,
                                      const int&      id,
                                      const uint32_t& length)
{
    // Set global heap signature
    const char global_heap_signature[] = {'G', 'C', 'O', 'L'};
 
    // Initialise string
    std::string string;
 
    // Convert address into string
    std::string str_address = gammalib::str(address);
 
    // Get "ROOT" entry in file
    GXmlElement* root = m_xml.element("entry");
 
    // Initialise pointer to group heap
    GXmlElement* global_heap = NULL;
 
    // Search all global heaps for the one with the appropriate address
    int elements = root->elements("global_heap");
    for (int i = 0; i < elements; ++i) {
        if (root->element("global_heap", i)->attribute("address") == str_address) {
            global_heap = root->element("global_heap", i);
        }
    }
 
    // If no group heap with appropriate address was found then create one
    // and read collection information from file
    if (global_heap == NULL) {
 
        // Create group heap
        global_heap = root->append("global_heap");
 
        // Store address as attribute
        global_heap->attribute("address", str_address);
 
        // Get current file pointer position
        uint64_t current_pos = std::ftell(fptr);
 
        // Move file pointer position to group heap
        std::fseek(fptr, address, SEEK_SET);
 
        // Read group heap signature
        std::string signature = fread_string(fptr, 4);
 
        // Check signature to verify that we point indeed to a group heap
        if (std::strncmp(signature.c_str(), &global_heap_signature[0], 4) != 0) {
            std::string msg = "No \"GCOL\" signature found at address "+
                              hexaddress(address)+" in HDF5 file. "
                              "Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_GLOBAL_HEAP_STRING, msg);
        }
 
        // Read version
        uint32_t version = fread_uint32(fptr, 1);
 
        // Check version
        if (version != 1) {
            std::string msg = "Global heap version "+gammalib::str(version)+
                              " is not 1. Please specify a valid HDF5 file.";
            throw GException::invalid_value(G_GLOBAL_HEAP_STRING, msg);
        }
 
        // Read fields
        fread_zero(fptr, 3);
        uint64_t size = fread_uint64(fptr, m_length);
 
        // Store attributes
        global_heap->attribute("version", gammalib::str(version));
        global_heap->attribute("size", gammalib::str(size));
 
        // Set object header size
        int object_header_size = 8 + m_length;
 
        // Loop over global heap collection without knowing a priori how many
        // objects are present in the heap
        int read = 0;
        while (true) {
 
            // Compute remaining Bytes
            int remaining = size - read;
 
            // If there is not enough space for an object header then
            // break
            if (remaining < object_header_size) {
                break;
            }
 
            // Read object header
            int index     = fread_int(fptr, 2);
            int ref_count = fread_int(fptr, 2);
            fread_zero(fptr, 4);
            uint64_t object_size = fread_uint64(fptr, m_length);
            read += object_header_size;
 
            // If the object header index is zero we reached the end
            // and break
            if (index == 0) {
                break;
            }
 
            // Read object data as string
            std::string value = fread_string(fptr, object_size);
            read += object_size;
 
            // Read additional storage space to round up to 8 bytes
            int padded = pad2eight(object_size) - object_size;
            if (padded > 0) {
                fread_data(fptr, padded);
                read += padded;
            }
 
            // Create object element
            GXmlElement* object = global_heap->append("object");
 
            // Store attributes
            object->attribute("index", gammalib::str(index));
            object->attribute("ref_count", gammalib::str(ref_count));
            object->attribute("object_size", gammalib::str(object_size));
            object->attribute("value", value);
 
        } // endwhile: looped over global heap collection
 
        // Recover file pointer position prior to reading group heap
        std::fseek(fptr, current_pos, SEEK_SET);
 
    } // endif: created global heap and filled it with collection information
 
    // Get number of objects in global heap
    int objects = global_heap->elements("object");
 
    // Initialise found flag
    bool found = false;
 
    // Loop over objects
    for (int i = 0; i < objects; ++i) {
 
        // Get pointer to object
        GXmlElement* object = global_heap->element("object", i);
 
        // If the heap object index was found then return string
        if (gammalib::toint(object->attribute("index")) == id) {
 
            // Check that length is consistent
            int object_length = gammalib::toint(object->attribute("object_size"));
            if (object_length != length) {
                std::string msg = "Length "+object->attribute("object_size")+
                                  " of object with index "+gammalib::str(id)+
                                  " in global heap at address "+
                                  hexaddress(address)+" in HDF5 file does "
                                  "not correspond to the expected length "+
                                  gammalib::str(length)+". "
                                  "Please specify a valid HDF5 file.";
                throw GException::invalid_value(G_GLOBAL_HEAP_STRING, msg);
            }
 
            // Extract string
            string = object->attribute("value");
 
            // Signal that string was found
            found = true;
 
            // Exit search loop
            break;
        }
 
    } // endfor: looped over collection objects
 
    // If object index was not found than throw an exception
    if (!found) {
        std::string msg = "No object with index "+gammalib::str(id)+
                          " found in global heap at address "+
                          hexaddress(address)+" in HDF5 file. "
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_GLOBAL_HEAP_STRING, msg);
    }
 
    // Return string
    return string;
}
 
 
/***********************************************************************//**
 * @brief Show debug message
 *
 * @param[in] msg Debug message.
 * @param[in] indent Number of spaces to indent
 ***************************************************************************/
void GHdf5::debug_msg(const std::string& msg, const int& indent)
{
    // Indent
    for (int i = 0; i < indent; ++i) {
        std::cout << " ";
    }
 
    // Print message
    std::cout << msg << std::endl;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Format address for debugging
 *
 * @param[in] address Address.
 * @return String with address hexadecimal format.
 *
 * Returns a string with the @p address in hexadecimal format.
 ***************************************************************************/
std::string GHdf5::hexaddress(const uint64_t& address)
{
    // Set address string
    std::string string = "0x" + gammalib::str(address, "%x");
 
    // Return string
    return string;
}
 
 
/***********************************************************************//**
 * @brief Extract string from position until next null termination
 *
 * @param[in] string String.
 * @param[in] pos Start position in string.
 * @return Extracted string.
 *
 * Returns extracted string.
 ***************************************************************************/
std::string GHdf5::extract(const std::string& string, const int& pos)
{
    // Initialise result
    std::string result = "";
 
    // Get maximum string position
    int n = string.length();
 
    // Loop over string until either a \0 or the end of the string is reached
    for (int i = pos; i < n; ++i) {
        if (string[i] != '\0') {
            result += string[i];
        }
        else {
            break;
        }
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Pad number of Bytes to 8 Bytes
 *
 * @param[in] nbytes Number of Bytes.
 * @return Padded number of Bytes to 8 Bytes.
 *
 * Returns a padded number of Bytes so that @p nbytes becomes a multiple of
 * 8 Bytes.
 ***************************************************************************/
int GHdf5::pad2eight(const int& nbytes)
{
    // Initialise result
    int result = nbytes;
 
    // If nbytes is not a multiple of 8 then pad it to a multiple of 8
    if (result % 8 != 0) {
        result = ((result / 8) + 1) * 8;
    }
 
    // Return result
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                       Public HDF5 support functions                     =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Read data from HDF5 file
 *
 * @param[in] fptr File pointer.
 * @param[in] nbytes Number of Bytes to read (>0).
 * @returns String holding the data that was read.
 *
 * @exception GException::invalid_argument
 *            Number of Bytes not positive
 * @exception GException::runtime_error
 *            Attempt to read beyond end of file
 *
 * Reads @p nbytes from the current position of the file pointer @p fptr into
 * a string.
 ***************************************************************************/
std::string gammalib::hdf5::fread_data(FILE* fptr, const int& nbytes)
{
    // Check that nbytes is not zero
    if (nbytes < 1) {
        std::string msg = "Non positive number "+gammalib::str(nbytes)+
                          " of Bytes specified. Please specify a number "
                          "larger than 0.";
        throw GException::invalid_argument(G_FREAD_DATA1, msg);
    }
 
    // Initialise string
    std::string string;
 
    // Allocate string space
    string.resize(nbytes);
 
    // Read Bytes from file
    const std::size_t n = std::fread(&string[0], 1, nbytes, fptr);
 
    // Check number of Bytes read
    if (n != nbytes) {
        std::string msg = "Attempt to read beyond the end of file. Read "
                          "only "+gammalib::str(n)+" of "+gammalib::str(nbytes)+
                          " Bytes from file.";
        throw GException::runtime_error(G_FREAD_DATA1, msg);
    }
 
    // Return string
    return (string);
}
 
 
/***********************************************************************//**
 * @brief Read data
 *
 * @param[in] fptr File pointer.
 * @param[in] dataspace Pointer to "Dataspace" message XML element.
 * @param[in] datatype Pointer to "Datatype" message XML element.
 * @param[in] datalayout Pointer to "DataStorageLayout" message XML element.
 * @param[in] datafilter Pointer to "DataStorageFilterPipeline" message XML element.
 * @return String holding the data.
 *
 * @exception GException::invalid_argument
 *            No "elements" attribute in dataspace
 *            No "size" attribute in datatype
 *            No "address" attribute in datalayout
 *            No "class" attribute in datalayout
 * @excaption GException::runtime_error
 *            Attempt to read beyond end of file
 *
 * Reads entire data structure and stores the Bytes into a string.
 *
 * The address from where to read the data is extracted from @p datalayout.
 * The current position of the file pointer @p fptr is not modified by the
 * function. The number of Bytes to read is extracted from @p dataspace and
 * @p datatype.
 *
 * If @p datafilter is not NULL, a data filter is applied to the data before
 * returning them to the client.
 ***************************************************************************/
std::string gammalib::hdf5::fread_data(FILE*              fptr,
                                       const GXmlElement* dataspace,
                                       const GXmlElement* datatype,
                                       const GXmlElement* datalayout,
                                       const GXmlElement* datafilter)
{
    // Initialise string
    std::string string;
 
    // Check for required attributes
    if (!dataspace->has_attribute("elements")) {
        std::string msg = "\"Dataspace\" message is missing \"elements\""
                          " attribute";
        throw GException::invalid_argument(G_FREAD_DATA2, msg);
    }
    if (!datatype->has_attribute("size")) {
        std::string msg = "\"Datatype\" message is missing \"size\""
                          " attribute";
        throw GException::invalid_argument(G_FREAD_DATA2, msg);
    }
    if (!datalayout->has_attribute("address")) {
        std::string msg = "\"DataStorageLayout\" message is missing \"address\""
                          " attribute";
        throw GException::invalid_argument(G_FREAD_DATA2, msg);
    }
    if (!datalayout->has_attribute("class")) {
        std::string msg = "\"DataStorageLayout\" message is missing \"class\""
                          " attribute";
        throw GException::invalid_argument(G_FREAD_DATA2, msg);
    }
 
    // Handle chunked data layout
    if (datalayout->attribute("class") == "2") {
 
        //TODO: Implement reading of chunked data
 
    }
 
    // ... handle other layouts
    else {
 
        // Get number of elements in dataspace
        int elements = gammalib::toint(dataspace->attribute("elements"));
 
        // Continue only if there are elements to read
        if (elements > 0) {
 
            // Get size of each element in Bytes
            int size = gammalib::toint(datatype->attribute("size"));
 
            // Continue only if the size if positive
            if (size > 0) {
 
                // Get address from where data should be read
                uint64_t address = gammalib::toulonglong(datalayout->attribute("address"));
 
                // Determine number of Bytes to read
                int nbytes = elements * size;
 
                // Allocate string space
                string.resize(nbytes);
 
                // Get current file pointer address
                uint64_t current_address = std::ftell(fptr);
 
                // Move file pointer to address
                std::fseek(fptr, address, SEEK_SET);
 
                // Read Bytes from file
                const std::size_t n = std::fread(&string[0], 1, nbytes, fptr);
 
                // Check number of Bytes read
                if (n != nbytes) {
                    std::string msg = "Attempt to read beyond the end of file. Read "
                                      "only "+gammalib::str(n)+" of "+gammalib::str(nbytes)+
                                      " Bytes from file.";
                    throw GException::runtime_error(G_FREAD_DATA2, msg);
                }
 
                // Restore file pointer to initial address
                std::fseek(fptr, current_address, SEEK_SET);
 
            } // endif: element size was positive
 
        } // endif: there were elements to read
 
        // Optionally apply data filter
        if (datafilter != NULL) {
            string = data_filter(string, datafilter);
        }
 
    } // endelse: handles non-chunked layouts
 
    // Return string
    return string;
}
 
 
/***********************************************************************//**
 * @brief Read data chunk
 *
 * @param[in] fptr File pointer.
 * @param[in] chunk Pointer to "key" XML element.
 * @param[in] dataspace Pointer to "Dataspace" message XML element.
 * @param[in] datatype Pointer to "Datatype" message XML element.
 * @param[in] datalayout Pointer to "DataStorageLayout" message XML element.
 * @param[in] datafilter Pointer to "DataStorageFilterPipeline" message XML element.
 * @return String holding the data chunk.
 *
 * @exception GException::invalid_argument
 *            No "address" attribute in chunk key
 *            No "size" attribute in chunk key
 *            No "class" attribute in datalayout
 *            "DataStorageLayout" is not chunked
 * @exception GException::runtime_error
 *            Attempt to read beyond end of file
 *
 * Reads a data chunk and stores the Bytes into a string.
 *
 * The address from where to read the data and the number of Bytes to read
 * is extracted from @p chunk.
 *
 * The current position of the file pointer @p fptr is not modified by the
 * function.
 *
 * If @p datafilter is not NULL, a data filter is applied to the data before
 * returning them to the client.
 ***************************************************************************/
std::string gammalib::hdf5::fread_data_chunk(FILE*              fptr,
                                             const GXmlElement* chunk,
                                             const GXmlElement* dataspace,
                                             const GXmlElement* datatype,
                                             const GXmlElement* datalayout,
                                             const GXmlElement* datafilter)
{
    // Initialise string
    std::string string;
 
    // Check for required attributes
    if (!chunk->has_attribute("address")) {
        std::string msg = "Chunk key message is missing \"address\""
                          " attribute";
        throw GException::invalid_argument(G_FREAD_DATA_CHUNK, msg);
    }
    if (!chunk->has_attribute("size")) {
        std::string msg = "Chunk key message is missing \"size\""
                          " attribute";
        throw GException::invalid_argument(G_FREAD_DATA_CHUNK, msg);
    }
    if (!datalayout->has_attribute("class")) {
        std::string msg = "\"DataStorageLayout\" message is missing \"class\""
                          " attribute";
        throw GException::invalid_argument(G_FREAD_DATA_CHUNK, msg);
    }
 
    // Check that ataStorageLayout is chunked
    if (datalayout->attribute("class") != "2") {
        std::string msg = "\"DataStorageLayout\" is " +
                          datalayout->attribute("class") + " which is not "
                          "chunked. Please specify a chunked layout.";
        throw GException::invalid_argument(G_FREAD_DATA_CHUNK, msg);
    }
 
    // Get address and size
    uint64_t address = gammalib::toulonglong(chunk->attribute("address"));
    int           size    = gammalib::toint(chunk->attribute("size"));
 
    // Allocate string space
    string.resize(size);
 
    // Get current file pointer address
    uint64_t current_address = std::ftell(fptr);
 
    // Move file pointer to address
    std::fseek(fptr, address, SEEK_SET);
 
    // Read Bytes from file
    const std::size_t n = std::fread(&string[0], 1, size, fptr);
 
    // Check number of Bytes read
    if (n != size) {
        std::string msg = "Attempt to read beyond the end of file. Read "
                          "only "+gammalib::str(n)+" of "+gammalib::str(size)+
                          " Bytes from file.";
        throw GException::runtime_error(G_FREAD_DATA_CHUNK, msg);
    }
 
    // Restore file pointer to initial address
    std::fseek(fptr, current_address, SEEK_SET);
 
    // Optionally apply data filter
    if (datafilter != NULL) {
        string = data_filter(string, datafilter);
    }
 
    // Return string
    return string;
}
 
 
/***********************************************************************//**
 * @brief Read data as integer from HDF5 file
 *
 * @param[in] fptr File pointer.
 * @param[in] nbytes Number of Bytes to read (1...4).
 * @returns Integer.
 *
 * @exception GException::out_of_range
 *            Number of Bytes not in valid range
 *            GException::runtime_error
 *            Attempt to read beyond end of file
 ***************************************************************************/
int gammalib::hdf5::fread_int(FILE* fptr, const int& nbytes)
{
    // Initialise buffer for reading
    char buffer[4];
 
    // Check that nbytes is in valid range
    if (nbytes < 1 || nbytes > 4) {
        throw GException::out_of_range(G_FREAD_INT, "Number of bytes", nbytes, 4);
    }
 
    // Read Bytes from file
    const std::size_t n = std::fread(&buffer[0], 1, nbytes, fptr);
 
    // Check number of Bytes read
    if (n != nbytes) {
        std::string msg = "Attempt to read beyond the end of file; read "
                          "only "+gammalib::str(n)+" of "+gammalib::str(nbytes)+
                          " Bytes from file.";
        throw GException::runtime_error(G_FREAD_INT, msg);
    }
 
    // Initialise integer
    int integer = 0;
 
    // Convert buffer into integer
    for (int i = 0, shift = 0; i < nbytes; ++i, shift += 8) {
        integer += (unsigned char)buffer[i] << shift;
    }
 
    // Return integer
    return (integer);
}
 
 
/***********************************************************************//**
 * @brief Read data as unsiged 32 Bit integer from HDF5 file
 *
 * @param[in] fptr File pointer.
 * @param[in] nbytes Number of Bytes to read (1...4).
 * @returns Unsiged 32 Bit integer.
 *
 * @exception GException::out_of_range
 *            Number of Bytes not in valid range
 *            GException::runtime_error
 *            Attempt to read beyond end of file
 ***************************************************************************/
uint32_t gammalib::hdf5::fread_uint32(FILE* fptr, const int& nbytes)
{
    // Initialise buffer for reading
    char buffer[4];
 
    // Check that nbytes is in valid range
    if (nbytes < 1 || nbytes > 4) {
        throw GException::out_of_range(G_FREAD_UINT32, "Number of bytes", nbytes, 8);
    }
 
    // Read Bytes from file
    const std::size_t n = std::fread(&buffer[0], 1, nbytes, fptr);
 
    // Check number of Bytes read
    if (n != nbytes) {
        std::string msg = "Attempt to read beyond the end of file; read "
                          "only "+gammalib::str(n)+" of "+gammalib::str(nbytes)+
                          " Bytes from file.";
        throw GException::runtime_error(G_FREAD_UINT32, msg);
    }
 
    // Initialise unsigned 32 Bit integer
    uint32_t integer = 0;
 
    // Convert buffer into unsigned 32 Bit integer
    for (int i = 0, shift = 0; i < nbytes; ++i, shift += 8) {
        integer += (unsigned char)buffer[i] << shift;
    }
 
    // Return unsigned 32 Bit integer
    return (integer);
}
 
 
/***********************************************************************//**
 * @brief Read data as unsiged 64 Bit integer from HDF5 file
 *
 * @param[in] fptr File pointer.
 * @param[in] nbytes Number of Bytes to read (1...8).
 * @returns Unsiged 64 Bit integer.
 *
 * @exception GException::out_of_range
 *            Number of Bytes not in valid range
 *            GException::runtime_error
 *            Attempt to read beyond end of file
 ***************************************************************************/
uint64_t gammalib::hdf5::fread_uint64(FILE* fptr, const int& nbytes)
{
    // Initialise buffer for reading
    char buffer[8];
 
    // Check that nbytes is in valid range
    if (nbytes < 1 || nbytes > 8) {
        throw GException::out_of_range(G_FREAD_UINT64, "Number of bytes", nbytes, 8);
    }
 
    // Read Bytes from file
    const std::size_t n = std::fread(&buffer[0], 1, nbytes, fptr);
 
    // Check number of Bytes read
    if (n != nbytes) {
        std::string msg = "Attempt to read beyond the end of file; read "
                          "only "+gammalib::str(n)+" of "+gammalib::str(nbytes)+
                          " Bytes from file.";
        throw GException::runtime_error(G_FREAD_UINT64, msg);
    }
 
    // Initialise unsigned 64 Bit integer
    uint64_t integer = 0;
 
    // Convert buffer into unsigned 64 Bit integer
    for (int i = 0, shift = 0; i < nbytes; ++i, shift += 8) {
        integer += (unsigned char)buffer[i] << shift;
    }
 
    // Return unsigned 64 Bit integer
    return (integer);
}
 
 
/***********************************************************************//**
 * @brief Read string data from HDF5 file
 *
 * @param[in] fptr File pointer.
 * @param[in] nbytes Number of Bytes to read (>0).
 * @returns String.
 *
 * @exception GException::invalid_argument
 *            Number of Bytes not positive
 * @exception GException::runtime_error
 *            Attempt to read beyond end of file
 *
 * Reads a string from the current position of the file pointer @p fptr into
 * a string. The method always reads @p nbytes from the file and advances
 * the file pointer by that amount. If a NULL termination character is
 * encountered the string is truncated to the segement before that character.
 ***************************************************************************/
std::string gammalib::hdf5::fread_string(FILE* fptr, const int& nbytes)
{
    // Check that nbytes is not zero
    if (nbytes < 1) {
        std::string msg = "Non positive number "+gammalib::str(nbytes)+
                          " of Bytes specified. Please specify a number "
                          "larger than 0.";
        throw GException::invalid_argument(G_FREAD_STRING, msg);
    }
 
    // Initialise string
    std::string string;
 
    // Allocate string space
    string.resize(nbytes);
 
    // Read Bytes from file
    const std::size_t n = std::fread(&string[0], 1, nbytes, fptr);
 
    // Check number of Bytes read
    if (n != nbytes) {
        std::string msg = "Attempt to read beyond the end of file; read "
                          "only "+gammalib::str(n)+" of "+gammalib::str(nbytes)+
                          " Bytes from file.";
        throw GException::runtime_error(G_FREAD_STRING, msg);
    }
 
    // Truncate string at termination character
    for (int i = 0; i < string.length(); ++i) {
        if (string[i] == '\0') {             // found
            if (i > 0) {                     // string is not empty
                string = string.substr(0,i); // truncate string
            }
            else {
                string.clear();              // clear string
            }
            break;                           // exit loop
        }
    } // endfor: looped over bytes
 
    // Return string
    return (string);
}
 
 
/***********************************************************************//**
 * @brief Read data with any data type as a string
 *
 * @param[in] fptr File pointer.
 * @param[in] datatype Datatype.
 * @return Data converted into a string.
 *
 * @exception GException::runtime_error
 *            Attempt reading beyond end of file
 *
 * Reads one element of any datatype and return it as a string.
 ***************************************************************************/
std::string gammalib::hdf5::fread_data_as_string(FILE* fptr, const GXmlElement* datatype)
{
    // Initialise string
    std::string string;
 
    // Get datatype size and class
    int size = gammalib::toint(datatype->attribute("size"));
    int cls  = gammalib::toint(datatype->attribute("class"));
 
    // Allocate string space
    string.resize(size);
 
    // Read Bytes from file
    const std::size_t n = std::fread(&string[0], 1, size, fptr);
 
    // Check number of Bytes read
    if (n != size) {
        std::string msg = "Attempt to read beyond the end of file; read "
                          "only "+gammalib::str(n)+" of "+gammalib::str(size)+
                          " Bytes from file.";
        throw GException::runtime_error(G_FREAD_DATA_AS_STRING, msg);
    }
 
    // Convert buffer into string
    string = data_to_string(string.c_str(), datatype);
 
    // Return string
    return string;
}
 
 
/***********************************************************************//**
 * @brief Read zero data from HDF5 file
 *
 * @param[in] fptr File pointer.
 * @param[in] nbytes Number of Bytes to read (>0).
 *
 * @exception GException::out_of_range
 *            Number of Bytes not positive
 *            GException::runtime_error
 *            Attempt to read beyond end of file
 ***************************************************************************/
void gammalib::hdf5::fread_zero(FILE* fptr, const int& nbytes)
{
    // Check that nbytes is not zero
    if (nbytes < 1) {
        std::string msg = "Non positive number "+gammalib::str(nbytes)+
                          " of Bytes specified. Please specify a number "
                          "larger than 0.";
        throw GException::invalid_argument(G_FREAD_ZERO, msg);
    }
 
    // Initialise buffer for reading
    char* buffer = new char(nbytes);
 
    // Get current file pointer position
    uint64_t pos = std::ftell(fptr);
 
    // Read Bytes from file
    const std::size_t n = std::fread(&buffer[0], 1, nbytes, fptr);
 
    // Check number of Bytes read
    if (n != nbytes) {
        std::string msg = "Attempt to read beyond the end of file; read "
                          "only "+gammalib::str(n)+" of "+gammalib::str(nbytes)+
                          " Bytes from file.";
        throw GException::runtime_error(G_FREAD_ZERO, msg);
    }
 
    // Check that all Bytes are zero
    for (int i = 0; i < nbytes; ++i) {
        if ((unsigned char)buffer[i] != 0) {
            std::string msg = "Non-zero Byte "+gammalib::str((int)buffer[i])+" encountered "
                              "at position "+gammalib::str(pos+i)+" in file.";
            throw GException::runtime_error(G_FREAD_ZERO, msg);
        }
    }
 
    // Free buffer
    delete buffer;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Convert data into integer value
 *
 * @param[in] data Pointer to data.
 * @param[in] datatype Pointer to datatype XML element.
 * @return Data converted into integer value.
 *
 * @exception GException::feature_not_implemented
 *            Conversion of specified datatype not yet supported
 *            GException::invalid_value
 *            Conversion of datatype not supported
 *
 * Converts one @p data element of given @p datatype into an integer value.
 * The method does not check for overflow. Fixed-point data are assumed to
 * be stored in little endian order.
 ***************************************************************************/
int gammalib::hdf5::data_to_int(const char*        data,
                                const GXmlElement* datatype)
{
    // Initialise integer value
    int value = 0;
 
    // Get datatype class
    int size = gammalib::toint(datatype->attribute("size"));
    int cls  = gammalib::toint(datatype->attribute("class"));
 
    // Convert data according to class
    switch (cls) {
 
    // Fixed point
    case 0: {
        // Get order
        int order = gammalib::toint(datatype->attribute("order"));
 
        // Handle little-endian
        if (order == 0) {
            for (int i = 0, shift = 0; i < size; ++i, shift += 8) {
                value += (int)((unsigned char)data[i] << shift);
            }
        }
 
        // Handle big-endian
        else {
            for (int i = 0, shift = 0; i < size; ++i, shift += 8) {
                value += (int)((unsigned char)data[size-i-1] << shift);
            }
        }
        break;
    }
 
    // Floating point
    case 1: {
        if (size == 4) {
            value = (int)(*((const float*)data));
        }
        else if (size == 8) {
            value = (int)(*((const double*)data));
        }
        else {
            std::string msg = "Conversion of floating point value of "
                              "size "+gammalib::str(size)+" to integer "
                              "value is not supported.";
            throw GException::invalid_value(G_DATA_TO_INT, msg);
        }
        break;
    }
 
    // Other cases
    default: {
        std::string msg = "Conversion of \""+classname(cls)+"\" to integer "
                          "value is not yet implemented.";
        throw GException::feature_not_implemented(G_DATA_TO_INT, msg);
        break;
    }
 
    } // endswitch: handle type dependent data
 
    // Return value
    return value;
}
 
 
/***********************************************************************//**
 * @brief Convert data into double precision floating point value
 *
 * @param[in] data Pointer to data.
 * @param[in] datatype Pointer to datatype XML element.
 * @return Data converted into double precision floating point value.
 *
 * @exception GException::feature_not_implemented
 *            Conversion of specified datatype not yet supported
 *            GException::invalid_value
 *            Conversion of datatype not supported
 *
 * Converts one @p data element of given @p datatype into a double precision
 * floating point value.
 ***************************************************************************/
double gammalib::hdf5::data_to_double(const char*        data,
                                      const GXmlElement* datatype)
{
    // Initialise double precision floating point value
    double value = 0.0;
 
    // Get datatype class
    int size = gammalib::toint(datatype->attribute("size"));
    int cls  = gammalib::toint(datatype->attribute("class"));
 
    // Convert data according to class
    switch (cls) {
 
    // Fixed point
    case 0: {
        // Get order
        int order = gammalib::toint(datatype->attribute("order"));
 
        // Handle little-endian
        if (order == 0) {
            for (int i = 0, shift = 0; i < size; ++i, shift += 8) {
                value += (double)((unsigned char)data[i] << shift);
            }
        }
 
        // Handle big-endian
        else {
            for (int i = 0, shift = 0; i < size; ++i, shift += 8) {
                value += (double)((unsigned char)data[size-i-1] << shift);
            }
        }
        break;
    }
 
    // Floating point
    case 1: {
        if (size == 4) {
            value = (double)(*((const float*)data));
        }
        else if (size == 8) {
            value = *((const double*)data);
        }
        else {
            std::string msg = "Conversion of floating point value of "
                              "size "+gammalib::str(size)+" to double "
                              "precision value is not supported.";
            throw GException::invalid_value(G_DATA_TO_DOUBLE, msg);
        }
        break;
    }
 
    // Other cases
    default: {
        std::string msg = "Conversion of \""+classname(cls)+"\" to double "
                          "precision value is not yet implemented.";
        throw GException::feature_not_implemented(G_DATA_TO_DOUBLE, msg);
        break;
    }
 
    } // endswitch: handle type dependent data
 
    // Return value
    return value;
}
 
 
/***********************************************************************//**
 * @brief Convert data into string
 *
 * @param[in] data Pointer to data.
 * @param[in] datatype Pointer to datatype XML element.
 * @return Data converted into string.
 *
 * @exception GException::feature_not_implemented
 *            Conversion of specified datatype not yet supported
 *
 * Converts one @p data element of given @p datatype into a character string.
 ***************************************************************************/
std::string gammalib::hdf5::data_to_string(const char*        data,
                                           const GXmlElement* datatype)
{
    // Initialise string
    std::string string;
 
    // Get datatype class
    int size = gammalib::toint(datatype->attribute("size"));
    int cls  = gammalib::toint(datatype->attribute("class"));
 
    // Convert data according to class
    switch (cls) {
 
    // Fixed point
    case 0: {
 
        // Get class specific attributes
        int sign = gammalib::toint(datatype->attribute("signed"));
 
        // Handle signed case
        if (sign == 1) {
 
            // Initialise signed 64 Bit integer
            int64_t value = 0;
 
            // Convert buffer into signed 64 Bit integer (assumes little-endian)
            for (int i = 0, shift = 0; i < size; ++i, shift += 8) {
                value += (unsigned char)data[i] << shift;
            }
 
            // Set string
            string = gammalib::str(value);
 
        } // endif: fixed point was signed
 
        // Handle unsigned case
        else {
 
            // Initialise unsigned 64 Bit integer
            uint64_t value = 0;
 
            // Convert buffer into unsigned 64 Bit integer (assumes little-endian)
            for (int i = 0, shift = 0; i < size; ++i, shift += 8) {
                value += (unsigned char)data[i] << shift;
            }
 
            // Set string
            string = gammalib::str(value);
 
        } // endelse: fixed point was unsigned
        break;
    }
 
    // Floating point
    /*
    case 1: {
        break;
    }
    */
 
    // Enumerated
    case 8: {
 
        // Get class specific attributes
        int number = gammalib::toint(datatype->attribute("number"));
 
        // Get pointer to basetype
        const GXmlElement* basetype = datatype->element("basetype");
 
        // Convert buffer into string
        string = data_to_string(data, basetype);
        break;
    }
 
    // Other cases
    default: {
        std::string msg = "Conversion of \""+classname(cls)+"\" to string "
                          "not yet implemented.";
        throw GException::feature_not_implemented(G_DATA_TO_STRING, msg);
        break;
    }
 
    } // endswitch: handle type dependent data
 
    // Return string
    return string;
}
 
 
/***********************************************************************//**
 * @brief Filter data
 *
 * @param[in] data Data string.
 * @param[in] datafilter Pointer to datafilter XML element.
 * @return Filtered data string.
 *
 * @exception GException::invalid_argument
 *            No "filters" attribute in @p datafilter
 * @exception GException::invalid_value
 *            Unxpected number of "filter" elements
 * @exception GException::feature_not_implemented
 *            Requested filter is not yet implemented
 *
 * Filter data according to specified @p datafilter.
 ***************************************************************************/
std::string gammalib::hdf5::data_filter(const std::string& data,
                                        const GXmlElement* datafilter)
{
    // Initialise string
    std::string string;
 
    // Check for required attributes
    if (!datafilter->has_attribute("filters")) {
        std::string msg = "\"DataStorageFilterPipeline\" message is missing "
                          "\"filters\" attribute";
        throw GException::invalid_argument(G_DATA_FILTER, msg);
    }
 
    // Get number of filters
    int filters = gammalib::toint(datafilter->attribute("filters"));
 
    // Check for corresponding number of filter elements
    if (datafilter->elements("filter") != filters) {
        std::string msg = gammalib::str(datafilter->elements("filter"))+
                          "filters found in \"DataStorageFilterPipeline\" "
                          "but attribute \"filters\"="+gammalib::str(filters)+
                          ". Please specify a valid HDF5 file.";
        throw GException::invalid_value(G_DATA_FILTER, msg);
    }
 
    // Loop over number of filters
    for (int i = 0; i < filters; ++i) {
 
        // Get filter XML element
        const GXmlElement* filter = datafilter->element("filter");
 
        // Check for required attributes
        if (!filter->has_attribute("id")) {
            std::string msg = "Filter element "+gammalib::str(i)+" is "
                              "missing \"id\" attribute";
            throw GException::invalid_argument(G_DATA_FILTER, msg);
        }
 
        // Get identifier
        int id = gammalib::toint(filter->attribute("id"));
 
        // Apply filter according to ID
        switch (id) {
        // Bitshuffle
        case 32008: {
            string = data_filter_bitshuffle(data, filter);
            break;
        }
        default: {
            std::string msg = "Requested data filter with id="+
                              gammalib::str(id)+" is not yet implemented.";
            throw GException::feature_not_implemented(G_DATA_FILTER, msg);
            break;
        }
        }
 
    } // endfor: looped over data filters
 
    // Return string
    return string;
}
 
 
/***********************************************************************//**
 * @brief Filter data according to Bitshuffle algorithm
 *
 * @param[in] data Data string.
 * @param[in] filter Pointer to bitshuffle filter XML element.
 * @return Filtered data string.
 *
 * @exception GException::invalid_argument
 *            Invalid number of client data values found in @p filter
 * @exception GException::feature_not_implemented
 *            Unsupported Bitshuffle algorithm specified in @p filter
 *
 * Filter data according to Bitshuffle algorithm with parameters specified
 * by @p filter.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the bshuf_h5_filter() function in the file bshuf_h5filter.c.
 ***************************************************************************/
std::string gammalib::hdf5::data_filter_bitshuffle(const std::string& data,
                                                   const GXmlElement* filter)
{
    // Initialise string
    std::string string;
 
    // Debug bitshuffle code
    #if defined(G_BITSHUFFLE_DEBUG)
    std::cout << "data_filter_bitshuffle" << std::endl;
    std::cout << "======================" << std::endl;
    #endif
 
    // Get client data values
    int client_values = filter->elements("data");
 
    // Check client data values
    if (client_values < 3) {
        std::string msg = "Only "+gammalib::str(client_values)+
                          " client data values found in filter, but at least "
                          "3 values are required for the bitshuffle filter. "
                          "Please specify a valid HDF5 file.";
        throw GException::invalid_argument(G_DATA_FILTER_BITSHUFFLE, msg);
    }
 
    // Extract element size of typed data
    size_t element_size = filter->element("data", 2)->integer();
 
    // Optionally extract block size
    size_t block_size = (client_values > 3) ? filter->element("data", 3)->integer() : 0;
 
    // If block size is zero then use default block size
    if (block_size == 0) {
        block_size = 8192 / element_size;                   // Target a block size of 8192
        block_size = (block_size / 8) * 8;                  // Block size must be a multiple of 8
        block_size = (block_size > 128) ? block_size : 128; // Block size should be at least 128
    }
 
    // Optionally extract compression algorithm
    int compression = (client_values > 4) ? filter->element("data", 4)->integer() : 0;
 
    // Check compression algorithm
    if (compression != 2) {
        std::string msg = "Unsupported bitshuffle compression algorithm "+
                          gammalib::str(compression)+" encountered. "
                          "Only algorithm=2 (LZ4 compression) is so far supported.";
        throw GException::feature_not_implemented(G_DATA_FILTER_BITSHUFFLE, msg);
    }
 
    // Set pointer to input buffer
    char* input = (char*)data.c_str();
 
    // Handle dependent on compression type
    switch(compression) {
 
    // No compression
    case 0: {
        //
        size_t nbytes_uncompressed = data.length();
        size_t buf_size_out        = data.length();
        size_t size                = nbytes_uncompressed / element_size;
 
        // Allocate string space for uncompressed data
        string.resize(buf_size_out);
 
        // Set pointer to output buffer
        char* output = (char*)string.c_str();
 
        // Bit unshuffle.
        //err = bshuf_bitunshuffle(input, output, size, element_size, block_size);
        break;
    }
 
    // LZ4 compression
    case 2: {
        // Read number of uncompressed Bytes and block size from data
        size_t nbytes_uncompressed = bitshuffle_read_uint64(input);
        size_t block_size          = bitshuffle_read_uint32((const char*)input + 8) / element_size;
        input                     += 12;
 
        // Compute number of elements
        size_t elements = nbytes_uncompressed / element_size;
 
        // Allocate string space for uncompressed data
        string.resize(nbytes_uncompressed);
 
        // Set pointer to output buffer
        char* output = (char*)string.c_str();
 
        // Debug bitshuffle code
        #if defined(G_BITSHUFFLE_DEBUG)
        std::cout << "LZ4 compression:";
        std::cout << " nbytes_uncompressed=" << nbytes_uncompressed;
        std::cout << ", block_size=" << block_size;
        std::cout << ", elements=" << elements;
        std::cout << ", element_size=" << element_size << std::endl;
        #endif
 
        // Decompress and unshuffle
        bitshuffle_decompress(input, output, elements, element_size, block_size);
        break;
    }
    
    } // endswitch: handled according to compression type
 
    // Return string
    return string;
}
 
 
/***********************************************************************//**
 * @brief Return specific occurence of header message of given type
 *
 * @param[in] header Pointer to XML header element.
 * @param[in] type Message type.
 * @param[in] number Occurence of message.
 * @return Pointer to XML element of specific occurence of message type.
 *
 * GException::invalid_value
 *             No @p type message found in header.
 *
 * Returns a XML pointer to the ocurrence @p number of a message of specified
 * @p type.
 ***************************************************************************/
const GXmlElement* gammalib::hdf5::xml_msg_type(const GXmlElement* header,
                                                const std::string& type,
                                                const int&         number)
{
    // Initialise result element
    const GXmlElement* element = NULL;
 
    // Get number of messages
    int messages = header->elements("message");
 
    // Initialise occurences
    int occurences = 0;
 
    // Loop over messages
    for (int i = 0; i < messages; ++i) {
 
        // Get message element
        const GXmlElement* message = header->element("message", i);
 
        // Search the element for requested type
        if (message->has_attribute("type")) {
            if (message->attribute("type") == type) {
                if (occurences == number) {
                    element = message;
                    break;
                }
                else {
                    occurences++;
                }
            }
        }
 
    } // endfor: looped over messages
 
    // Throw an exception if no entry was found
    if (element == NULL) {
        if (number == 0) {
            std::string msg = "No \""+type+"\" message found in HDF5 object header.";
            throw GException::invalid_value(G_XML_MSG_TYPE, msg);
        }
        else {
            std::string msg = gammalib::order(number)+" occurence of \""+
                              type+"\" message not found in HDF5 object header.";
            throw GException::invalid_value(G_XML_MSG_TYPE, msg);
        }
    }
 
    // Return element
    return element;
}
 
 
/***********************************************************************//**
 * @brief Checks if specific occurence of header message of given type exists
 *
 * @param[in] header Pointer to XML header element.
 * @param[in] type Message type.
 * @param[in] number Occurence of message.
 * @return True if ith occurence of message exists, false otherwise.
 *
 * Returns true if the occurence @p number of header message of @p type
 * exists.
 ***************************************************************************/
bool gammalib::hdf5::xml_has_msg_type(const GXmlElement* header,
                                      const std::string& type,
                                      const int&         number)
{
    // Initialise result element
    const GXmlElement* element = NULL;
 
    // Get number of messages
    int messages = header->elements("message");
 
    // Initialise occurences
    int occurences = 0;
 
    // Loop over messages
    for (int i = 0; i < messages; ++i) {
 
        // Get message element
        const GXmlElement* message = header->element("message", i);
 
        // Search the element for with specified type
        if (message->has_attribute("type")) {
            if (message->attribute("type") == type) {
                if (occurences == number) {
                    element = message;
                    break;
                }
                else {
                    occurences++;
                }
            }
        }
 
    } // endfor: looped over messages
 
    // Return existence
    return (element != NULL);
}
 
/***********************************************************************//**
 * @brief Return number of "Attribute" messages in header
 *
 * @param[in] header Pointer to XML header element.
 * @param[in] name Name of "Attribute" message.
 * @return Number of attributes in message header.
 *
 * Returns the number of message elements of type "Attribute" in the specified
 * header element. If @p name is not empty, the method will only count the
 * attributes of which the "name" attribute equals to @p name.
 ***************************************************************************/
int gammalib::hdf5::xml_msg_attributes(const GXmlElement* header,
                                       const std::string& name)
{
    // Initialise number of attributes
    int number = 0;
 
    // Get number of messages
    int messages = header->elements("message");
 
    // Signal that name should be checked
    bool check_name = (name.length() > 0);
 
    // Loop over messages
    for (int i = 0; i < messages; ++i) {
 
        // Get message element
        const GXmlElement* message = header->element("message", i);
 
        // Count number of "Attribute" types
        if (message->has_attribute("type")) {
            if (message->attribute("type") == "Attribute") {
                if (check_name) {
                    if (message->has_attribute("name")) {
                        if (message->attribute("name") == name) {
                            number++;
                        }
                    }
                }
                else {
                 number++;
                }
            }
        }
 
    } // endfor: looped over messages
 
    // Return number of attributes
    return number;
}
 
 
/***********************************************************************//**
 * @brief Return "Attribute" message XML element with specified name
 *
 * @param[in] header Pointer to XML header element.
 * @param[in] name Name of "Attribute" message.
 * @param[in] number Number of occurences to return (0=return first occurence).
 * @return Pointer to "Attribute" message XML element with specified @p name.
 *
 * Returns a XML pointer to a message element of type "Attribute" and with
 * the specified @p name.
 *
 * If several "Attribute" messages have the same @p name, the first of these
 * messages is returned.
 ***************************************************************************/
const GXmlElement* gammalib::hdf5::xml_msg_attribute(const GXmlElement* header,
                                                     const std::string& name,
                                                     const int&         number)
{
    // Initialise result element
    const GXmlElement* element = NULL;
 
    // Get number of messages
    int messages = header->elements("message");
 
    // Initialise occurences
    int occurences = 0;
 
    // Loop over messages
    for (int i = 0; i < messages; ++i) {
 
        // Get message element
        const GXmlElement* message = header->element("message", i);
 
        // Search the element for "Attribute" with specified name
        if (message->has_attribute("type")) {
            if (message->attribute("type") == "Attribute") {
                if (message->has_attribute("name")) {
                    if (message->attribute("name") == name) {
                        if (occurences == number) {
                            element = message;
                            break;
                        }
                        else {
                            occurences++;
                        }
                    }
                }
            }
        }
 
    } // endfor: looped over messages
 
    // Throw an exception if no entry was found
    if (element == NULL) {
        if (number == 0) {
            std::string msg = "No \"Attribute\" message with attribute name=\""+
                              name+"\" found in HDF5 metadata.";
            throw GException::invalid_value(G_XML_MSG_ATTRIBUTE, msg);
        }
        else {
            std::string msg = gammalib::order(number)+" occurence of "
                              "\"Attribute\" message with attribute name=\""+
                              name+"\" not found in HDF5 metadata.";
            throw GException::invalid_value(G_XML_MSG_ATTRIBUTE, msg);
        }
    }
 
    // Return element
    return element;
}
 
 
/***********************************************************************//**
 * @brief Checks if "Attribute" message XML element with specified name exists
 *
 * @param[in] header Pointer to XML header element.
 * @param[in] name Name of "Attribute" message.
 * @param[in] number Number of occurences to return (0=return first occurence).
 * @return True if message exists, false otherwise.
 *
 * Returns true if an "Attribute" message with @p name exists.
 ***************************************************************************/
bool gammalib::hdf5::xml_has_msg_attribute(const GXmlElement* header,
                                           const std::string& name,
                                           const int&         number)
{
    // Initialise result element
    const GXmlElement* element = NULL;
 
    // Get number of messages
    int messages = header->elements("message");
 
    // Initialise occurences
    int occurences = 0;
 
    // Loop over messages
    for (int i = 0; i < messages; ++i) {
 
        // Get message element
        const GXmlElement* message = header->element("message", i);
 
        // Search the element for "Attribute" with specified name
        if (message->has_attribute("type")) {
            if (message->attribute("type") == "Attribute") {
                if (message->has_attribute("name")) {
                    if (message->attribute("name") == name) {
                        if (occurences == number) {
                            element = message;
                            break;
                        }
                        else {
                            occurences++;
                        }
                    }
                }
            }
        }
 
    } // endfor: looped over messages
 
    // Return existence
    return (element != NULL);
}
 
 
/***********************************************************************//**
 * @brief Return "Attribute" message key value from attribute with specified name
 *
 * @param[in] header Pointer to XML header element.
 * @param[in] name Name of "Attribute" message.
 * @param[in] key Key of "Attribute" message.
 * @return String with value of @p key in "Attribute" message with @p name.
 *
 * Returns the value of an attribute @p key in an "Attribute" message with
 * the specified @p name as a string. If the attribute or the key was not
 * found then return an empty string.
 *
 * If several "Attribute" messages have the same @p name and p@ key, the
 * first corresponding value is returned.
 ***************************************************************************/
std::string gammalib::hdf5::xml_msg_attribute(const GXmlElement* header,
                                              const std::string& name,
                                              const std::string& key)
{
    // Initialise value string
    std::string value;
 
    // Get number of messages
    int messages = header->elements("message");
 
    // Loop over messages
    for (int i = 0; i < messages; ++i) {
 
        // Get message element
        const GXmlElement* message = header->element("message", i);
 
        // Search the element for "Attribute" with specified name and
        // specified key and return it's value
        if (message->has_attribute("type")) {
            if (message->attribute("type") == "Attribute") {
                if (message->has_attribute("name")) {
                    if (message->attribute("name") == name) {
                        if (message->has_attribute(key)) {
                            value = message->attribute(key);
                            break;
                        }
                    }
                }
            }
        }
 
    } // endfor: looped over all messages
 
    // Return value
    return value;
}
 
 
/***********************************************************************//**
 * @brief Convert datatype class into classname
 *
 * @param[in] cls Datatype class.
 * @return Classname.
 *
 * Return class name from datatype class.
 ***************************************************************************/
std::string gammalib::hdf5::classname(const int& cls)
{
    // Initialise classname
    std::string classname;
 
    // Convert datatype class into string
    switch (cls) {
    case 0:
        classname = "Fixed point";
        break;
    case 1:
        classname = "Floating point";
        break;
    case 2:
        classname = "Time";
        break;
    case 3:
        classname = "String";
        break;
    case 4:
        classname = "Bit field";
        break;
    case 5:
        classname = "Opaque";
        break;
    case 6:
        classname = "Compound";
        break;
    case 7:
        classname = "Reference";
        break;
    case 8:
        classname = "Enumerated";
        break;
    case 9:
        classname = "Variable-Length";
        break;
    case 10:
        classname = "Array";
        break;
    case 11:
        classname = "Complex";
        break;
    default:
        break;
    }
 
    // Return classname
    return classname;
}
 
 
/*==========================================================================
 =                                                                         =
 =                       Private HDF5 support functions                    =
 =                                                                         =
 ==========================================================================*/
 
 
/***********************************************************************//**
 * @brief Decompress and bitshuffle data
 *
 * @param[in] input Pointer to input buffer.
 * @param[out] out Pointer to output buffer (@p elements * @p element_size bytes).
 * @param[in] elements Number of elements in input.
 * @param[in] element_size Element size of typed data.
 * @param[in] block_size Process in blocks of this many elements.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the bshuf_decompress_lz4() and bshuf_decompress_lz4_block()
 * functions in the bitshuffle.c file and the bshuf_blocked_wrap_fun()
 * function in the bitshuffle_core.c file.
 ***************************************************************************/
void bitshuffle_decompress(char*         input,
                           char*         output,
                           const size_t& elements,
                           const size_t& element_size,
                           size_t&       block_size)
{
    // Initialise cumulative count
    uint64_t cum_count = 0;
 
    // Determine number of blocks
    int nblocks = elements / block_size;
 
    // Debug bitshuffle code
    #if defined(G_BITSHUFFLE_DEBUG)
    std::cout << "bitshuffle_decompress" << std::endl;
    std::cout << "=====================" << std::endl;
    std::cout << "nblocks=" << nblocks << std::endl;
    #endif
 
    // Loop over blocks, including the last block
    for (int i = 0; i < nblocks+1; ++i) {
 
        // Determine size of this block, where each block except of the
        // last block have the same user specified size
        size_t this_block_size = block_size;
        if (i == nblocks) {
            this_block_size = elements % block_size;
            this_block_size = this_block_size - this_block_size % 8;
        }
 
        // Continue only if size of this block is non-zero
        if (this_block_size != 0) {
 
            // Read number of Bytes to decompress from header
            uint32_t nbytes_from_header = bitshuffle_read_uint32(input);
 
            // Set buffer size
            int buffer_size = this_block_size * element_size;
 
            // Allocate temporary memory
            char* buffer = new char[buffer_size];
 
            // LZ4 decompression
            int nbytes = decompress_lz4((const char*)input + 4,
                                        buffer,
                                        nbytes_from_header,
                                        buffer_size);
 
            // Debug bitshuffle code
            #if defined(G_BITSHUFFLE_DEBUG)
            std::cout << "block=" << i;
            std::cout << ", nbytes_from_header=" << nbytes_from_header;
            std::cout << ", block_size=" << this_block_size;
            std::cout << ", buffer_size=" << buffer_size;
            std::cout << ", nbytes=" << nbytes;
            std::cout << std::endl;
            #endif
 
            // Check number of decompressed Bytes
            if (nbytes != buffer_size) {
                std::string msg = "Invalid number of "+gammalib::str(nbytes)+
                                  " decompressed Bytes, expected "+
                                  gammalib::str(buffer_size)+" Bytes.";
                throw GException::runtime_error(G_FREAD_DATA1, msg);
            }
 
            // Unshuffle elements
            bitshuffle_elements(buffer, output, this_block_size, element_size);
 
            // Increment input and output pointers
            input  += nbytes_from_header + 4;
            output += this_block_size * element_size;
 
            // Cumulate counts
            cum_count += nbytes_from_header + 4;
 
            // Free buffer
            delete [] buffer;
 
        } // endif: size of this block was not zero
 
    } // endfor: looped over blocks
 
    // Compute number of leftover Bytes
    size_t leftover_bytes = elements % 8 * element_size;
 
    // Debug bitshuffle code
    #if defined(G_BITSHUFFLE_DEBUG)
    std::cout << "leftover_bytes=" << leftover_bytes;
    std::cout << std::endl;
    #endif
 
    // Copy over leftover Bytes
    //memcpy(output, input, leftover_bytes);
    for (int i = 0; i < leftover_bytes; ++i) {
        *output++ = *input++;
    }
 
    // Cumulate counts
    cum_count += leftover_bytes;
 
    // Debug bitshuffle code
    #if defined(G_BITSHUFFLE_DEBUG)
    std::cout << "cum_count=" << cum_count;
    std::cout << std::endl;
    #endif
 
    // Return
    return;
 
}
 
 
/***********************************************************************//**
 * @brief Read a 64 bit unsigned integer from a buffer big endian order
 *
 * @param[in] buffer Buffer.
 * @return 64 bit unsigned integer.
 *
 * Reads a 64 bit unsigned integer stored in big endian order from a Byte
 * buffer.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the bshuf_read_uint64_BE() function in the bitshuffle.c file.
 ***************************************************************************/
uint64_t bitshuffle_read_uint64(const char* buffer)
{
    // Set constants
    const uint64_t pow28 = 256;
 
    // Cast buffer to unsigned Byte
    uint8_t* b = (uint8_t*)buffer;
 
    // Initialise conversion
    uint64_t num = 0;
    uint64_t cp  = 1;
 
    // Loop over 8 Bytes backwards
    for (int i = 7; i >= 0; i--) {
        num += b[i] * cp;
        cp  *= pow28;
    }
 
    // Returns result
    return num;
}
 
 
/***********************************************************************//**
 * @brief Read a 32 bit unsigned integer from a buffer big endian order
 *
 * @param[in] buffer Buffer.
 * @return 32 bit unsigned integer.
 *
 * Reads a 32 bit unsigned integer stored in big endian order from a Byte
 * buffer.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the bshuf_read_uint32_BE() function in the bitshuffle.c file.
 ***************************************************************************/
uint32_t bitshuffle_read_uint32(const char* buffer)
{
    // Set constants
    const uint32_t pow28 = 256;
 
    // Cast buffer to unsigned Byte
    uint8_t* b = (uint8_t*)buffer;
 
    // Initialise conversion
    uint32_t num = 0;
    uint32_t cp  = 1;
 
    // Loop over 4 Bytes backwards
    for (int i = 3; i >= 0; i--) {
        num += b[i] * cp;
        cp  *= pow28;
    }
 
    // Return result
    return num;
}
 
 
/***********************************************************************//**
 * @brief Untranspose bits within elements
 *
 * @param[in] input Pointer to input buffer.
 * @param[out] output Pointer to output buffer (@p elements * @p element_size bytes).
 * @param[in] elements Number of elements in input.
 * @param[in] element_size Element size of typed data.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the bshuf_untrans_bit_elem_scal() function in the file
 * bitshuffle_core.c.
 ***************************************************************************/
void bitshuffle_elements(const char*   input,
                         char*         output,
                         const size_t& elements,
                         const size_t& element_size)
{
    // Compute number of Bytes on output
    int nbytes = elements * element_size;
 
    // Allocate temporary buffer
    char* buffer = new char[nbytes];
 
    // Compute number of Bytes in row
    int nbyte_row = elements / 8;
 
    // For data organized into a row for each bit (8 * element_size rows),
    // transpose the bytes
    for (int i = 0; i < element_size; ++i) {
        int i8 = i * 8;
        for (int j = 0; j < nbyte_row; ++j) {
            int joffset = j * 8 * element_size + i8;
            for (int k = 0; k < 8; ++k) {
                buffer[joffset + k] = input[(i8 + k) * nbyte_row + j];
            }
        }
    }
 
    // Shuffle bits within the bytes of eight element blocks. The shuffling
    // depends on the storage order. First handle little endian...
    if (gammalib::little_endian()) {
        for (int i = 0; i < 8 * element_size; i += 8) {
            for (int j = 0; j + 8 * element_size - 1 < nbytes; j += 8 * element_size) {
                uint64_t x = *((uint64_t*)&buffer[j + i]);
                uint64_t t = (x ^ (x >> 7)) & 0x00AA00AA00AA00AALL;
                x          = x ^ t ^ (t << 7);
                t          = (x ^ (x >> 14)) & 0x0000CCCC0000CCCCLL;
                x          = x ^ t ^ (t << 14);
                t          = (x ^ (x >> 28)) & 0x00000000F0F0F0F0LL;
                x          = x ^ t ^ (t << 28);
                for (int k = 0; k < 8; ++k) {
                    *((uint8_t*)&output[j + i / 8 + k * element_size]) = x;
                    x                                                  = x >> 8;
                }
            }
        }
    } // endif: integers are stored little endian
 
    // ... otherwise handle big endian
    else {
        for (int i = 0; i < 8 * element_size; i += 8) {
            for (int j = 0; j + 8 * element_size - 1 < nbytes; j += 8 * element_size) {
                uint64_t x = *((uint64_t*)&buffer[j + i]);
                uint64_t t = (x ^ (x >> 9)) & 0x0055005500550055LL;
                x          = x ^ t ^ (t << 9);
                t          = (x ^ (x >> 18)) & 0x0000333300003333LL;
                x          = x ^ t ^ (t << 18);
                t          = (x ^ (x >> 36)) & 0x000000000F0F0F0FLL;
                x          = x ^ t ^ (t << 36);
                for (int k = 0; k < 8; ++k) {
                    *((uint8_t*)&output[j + i / 8 + (7-k) * element_size]) = x;
                    x                                                      = x >> 8;
                }
            }
        }
    } // endelse: integers are stored big endian
 
    // Free temporary buffer
    delete [] buffer;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Decompress data using the LZ4 algorithm
 *
 * @param[in] src Pointer to compressed data.
 * @param[out] dest Pointer to decompressed data (must be already allocated).
 * @param[in] srcSize Size of compressed data.
 * @param[in] outputSize Maximum size of allocated decompressed data.
 * @returns Number of decompressed Bytes.
 *
 * Decompress data using the LZ4 algorithm.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the LZ4_decompress_safe() function in the file lz4.c.
 ***************************************************************************/
int decompress_lz4(const char*   src,
                   char*         dst,
                   const size_t& srcSize,
                   const size_t& outputSize)
{
    // Set constants
    const int      minmatch                 = 4;
    const int      wildcopylength           = 8;
    const int      lastliterals             = 5;
    const int      mflimit                  = 12;
    const int      match_safeguard_distance = ((2*wildcopylength) - minmatch);
    const int      fastloop_safe_distance   = 64;
    const int      ml_bits                  = 4;
    const int      ml_mask                  = ((1U<<ml_bits)-1);
    const int      run_bits                 = (8-ml_bits);
    const int      run_mask                 = ((1U<<run_bits)-1);
    const unsigned inc32table[8]            = {0, 1, 2,  1,  0,  4, 4, 4};
    const int      dec64table[8]            = {0, 0, 0, -1, -4,  1, 2, 3};
 
    // Derive variables from arguments
    const uint8_t* lowPrefix = (uint8_t*)dst;
 
    // Set start and end Byte pointers to compressed data
    const uint8_t* ip   = (const uint8_t*)src;
    const uint8_t* iend = ip + srcSize;
 
    // Set start and end Byte pointers to uncompressed data
    uint8_t* op   = (uint8_t*)dst;
    uint8_t* oend = op + outputSize;
    uint8_t* cpy;
 
    // Set up the "end" pointers for the shortcut
    const uint8_t* shortiend = iend - 14 /*maxLL*/ - 2 /*offset*/;
    const uint8_t* shortoend = oend - 14 /*maxLL*/ - 18 /*maxML*/;
 
    // Allocate working variables
    const uint8_t* match;
    size_t              offset;
    unsigned            token;
    size_t              length;
 
    // Use fast loop
    if ((oend - op) >= fastloop_safe_distance) {
 
        // Decode sequences as long as output < oend-fastloop_safe_distance
        while (true) {
 
            // Main fastloop assertion: We can always wildcopy fastloop_safe_distance
            //assert(oend - op >= fastloop_safe_distance);
            //assert(ip < iend);
 
            // Get next token and literal length
            token  = *ip++;
            length = token >> ml_bits;
 
            // Decode literal length
            if (length == run_mask) {
 
                // Add length
                length += read_variable_length(&ip, iend-run_mask, true);
 
                // Check for overflow
                if (((uintptr_t)(op)+length < (uintptr_t)(op)) ||
                    ((uintptr_t)(ip)+length < (uintptr_t)(ip))) {
                    std::string msg = "Overflow detection.";
                    throw GException::invalid_value(G_DECOMPRESS_LZ4, msg);
                }
 
                // If we're close to the end of the buffers then use safe literal
                // copy
                //LZ4_STATIC_ASSERT(mflimit >= wildcopylength);
                if ((op+length > oend-32) || (ip+length > iend-32)) {
                    goto safe_literal_copy;
                }
 
                // ... otherwise copy literals
                memory_beyond32(ip, op, op+length);
 
                // Increment pointers
                ip += length;
                op += length;
 
            } // endif: decoded literal
 
            // ... otherwise if at least 16 Bytes remain in input buffer then
            // copy 16-Bytes stripe. We don't need to check oend, since we check
            // it once for each loop below.
            else if (ip <= iend-(16 + 1)) {
 
                // Copy 16-Byte stripe
                std::memcpy(op, ip, 16);
 
                // Increment pointers
                ip += length;
                op += length;
 
            } // endif: copied 16-Byte stripe
 
            // ... otherwise use safe literal copy
            else {
                goto safe_literal_copy;
            }
 
            // Get offset
            offset = read_uint16(ip);
            ip    += 2;
            match  = op - offset;
            //assert(match <= op);  /* overflow check */
 
            // Get matchlength
            length = token & ml_mask;
 
            //
            if (length == ml_mask) {
 
                // Add length
                length += read_variable_length(&ip, iend - lastliterals + 1, false);
                length += minmatch;
 
                // Check for overflow
                if ((uintptr_t)(op)+length < (uintptr_t)(op)) {
                    std::string msg = "Overflow detection.";
                    throw GException::invalid_value(G_DECOMPRESS_LZ4, msg);
                }
 
                //
                if (op + length >= oend - fastloop_safe_distance) {
                    goto safe_match_copy;
                }
                
            } // endif
 
            // ...
            else {
                length += minmatch;
                if (op + length >= oend - fastloop_safe_distance) {
                    //DEBUGLOG(7, "moving to safe_match_copy (ml==%u)", (unsigned)length);
                    goto safe_match_copy;
                }
 
                /* Fastpath check: skip memory_beyond32 when true */
                if ((match >= lowPrefix)) {
                    if (offset >= 8) {
                        //assert(match >= lowPrefix);
                        //assert(match <= op);
                        //assert(op + 18 <= oend);
                        std::memcpy(op,    match,    8);
                        std::memcpy(op+8,  match+8,  8);
                        std::memcpy(op+16, match+16, 2);
                        op += length;
                        continue;
                    }
                }
 
            } // endelse
 
            // ...
            if (match < lowPrefix) {
                std::string msg = "Position outside buffers.";
                throw GException::invalid_value(G_DECOMPRESS_LZ4, msg);
            }
 
            // Copy match within block
            cpy = op + length;
 
            //assert((op <= oend) && (oend-op >= 32));
            if (offset < 16) {
                memcpy_using_offset(match, op, cpy, offset);
            }
            else {
                memory_beyond32(match, op, cpy);
            }
 
            // Adopt wildcopy correction
            op = cpy;
 
        } // endwhile
 
    } // endif: used fast loop
 
    // Main loop that decompresses the remaining sequence where
    // output < fastloop_safe_distance
    while (true) {
    
        // Check
        //assert(ip < iend);
 
        // Get literal length token
        token  = *ip++;
        length = token >> ml_bits;
 
        // Two-stage shortcut for the most common case:
        // 1) If the literal length is 0..14, and there is enough space,
        //    enter the shortcut and copy 16 bytes on behalf of the literals
        //    (in the fast mode, only 8 bytes can be safely copied this way).
        // 2) Further if the match length is 4..18, copy 18 bytes in a similar
        //    manner; but we ensure that there's enough space in the output for
        //    those 18 bytes earlier, upon entering the shortcut (in other words,
        //    there is a combined check for both stages).
        if ((length != run_mask) && (ip < shortiend) & (op <= shortoend)) {
 
            // Copy the literals
            std::memcpy(op, ip, 16);
            op += length;
            ip += length;
 
            // Second stage: prepare for match copying, decode full info. If it
            // doesn't work out, the info won't be wasted.
            length = token & ml_mask;
            offset = read_uint16(ip);
            ip    += 2;
            match  = op - offset;
            //assert(match <= op); /* check overflow */
 
            // Do not deal with overlapping matches
            if ((length != ml_mask) && (offset >= 8) && (match >= lowPrefix)) {
                std::memcpy(op + 0, match + 0, 8);
                std::memcpy(op + 8, match + 8, 8);
                std::memcpy(op +16, match +16, 2);
                op += length + minmatch;
                // Both stages worked, load the next token
                continue;
            }
 
            // The second stage didn't work out, but the info is ready. Propel
            // it right to the point of match copying.
            goto _copy_match;
 
        } // endif: Two-stage shortcut
 
        // Decode literal length. Check for overflow.
        if (length == run_mask) {
            length += read_variable_length(&ip, iend-run_mask, true);
            if (((uintptr_t)(op)+length < (uintptr_t)(op)) ||
                ((uintptr_t)(ip)+length < (uintptr_t)(ip))) {
                std::string msg = "Overflow detection.";
                throw GException::invalid_value(G_DECOMPRESS_LZ4, msg);
            }
        }
 
safe_literal_copy:
        /* copy literals */
        cpy = op+length;
 
        //LZ4_STATIC_ASSERT(mflimit >= wildcopylength);
        if ((cpy>oend-mflimit) || (ip+length>iend-(2+1+lastliterals))) {
            /* We've either hit the input parsing restriction or the output parsing restriction.
             * In the normal scenario, decoding a full block, it must be the last sequence,
             * otherwise it's an error (invalid input or dimensions).
             * In partialDecoding scenario, it's necessary to ensure there is no buffer overflow.
             */
            /* We must be on the last sequence (or invalid) because of the parsing limitations
             * so check that we exactly consume the input and don't overrun the output buffer.
             */
            if ((ip+length != iend) || (cpy > oend)) {
                std::string msg = "Overrun output buffer.";
                throw GException::invalid_value(G_DECOMPRESS_LZ4, msg);
            }
 
            std::memmove(op, ip, length);
            ip += length;
            op += length;
            /* Necessarily EOF when !partialDecoding.
             * When partialDecoding, it is EOF if we've either
             * filled the output buffer or
             * can't proceed with reading an offset for following match.
             */
            if ((cpy == oend) || (ip >= (iend-2))) {
                break;
            }
        }
        else {
            memcpy_beyond8(ip, op, cpy);
            ip += length;
            op  = cpy;
        }
 
        // Get offset
        offset = read_uint16(ip);
        ip    += 2;
        match  = op - offset;
 
        // Get matchlength
        length = token & ml_mask;
 
_copy_match:
        if (length == ml_mask) {
 
            // Add variable length
            length += read_variable_length(&ip, iend - lastliterals + 1, false);
 
            // Check for overflow
            if ((uintptr_t)(op)+length < (uintptr_t)(op)) {
                std::string msg = "Overflow detection.";
                throw GException::invalid_value(G_DECOMPRESS_LZ4, msg);
            }
 
        }
        length += minmatch;
 
safe_match_copy:
        // Check if offset is outside buffers
        if (match < lowPrefix) {
            std::string msg = "Offset outside buffers.";
            throw GException::invalid_value(G_DECOMPRESS_LZ4, msg);
        }
 
        // Copy match within block
        cpy = op + length;
 
        //
        if (offset < 8) {
            *(uint32_t*)op = 0;
            op[0]          = match[0];
            op[1]          = match[1];
            op[2]          = match[2];
            op[3]          = match[3];
            match         += inc32table[offset];
            std::memcpy(op+4, match, 4);
            match -= dec64table[offset];
        }
        else {
            std::memcpy(op, match, 8);
            match += 8;
        }
        op += 8;
 
        // ...
        if (cpy > oend-match_safeguard_distance) {
 
            // Get ...
            uint8_t* oCopyLimit = oend - (wildcopylength-1);
 
            // Check ...
            if (cpy > oend-lastliterals) {
                std::string msg = "Last 5 Bytes were not uncompressed literals.";
                throw GException::invalid_value(G_DECOMPRESS_LZ4, msg);
            }
 
            //
            if (op < oCopyLimit) {
                memcpy_beyond8(match, op, oCopyLimit);
                match += oCopyLimit - op;
                op     = oCopyLimit;
            }
            while (op < cpy) {
                *op++ = *match++;
            }
            
        } // endif
 
        else {
            std::memcpy(op, match, 8);
            if (length > 16) {
                memcpy_beyond8(match+8, op+8, cpy);
            }
        } // endelse
 
        // Apply wildcopy correction
        op = cpy;
 
    } // endwhile
 
    // Return number of decompressed Bytes
    return (int)(((char*)op)-dst);
}
 
 
/***********************************************************************//**
 * @brief Read the variable-length literal or match length
 *
 * @param[in] ip Input pointer.
 * @param[in] ilimit Limiting position.
 * @param[in] initial_check Performs initial check.
 * @return Length literal or match length.
 *
 * @exception GException::invalid_argument
 *            Initial check revealed that read limit was reached
 * @exception GException::invalid_value
 *            Read limit was reached or accumulator overflow was detected
 *
 * @p ilimit specifies the position after which, if length is not decoded,
 * the input is necessarily corrupted. In this case an exception is thrown.
 *
 * If @p initial_check is true then a check @p ip < @p ilimit is performed,
 * and if the check is false, an exception is thrown.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the read_variable_length() function in the file lz4.c.
 ***************************************************************************/
size_t read_variable_length(const uint8_t** ip,
                            const uint8_t*  ilimit,
                            const bool&     initial_check)
{
    // Initialise length
    size_t length = 0;
 
    // If initial check is requested then throw an exception if the read
    // limit was reached
    if (initial_check && (*ip) >= ilimit) {
        std::string msg = "Initial check revealed that read limit was reached.";
        throw GException::invalid_argument(G_READ_VARIABLE_LENGTH, msg);
    }
 
    // ... otherwise determine length
    else {
 
        // Allocate size
        size_t s;
 
        // Loop until size is no longer 255
        do {
            // Get size, add it to length and increment input pointer
            s       = **ip;
            length += s;
            (*ip)++;
 
            // Throw an exception if read limit was reached
            if ((*ip) > ilimit) {
                std::string msg = "Read limit was reached.";
                throw GException::invalid_value(G_READ_VARIABLE_LENGTH, msg);
            }
 
            // Stop with error if accumulator overflow detected
            if ((sizeof(length) < 8) && length > ((size_t)(-1)/2)) {
                std::string msg = "Accumulator overflow detected.";
                throw GException::invalid_value(G_READ_VARIABLE_LENGTH, msg);
            }
 
        } while (s == 255);
 
    } // endelse: determined length
 
    // Return
    return length;
}
 
 
/***********************************************************************//**
* @brief Read unsigned 16 Bit integer stored as little endian
*
* @param[in] memPtr Input pointer.
* @return Unsigned 16 Bit integer.
*
* Reads an unsigend 16 Bit integer from the current memory location that is
* stored in little endian order. The code handles both little and big endian
* architectures.
***************************************************************************/
uint16_t read_uint16(const void* memPtr)
{
    // Allocate result
    uint16_t result;
 
    // If architecture is little endian then just store value in memory
    // location
    if (gammalib::little_endian()) {
        result = *(const uint16_t*)(memPtr);
    }
 
    // ... otherwise switch the Byte order
    else {
        const uint8_t* p = (const uint8_t*)memPtr;
        result = (uint16_t)((uint16_t)p[0] | (p[1]<<8));
    }
 
    // Return result
    return result;
}
 
 
/***********************************************************************//**
 * @brief Memory copy
 *
 * @param[in] srcPtr Input pointer.
 * @param[in] dstPtr Pointer to start of output.
 * @param[in] dstEnd Pointer to end of output
 * @param[in] offset .
 *
 * @exception
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the LZ4_memcpy_using_offset() function in the file lz4.c.
 ***************************************************************************/
void memcpy_using_offset(const uint8_t* srcPtr,
                         uint8_t*       dstPtr,
                         uint8_t*       dstEnd,
                         const size_t&  offset)
{
    // Check that destination end pointer is at least 4 Bytes behind start
    // pointer
    if (dstEnd < dstPtr + 4) {
        std::string msg = "Destination end pointer "+
                          gammalib::str((uint64_t)dstPtr)+
                          " is not at least 4 Bytes behind start pointer "+
                          gammalib::str((uint64_t)srcPtr)+".";
        throw GException::invalid_argument(G_MEMCPY_USING_OFFSET, msg);
    }
 
    // Buffer
    uint8_t v[8];
 
    switch (offset) {
    case 1:
        std::memset(v, *srcPtr, 8);
        break;
    case 2:
        std::memcpy(v,     srcPtr, 2);
        std::memcpy(&v[2], srcPtr, 2);
        std::memcpy(&v[4], v,      4);
        break;
    case 4:
        std::memcpy(v,     srcPtr, 4);
        std::memcpy(&v[4], srcPtr, 4);
        break;
    default: {
        static const unsigned inc32table[8] = {0, 1, 2,  1,  0,  4, 4, 4};
        static const int      dec64table[8] = {0, 0, 0, -1, -4,  1, 2, 3};
        //assert(srcPtr + offset == dstPtr);
        if (offset < 8) {
            *(uint32_t*)dstPtr = 0;
            dstPtr[0] = srcPtr[0];
            dstPtr[1] = srcPtr[1];
            dstPtr[2] = srcPtr[2];
            dstPtr[3] = srcPtr[3];
            srcPtr += inc32table[offset];
            std::memcpy(dstPtr+4, srcPtr, 4);
            srcPtr -= dec64table[offset];
            dstPtr += 8;
        } else {
            std::memcpy(dstPtr, srcPtr, 8);
            dstPtr += 8;
            srcPtr += 8;
        }
 
        memcpy_beyond8(srcPtr, dstPtr, dstEnd);
 
        return;
    }
    }
 
    std::memcpy(dstPtr, v, 8);
    dstPtr += 8;
    while (dstPtr < dstEnd) {
        std::memcpy(dstPtr, v, 8);
        dstPtr += 8;
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Memory copy allowing to overwrite up to 8 Bytes beyond the end
 *
 * @param[in] srcPtr Pointer to start of source.
 * @param[in] dstPtr Pointer to start of destination.
 * @param[in] dstEnd Pointer to end of destination.
 *
 * This method is a customised variant of memcpy, which can overwrite up
 * to 8 Bytes beyond @p dstEnd. This method copies 8 Bytes a time using
 * the std::memcpy() function.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the LZ4_wildCopy8() function in the file lz4.c.
 ***************************************************************************/
void memcpy_beyond8(const uint8_t* srcPtr,
                    uint8_t*       dstPtr,
                    uint8_t*       dstEnd)
{
    // Perform memory copy
    do {
        std::memcpy(dstPtr, srcPtr, 8);
        dstPtr += 8;
        srcPtr += 8;
    } while (dstPtr < dstEnd);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Memory copy allowing to overwrite up to 32 Bytes beyond the end
 *
 * @param[in] srcPtr Pointer to start of source.
 * @param[in] dstPtr Pointer to start of destination.
 * @param[in] dstEnd Pointer to end of destination.
 *
 * This method is a customised variant of memcpy, which can overwrite up
 * to 32 Bytes beyond @p dstEnd. This method copies 16 Bytes a time using
 * the std::memcpy() function.
 *
 * The method was inspired by https://github.com/kiyo-masui/bitshuffle and
 * specifically the LZ4_wildCopy8() function in the file lz4.c.
 ***************************************************************************/
void memory_beyond32(const uint8_t* srcPtr, uint8_t* dstPtr, uint8_t* dstEnd)
{
    // Perform memory copy (in chunks of 16 Bytes)
    do {
        std::memcpy(dstPtr, srcPtr, 16);
        dstPtr += 16;
        srcPtr += 16;
        std::memcpy(dstPtr, srcPtr, 16);
        dstPtr += 16;
        srcPtr += 16;
    } while (dstPtr < dstEnd);
 
    // Return
    return;
}