HDF5Utils.hxx

/*
 * This file is part of the statismo library.
 *
 * Author: Marcel Luethi (marcel.luethi@unibas.ch)
 *
 * Copyright (c) 2011 University of Basel
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 * Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * Redistributions in binary form must reproduce the above copyright
 * notice, this list of conditions and the following disclaimer in the
 * documentation and/or other materials provided with the distribution.
 *
 * Neither the name of the project's author nor the names of its
 * contributors may be used to endorse or promote products derived from
 * this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
 * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
 * TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *
 */

#ifndef __HDF5_UTILS_CXX
#define __HDF5_UTILS_CXX

#include "HDF5Utils.h"

#include <algorithm>
#include <fstream>
#include <iostream>
#include <iterator>

#include "CommonTypes.h"
#include "Exceptions.h"
#include "H5Cpp.h"

namespace statismo {

inline
H5::H5File
HDF5Utils::openOrCreateFile(const std::string filename) {

    // check if file exists
    std::ifstream ifile(filename.c_str());
    H5::H5File file;

    if (!ifile) {
        // create it
        file = H5::H5File( filename.c_str(), H5F_ACC_EXCL);
    } else {
        // open it
        file = H5::H5File( filename.c_str(), H5F_ACC_RDWR);
    }
    return file;
}



inline
H5::Group
HDF5Utils::openPath(H5::H5File& file, const std::string& path, bool createPath) {
    H5::Group group;

    // take the first part of the path
    size_t curpos = 1;
    size_t nextpos = path.find_first_of("/", curpos);
    H5::Group g = file.openGroup("/");

    std::string name = path.substr(curpos, nextpos-1);

    while (curpos != std::string::npos && name != "") {

        if (existsObjectWithName(g, name)) {
            g = g.openGroup(name);
        } else {
            if (createPath) {
                g = g.createGroup(name);
            } else {
                std::string msg = std::string("the path ") +path +" does not exist";
                throw StatisticalModelException(msg.c_str());
            }
        }

        curpos = nextpos+1;
        nextpos = path.find_first_of("/", curpos);
        if ( nextpos != std::string::npos )
            name = path.substr(curpos, nextpos-curpos);
        else
            name = path.substr(curpos);
    }

    return g;
}

template <class T>
inline
void HDF5Utils::readMatrixOfType(const H5::CommonFG& fg, const char* name, typename GenericEigenType<T>::MatrixType& matrix) {
    throw StatisticalModelException("Invalid type proided for writeMatrixOfType");
}

template <>
inline
void HDF5Utils::readMatrixOfType<unsigned int>(const H5::CommonFG& fg, const char* name, GenericEigenType<unsigned int>::MatrixType& matrix) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[2];
    ds.getSpace().getSimpleExtentDims(dims, NULL);

    // simply read the whole dataspace
    matrix.resize(dims[0], dims[1]);
    ds.read(matrix.data(), H5::PredType::NATIVE_UINT);
}

template <>
inline
void HDF5Utils::readMatrixOfType<float>(const H5::CommonFG& fg, const char* name, GenericEigenType<float>::MatrixType& matrix) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[2];
    ds.getSpace().getSimpleExtentDims(dims, NULL);

    // simply read the whole dataspace
    matrix.resize(dims[0], dims[1]);
    ds.read(matrix.data(), H5::PredType::NATIVE_FLOAT);
}

template <>
inline
void HDF5Utils::readMatrixOfType<double>(const H5::CommonFG& fg, const char* name, GenericEigenType<double>::MatrixType& matrix) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[2];
    ds.getSpace().getSimpleExtentDims(dims, NULL);

    // simply read the whole dataspace
    matrix.resize(dims[0], dims[1]);
    ds.read(matrix.data(), H5::PredType::NATIVE_DOUBLE);
}


inline
void HDF5Utils::readMatrix(const H5::CommonFG& fg, const char* name, MatrixType& matrix) {
    readMatrixOfType<ScalarType>(fg, name, matrix);
}


inline
void HDF5Utils::readMatrix(const H5::CommonFG& fg, const char* name, unsigned maxNumColumns, MatrixType& matrix) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[2];
    ds.getSpace().getSimpleExtentDims(dims, NULL);

    hsize_t nRows = dims[0]; // take the number of rows defined in the hdf5 file
    hsize_t nCols = std::min(dims[1], static_cast<hsize_t>(maxNumColumns)); // take the number of cols provided by the user

    hsize_t offset[2] = {0,0};   // hyperslab offset in the file
    hsize_t count[2];
    count[0] = nRows;
    count[1] =  nCols;

    H5::DataSpace dataspace = ds.getSpace();
    dataspace.selectHyperslab( H5S_SELECT_SET, count, offset );

    /* Define the memory dataspace. */
    hsize_t     dimsm[2];
    dimsm[0] = nRows;
    dimsm[1] = nCols;
    H5::DataSpace memspace( 2, dimsm );

    /* Define memory hyperslab. */
    hsize_t      offset_out[2] = {0, 0};       // hyperslab offset in memory
    hsize_t      count_out[2];        // size of the hyperslab in memory

    count_out[0]  = nRows;
    count_out[1] = nCols;
    memspace.selectHyperslab( H5S_SELECT_SET, count_out, offset_out );

    matrix.resize(nRows, nCols);
    ds.read(matrix.data(), H5::PredType::NATIVE_FLOAT, memspace, dataspace);

}

template <class T>
inline
H5::DataSet HDF5Utils::writeMatrixOfType(const H5::CommonFG& fg, const char* name, const typename GenericEigenType<T>::MatrixType& matrix) {
    throw StatisticalModelException("Invalid type proided for writeMatrixOfType");
}

template <>
inline
H5::DataSet HDF5Utils::writeMatrixOfType<unsigned int>(const H5::CommonFG& fg, const char* name, const GenericEigenType<unsigned int>::MatrixType& matrix) {
    // HDF5 does not like empty matrices.
    //
    if (matrix.rows() == 0 || matrix.cols() == 0) {
        throw StatisticalModelException("Empty matrix provided to writeMatrix");
    }

    hsize_t dims[2] = {matrix.rows(), matrix.cols()};
    H5::DataSet ds = fg.createDataSet( name, H5::PredType::NATIVE_UINT, H5::DataSpace(2, dims));
    ds.write( matrix.data(), H5::PredType::NATIVE_UINT );
    return ds;
}

template <>
inline
H5::DataSet HDF5Utils::writeMatrixOfType<float>(const H5::CommonFG& fg, const char* name, const GenericEigenType<float>::MatrixType& matrix) {
    // HDF5 does not like empty matrices.
    //
    if (matrix.rows() == 0 || matrix.cols() == 0) {
        throw StatisticalModelException("Empty matrix provided to writeMatrix");
    }

    hsize_t dims[2] = {matrix.rows(), matrix.cols()};
    H5::DataSet ds = fg.createDataSet( name, H5::PredType::NATIVE_FLOAT, H5::DataSpace(2, dims));
    ds.write( matrix.data(), H5::PredType::NATIVE_FLOAT );
    return ds;
}

template <>
inline
H5::DataSet HDF5Utils::writeMatrixOfType<double>(const H5::CommonFG& fg, const char* name, const GenericEigenType<double>::MatrixType& matrix) {
    // HDF5 does not like empty matrices.
    //
    if (matrix.rows() == 0 || matrix.cols() == 0) {
        throw StatisticalModelException("Empty matrix provided to writeMatrix");
    }

    hsize_t dims[2] = {matrix.rows(), matrix.cols()};
    H5::DataSet ds = fg.createDataSet( name, H5::PredType::NATIVE_DOUBLE, H5::DataSpace(2, dims));
    ds.write( matrix.data(), H5::PredType::NATIVE_DOUBLE );
    return ds;
}


inline
H5::DataSet HDF5Utils::writeMatrix(const H5::CommonFG& fg, const char* name, const MatrixType& matrix) {
    return writeMatrixOfType<ScalarType>(fg, name, matrix);
}


template <class T>
inline
void HDF5Utils::readVectorOfType(const H5::CommonFG& fg, const char* name, typename GenericEigenType<T>::VectorType& vector) {
    throw StatisticalModelException("Invalid type proided for readVectorOfType");
}

template <>
inline
void HDF5Utils::readVectorOfType<double>(const H5::CommonFG& fg, const char* name,  GenericEigenType<double>::VectorType& vector) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[1];
    ds.getSpace().getSimpleExtentDims(dims, NULL);
    vector.resize(dims[0], 1);
    ds.read(vector.data(), H5::PredType::NATIVE_DOUBLE);
}

template <>
inline
void HDF5Utils::readVectorOfType<float>(const H5::CommonFG& fg, const char* name,  GenericEigenType<float>::VectorType& vector) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[1];
    ds.getSpace().getSimpleExtentDims(dims, NULL);
    vector.resize(dims[0], 1);
    ds.read(vector.data(), H5::PredType::NATIVE_FLOAT);
}

template <>
inline
void HDF5Utils::readVectorOfType<int>(const H5::CommonFG& fg, const char* name,  GenericEigenType<int>::VectorType& vector) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[1];
    ds.getSpace().getSimpleExtentDims(dims, NULL);
    vector.resize(dims[0], 1);
    ds.read(vector.data(), H5::PredType::NATIVE_INT);
}

inline
void HDF5Utils::readVector(const H5::CommonFG& fg, const char* name, VectorType& vector) {
    readVectorOfType<ScalarType>(fg, name, vector);
}


inline
void HDF5Utils::readVector(const H5::CommonFG& fg, const char* name, unsigned maxNumElements, VectorType& vector) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[1];
    ds.getSpace().getSimpleExtentDims(dims, NULL);

    hsize_t nElements = std::min(dims[0], static_cast<hsize_t>(maxNumElements)); // take the number of rows defined in the hdf5 file

    hsize_t offset[1] = {0};   // hyperslab offset in the file
    hsize_t count[1];
    count[0] = nElements;

    H5::DataSpace dataspace = ds.getSpace();
    dataspace.selectHyperslab( H5S_SELECT_SET, count, offset );

    /* Define the memory dataspace. */
    hsize_t     dimsm[1];
    dimsm[0] = nElements;
    H5::DataSpace memspace( 1, dimsm );

    /* Define memory hyperslab. */
    hsize_t      offset_out[1] = {0};       // hyperslab offset in memory
    hsize_t      count_out[1];        // size of the hyperslab in memory

    count_out[0]  = nElements;
    memspace.selectHyperslab( H5S_SELECT_SET, count_out, offset_out );

    vector.resize(nElements);
    ds.read(vector.data(), H5::PredType::NATIVE_FLOAT, memspace, dataspace);
}




template <class T>
inline
H5::DataSet HDF5Utils::writeVectorOfType(const H5::CommonFG& fg, const char* name, const typename GenericEigenType<T>::VectorType& vector) {
    throw StatisticalModelException("Invalid type provided for writeVectorOfType");
}

template <>
inline
H5::DataSet HDF5Utils::writeVectorOfType<double>(const H5::CommonFG& fg, const char* name, const GenericEigenType<double>::VectorType& vector) {
    hsize_t dims[1] = {vector.size()};
    H5::DataSet ds = fg.createDataSet( name, H5::PredType::NATIVE_DOUBLE, H5::DataSpace(1, dims));
    ds.write( vector.data(), H5::PredType::NATIVE_DOUBLE );
    return ds;
}

template <>
inline
H5::DataSet HDF5Utils::writeVectorOfType<float>(const H5::CommonFG& fg, const char* name, const GenericEigenType<float>::VectorType& vector) {
    hsize_t dims[1] = {vector.size()};
    H5::DataSet ds = fg.createDataSet( name, H5::PredType::NATIVE_FLOAT, H5::DataSpace(1, dims));
    ds.write( vector.data(), H5::PredType::NATIVE_FLOAT );
    return ds;
}

template <>
inline
H5::DataSet HDF5Utils::writeVectorOfType<int>(const H5::CommonFG& fg, const char* name, const GenericEigenType<int>::VectorType& vector) {
    hsize_t dims[1] = {vector.size()};
    H5::DataSet ds = fg.createDataSet( name, H5::PredType::NATIVE_INT, H5::DataSpace(1, dims));
    ds.write( vector.data(), H5::PredType::NATIVE_INT );
    return ds;
}

inline
H5::DataSet HDF5Utils::writeVector(const H5::CommonFG& fg, const char* name, const VectorType& vector) {
    return writeVectorOfType<ScalarType>(fg, name, vector);
}


inline
H5::DataSet HDF5Utils::writeString(const H5::CommonFG& fg, const char* name, const std::string& s) {
    H5::StrType fls_type(H5::PredType::C_S1, s.length() + 1); // + 1 for trailing zero
    H5::DataSet ds = fg.createDataSet(name, fls_type, H5::DataSpace(H5S_SCALAR));
    ds.write(s, fls_type);
    return ds;
}


inline
std::string
HDF5Utils::readString(const H5::CommonFG& fg, const char* name) {
    H5std_string outputString;
    H5::DataSet ds = fg.openDataSet(name);
    ds.read(outputString, ds.getStrType());
    return outputString;
}

inline
void HDF5Utils::writeStringAttribute(const H5::H5Object& fg, const char* name, const std::string& s) {
    H5::StrType strdatatype(H5::PredType::C_S1, s.length() + 1 ); // + 1 for trailing 0
    H5::Attribute att = fg.createAttribute(name, strdatatype, H5::DataSpace(H5S_SCALAR));
    att.write(strdatatype, s);
    att.close();
}


inline
std::string
HDF5Utils::readStringAttribute(const H5::H5Object& fg, const char* name) {
    H5std_string outputString;

    H5::Attribute myatt_out = fg.openAttribute(name);
    myatt_out.read(myatt_out.getStrType(), outputString);
    return outputString;
}

inline
void HDF5Utils::writeIntAttribute(const H5::H5Object& fg, const char* name, int value) {
    H5::IntType int_type(H5::PredType::NATIVE_INT32);
    H5::DataSpace att_space(H5S_SCALAR);
    H5::Attribute att = fg.createAttribute(name, int_type, att_space );
    att.write( int_type, &value);
    att.close();
}

inline
int
HDF5Utils::readIntAttribute(const H5::H5Object& fg, const char* name) {
    H5::IntType fls_type(H5::PredType::NATIVE_INT32);
    int value = 0;
    H5::Attribute myatt_out = fg.openAttribute(name);
    myatt_out.read(fls_type, &value);
    return value;
}


inline
H5::DataSet HDF5Utils::writeInt(const H5::CommonFG& fg, const char* name, int value) {
    H5::IntType fls_type(H5::PredType::NATIVE_INT32); // 0 is a dummy argument
    H5::DataSet ds = fg.createDataSet(name, fls_type, H5::DataSpace(H5S_SCALAR));
    ds.write(&value, fls_type);
    return ds;
}

inline
int HDF5Utils::readInt(const H5::CommonFG& fg, const char* name) {
    H5::IntType fls_type(H5::PredType::NATIVE_INT32);
    H5::DataSet ds = fg.openDataSet( name );

    int value = 0;
    ds.read(&value, fls_type);
    return value;
}

inline
H5::DataSet HDF5Utils::writeFloat(const H5::CommonFG& fg, const char* name, float value) {
    H5::FloatType fls_type(H5::PredType::NATIVE_FLOAT); // 0 is a dummy argument
    H5::DataSet ds = fg.createDataSet(name, fls_type, H5::DataSpace(H5S_SCALAR));
    ds.write(&value, fls_type);
    return ds;
}

inline
float HDF5Utils::readFloat(const H5::CommonFG& fg, const char* name) {
    H5::FloatType fls_type(H5::PredType::NATIVE_FLOAT);
    H5::DataSet ds = fg.openDataSet( name );

    float value = 0;
    ds.read(&value, fls_type);
    return value;
}

inline
void HDF5Utils::getFileFromHDF5(const H5::CommonFG& fg, const char* name, const char* filename) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[1];
    ds.getSpace().getSimpleExtentDims(dims, NULL);
    std::vector<char> buffer(dims[0]);
    if(!buffer.empty()) ds.read(&buffer[0], H5::PredType::NATIVE_CHAR);

    typedef std::ostream_iterator<char> ostream_iterator;
    std::ofstream ofile(filename, std::ios::binary);
    if (!ofile) {
        std::string s= std::string("could not open file ") +filename;
        throw StatisticalModelException(s.c_str());
    }

    std::copy(buffer.begin(), buffer.end(), ostream_iterator(ofile));
    ofile.close();
}

inline
void
HDF5Utils::dumpFileToHDF5( const char* filename, const H5::CommonFG& fg, const char* name) {

    typedef std::istream_iterator<char> istream_iterator;

    std::ifstream ifile(filename, std::ios::binary);
    if (!ifile) {
        std::string s= std::string("could not open file ") +filename;
        throw StatisticalModelException(s.c_str());
    }

    std::vector<char> buffer;
    ifile >> std::noskipws;
    std::copy(istream_iterator(ifile), istream_iterator(), std::back_inserter(buffer));

    ifile.close();

    hsize_t dims[] = {buffer.size()};
    H5::DataSet ds = fg.createDataSet( name,  H5::PredType::NATIVE_CHAR, H5::DataSpace(1, dims));
    ds.write( &buffer[0], H5::PredType::NATIVE_CHAR );

}

template<typename T>
inline
void
HDF5Utils::readArray(const H5::CommonFG& fg, const char* name, std::vector<T> & array) {
    throw StatisticalModelException( "not implemented" );
}


template<typename T>
inline
H5::DataSet
HDF5Utils::writeArray(const H5::CommonFG& fg, const char* name, std::vector<T> const& array) {
    throw StatisticalModelException( "not implemented" );
}

template<>
inline
void
HDF5Utils::readArray(const H5::CommonFG& fg, const char* name, std::vector<int> & array) {
    H5::DataSet ds = fg.openDataSet( name );
    hsize_t dims[1];
    ds.getSpace().getSimpleExtentDims(dims, NULL);
    array.resize(dims[0]);
    ds.read( &array[0], H5::PredType::NATIVE_INT32);
}

template<>
inline
H5::DataSet
HDF5Utils::writeArray(const H5::CommonFG& fg, const char* name, std::vector<int> const& array) {
    hsize_t dims[1] = {array.size()};
    H5::DataSet ds = fg.createDataSet( name, H5::PredType::NATIVE_INT32, H5::DataSpace(1, dims));
    ds.write( &array[0], H5::PredType::NATIVE_INT32 );
    return ds;
}

inline
bool
HDF5Utils::existsObjectWithName(const H5::CommonFG& fg, const std::string& name) {
    for (hsize_t i = 0; i < fg.getNumObjs(); ++i) {
        std::string objname=  fg.getObjnameByIdx(i);
        if (objname == name) {
            return true;
        }
    }
    return false;
}

} //namespace statismo

#endif