hdf5io.hpp

Go to the documentation of this file.

//============================================================================
// vSMC/include/vsmc/utility/hdf5io.hpp
//----------------------------------------------------------------------------
//                         vSMC: Scalable Monte Carlo
//----------------------------------------------------------------------------
// Copyright (c) 2013-2015, Yan Zhou
// 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.
//
// 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 VSMC_UTILITY_HDF5IO_HPP
#define VSMC_UTILITY_HDF5IO_HPP

#include <vsmc/internal/common.hpp>
#include <hdf5.h>

namespace vsmc
{

namespace internal
{

template <MatrixOrder>
inline void hdf5io_matrix_dim(std::size_t, std::size_t, ::hsize_t *);

template <>
inline void hdf5io_matrix_dim<RowMajor>(
    std::size_t nrow, std::size_t ncol, ::hsize_t *dim)
{
    dim[0] = nrow;
    dim[1] = ncol;
}

template <>
inline void hdf5io_matrix_dim<ColMajor>(
    std::size_t nrow, std::size_t ncol, ::hsize_t *dim)
{
    dim[0] = ncol;
    dim[1] = nrow;
}

template <typename T>
class HDF5LoadDataPtr
{
    public:
    HDF5LoadDataPtr() : ptr_(nullptr) {}

    template <typename OutputIter>
    void set(std::size_t n, OutputIter)
    {
        data_.resize(n);
    }

    void set(std::size_t, T *ptr) { ptr_ = ptr; }

    T *get() { return ptr_ == nullptr ? data_.data() : ptr_; }

    bool is_raw_ptr() const { return ptr_ == nullptr; }

    private:
    T *ptr_;
    Vector<T> data_;
}; // class HDF5LoadDataPtr

template <typename T>
class HDF5StoreDataPtr
{
    public:
    HDF5StoreDataPtr() : ptr_(nullptr) {}

    template <typename InputIter>
    InputIter set(std::size_t n, InputIter first)
    {
        data_.resize(n);
        T *dst = data_.data();
        for (std::size_t i = 0; i != n; ++i, ++first)
            dst[i] = *first;

        return first;
    }

    T *set(std::size_t n, T *ptr)
    {
        ptr_ = ptr;
        return ptr + n;
    }

    const T *set(std::size_t n, const T *ptr)
    {
        ptr_ = ptr;
        return ptr + n;
    }

    const T *get() const { return ptr_ == nullptr ? data_.data() : ptr_; }

    bool is_raw_ptr() const { return ptr_ == nullptr; }

    private:
    const T *ptr_;
    Vector<T> data_;
}; // class HDF5StoreDataPtr

} // namespace vsmc::internal

template <typename>
inline ::hid_t hdf5io_datatype()
{
    return -1;
}

template <>
inline ::hid_t hdf5io_datatype<char>()
{
    return ::H5Tcopy(H5T_NATIVE_CHAR);
}

template <>
inline ::hid_t hdf5io_datatype<signed char>()
{
    return ::H5Tcopy(H5T_NATIVE_SCHAR);
}

template <>
inline ::hid_t hdf5io_datatype<unsigned char>()
{
    return ::H5Tcopy(H5T_NATIVE_UCHAR);
}

template <>
inline ::hid_t hdf5io_datatype<short>()
{
    return ::H5Tcopy(H5T_NATIVE_SHORT);
}

template <>
inline ::hid_t hdf5io_datatype<unsigned short>()
{
    return ::H5Tcopy(H5T_NATIVE_UCHAR);
}

template <>
inline ::hid_t hdf5io_datatype<int>()
{
    return ::H5Tcopy(H5T_NATIVE_INT);
}

template <>
inline ::hid_t hdf5io_datatype<unsigned int>()
{
    return ::H5Tcopy(H5T_NATIVE_UINT);
}

template <>
inline ::hid_t hdf5io_datatype<long>()
{
    return ::H5Tcopy(H5T_NATIVE_LONG);
}

template <>
inline ::hid_t hdf5io_datatype<unsigned long>()
{
    return ::H5Tcopy(H5T_NATIVE_ULONG);
}

template <>
inline ::hid_t hdf5io_datatype<long long>()
{
    return ::H5Tcopy(H5T_NATIVE_LLONG);
}

template <>
inline ::hid_t hdf5io_datatype<unsigned long long>()
{
    return ::H5Tcopy(H5T_NATIVE_ULLONG);
}

template <>
inline ::hid_t hdf5io_datatype<float>()
{
    return ::H5Tcopy(H5T_NATIVE_FLOAT);
}

template <>
inline ::hid_t hdf5io_datatype<double>()
{
    return ::H5Tcopy(H5T_NATIVE_DOUBLE);
}

template <>
inline ::hid_t hdf5io_datatype<long double>()
{
    return ::H5Tcopy(H5T_NATIVE_LDOUBLE);
}

inline ::hsize_t hdf5size(
    const std::string &file_name, const std::string &data_name)
{
    ::hid_t datafile =
        ::H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
    ::hid_t dataset = ::H5Dopen(datafile, data_name.c_str(), H5P_DEFAULT);
    ::hid_t dataspace = ::H5Dget_space(dataset);
    ::hid_t datatype = ::H5Dget_type(dataset);
    std::size_t n =
        static_cast<std::size_t>(H5Sget_simple_extent_npoints(dataspace));
    std::size_t b = H5Tget_size(datatype);
    std::size_t bytes = n * b;

    ::H5Tclose(datatype);
    ::H5Sclose(dataspace);
    ::H5Dclose(dataset);
    ::H5Fclose(datafile);

    return bytes;
}

template <typename T>
inline ::hsize_t hdf5size(
    const std::string &file_name, const std::string &data_name)
{
    return hdf5size(file_name, data_name) / sizeof(T);
}

template <typename T, typename OutputIter>
inline OutputIter hdf5load(const std::string &file_name,
    const std::string &data_name, OutputIter first)
{
    std::size_t n = hdf5size(file_name, data_name) / sizeof(T);
    internal::HDF5LoadDataPtr<T> data_ptr;
    data_ptr.set(n, first);
    T *data = data_ptr.get();

    ::hid_t datafile =
        ::H5Fopen(file_name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
    ::hid_t dataset = ::H5Dopen(datafile, data_name.c_str(), H5P_DEFAULT);
    ::hid_t datatype = ::H5Dget_type(dataset);
    ::H5Dread(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
    if (data_ptr.is_raw_ptr()) {
        first += n;
    } else {
        for (std::size_t i = 0; i != n; ++i, ++first)
            *first = data[i];
    }

    ::H5Tclose(datatype);
    ::H5Dclose(dataset);
    ::H5Fclose(datafile);

    return first;
}

inline void hdf5store_new(const std::string &file_name)
{
    ::hid_t datafile = ::H5Fcreate(
        file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    ::H5Fclose(datafile);
}

template <MatrixOrder Order, typename T, typename InputIter>
inline InputIter hdf5store_matrix(std::size_t nrow, std::size_t ncol,
    const std::string &file_name, const std::string &data_name,
    InputIter first, bool append = false)
{
    if (nrow == 0 || ncol == 0)
        return first;

    std::string dataset_name("/" + data_name);
    ::hsize_t dim[2];
    internal::hdf5io_matrix_dim<Order>(nrow, ncol, dim);
    internal::HDF5StoreDataPtr<T> data_ptr;
    InputIter last = data_ptr.set(nrow * ncol, first);
    const T *data = data_ptr.get();

    ::hid_t datafile;
    if (append) {
        datafile = ::H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
    } else {
        datafile = ::H5Fcreate(
            file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    }
    ::hid_t dataspace = ::H5Screate_simple(2, dim, nullptr);
    ::hid_t datatype = hdf5io_datatype<T>();
    ::hid_t dataset = ::H5Dcreate(datafile, dataset_name.c_str(), datatype,
        dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    ::H5Dwrite(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);

    ::H5Dclose(dataset);
    ::H5Tclose(datatype);
    ::H5Sclose(dataspace);
    ::H5Fclose(datafile);

    return last;
}

inline void hdf5store_list_empty(const std::string &file_name,
    const std::string &data_name, bool append = false)
{
    std::string group_name("/" + data_name);

    ::hid_t datafile;
    if (append) {
        datafile = ::H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
    } else {
        datafile = ::H5Fcreate(
            file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    }
    ::hid_t datagroup = ::H5Gcreate(
        datafile, group_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    ::H5Gclose(datagroup);
    ::H5Fclose(datafile);
}

template <typename T, typename InputIterIter, typename SInputIter>
inline void hdf5store_list(std::size_t nrow, std::size_t ncol,
    const std::string &file_name, const std::string &data_name,
    InputIterIter first, SInputIter sfirst, bool append = false)
{
    std::string group_name("/" + data_name);
    ::hsize_t dim[1] = {nrow};

    ::hid_t datafile;
    if (append) {
        datafile = ::H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
    } else {
        datafile = ::H5Fcreate(
            file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
    }
    ::hid_t datagroup = ::H5Gcreate(
        datafile, group_name.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);

    if (nrow != 0 && ncol != 0) {
        ::hid_t dataspace = ::H5Screate_simple(1, dim, nullptr);
        ::hid_t datatype = hdf5io_datatype<T>();
        internal::HDF5StoreDataPtr<T> data_ptr;
        for (std::size_t j = 0; j != ncol; ++j, ++first, ++sfirst) {
            data_ptr.set(nrow, *first);
            const T *data = data_ptr.get();
            std::string dataset_name(group_name + "/" + (*sfirst));
            ::hid_t dataset = ::H5Dcreate(datafile, dataset_name.c_str(),
                datatype, dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            ::H5Dwrite(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);
            ::H5Dclose(dataset);
        }
        ::H5Tclose(datatype);
        ::H5Sclose(dataspace);
    }

    ::H5Gclose(datagroup);
    ::H5Fclose(datafile);
}

template <typename T, typename InputIter>
inline void hdf5store_list_insert(std::size_t N, const std::string &file_name,
    const std::string &data_name, InputIter first, const std::string &vname)
{
    if (N == 0)
        return;

    std::string dataset_name("/" + data_name + "/" + vname);
    ::hsize_t dim[1] = {N};

    ::hid_t datafile = ::H5Fopen(file_name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
    ::hid_t dataspace = ::H5Screate_simple(1, dim, nullptr);
    ::hid_t datatype = hdf5io_datatype<T>();
    internal::HDF5StoreDataPtr<T> data_ptr;
    data_ptr.set(N, first);
    const T *data = data_ptr.get();
    ::hid_t dataset = ::H5Dcreate(datafile, dataset_name.c_str(), datatype,
        dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    ::H5Dwrite(dataset, datatype, H5S_ALL, H5S_ALL, H5P_DEFAULT, data);

    ::H5Dclose(dataset);
    ::H5Tclose(datatype);
    ::H5Sclose(dataspace);
    ::H5Fclose(datafile);
}

namespace internal
{

template <typename InputIter, typename... InputIters>
inline void hdf5store_list_insert_tuple(std::size_t nrow,
    const std::string &file_name, const std::string &data_name,
    const std::tuple<InputIter, InputIters...> &first, const std::string *sptr,
    std::integral_constant<std::size_t, 0>)
{
    using value_type =
        typename std::iterator_traits<typename std::tuple_element<0,
            std::tuple<InputIter, InputIters...>>::type>::value_type;
    HDF5StoreDataPtr<value_type> data_ptr;
    data_ptr.set(nrow, std::get<0>(first));
    const value_type *data = data_ptr.get();
    hdf5store_list_insert<value_type>(nrow, file_name, data_name, data, *sptr);
}

template <typename InputIter, typename... InputIters, std::size_t Pos>
inline void hdf5store_list_insert_tuple(std::size_t nrow,
    const std::string &file_name, const std::string &data_name,
    const std::tuple<InputIter, InputIters...> &first, const std::string *sptr,
    std::integral_constant<std::size_t, Pos>)
{
    using value_type =
        typename std::iterator_traits<typename std::tuple_element<Pos,
            std::tuple<InputIter, InputIters...>>::type>::value_type;
    HDF5StoreDataPtr<value_type> data_ptr;
    data_ptr.set(nrow, std::get<Pos>(first));
    const value_type *data = data_ptr.get();
    hdf5store_list_insert<value_type>(nrow, file_name, data_name, data, *sptr);
    hdf5store_list_insert_tuple(nrow, file_name, data_name, first, --sptr,
        std::integral_constant<std::size_t, Pos - 1>());
}

} // namespace vsmc::internal

template <typename SInputIter, typename InputIter, typename... InputIters>
inline void hdf5store_list(std::size_t nrow, const std::string &file_name,
    const std::string &data_name,
    const std::tuple<InputIter, InputIters...> &first, SInputIter sfirst,
    bool append = false)
{
    static constexpr std::size_t dim = sizeof...(InputIters) + 1;
    internal::HDF5StoreDataPtr<std::string> vnames;
    vnames.set(dim, sfirst);
    const std::string *sptr = vnames.get() + dim;
    hdf5store_list_empty(file_name, data_name, append);
    internal::hdf5store_list_insert_tuple(nrow, file_name, data_name, first,
        --sptr, std::integral_constant<std::size_t, dim - 1>());
}

namespace internal
{

template <typename IntType>
inline bool hdf5store_int(std::size_t n, IntType *r, std::false_type)
{
    if (sizeof(int) > sizeof(IntType))
        return true;

    bool flag = true;
    for (std::size_t i = 0; i != n; ++i) {
        if (r[i] > std::numeric_limits<int>::max VSMC_MNE()) {
            flag = false;
            break;
        }
    }

    return flag;
}

template <typename IntType>
inline bool hdf5store_int(std::size_t n, IntType *r, std::true_type)
{
    if (sizeof(int) > sizeof(IntType))
        return true;

    bool flag = true;
    for (std::size_t i = 0; i != n; ++i) {
        if (r[i] < std::numeric_limits<int>::min VSMC_MNE()) {
            flag = false;
            break;
        }
        if (r[i] > std::numeric_limits<int>::max VSMC_MNE()) {
            flag = false;
            break;
        }
    }

    return flag;
}

} // namespace vsmc::internal

template <typename T>
inline void hdf5store(const Sampler<T> &sampler, const std::string &file_name,
    const std::string &data_name, bool append = false)
{
    using size_type = typename Sampler<T>::size_type;

    std::size_t nrow = sampler.iter_size();

    if (nrow == 0)
        return;

    hdf5store_list_empty(file_name, data_name, append);

    std::size_t ncol_int = sampler.summary_header_size_int();
    Vector<std::string> header_int(ncol_int);
    Vector<size_type> data_int(nrow * ncol_int);
    sampler.summary_header_int(header_int.begin());
    sampler.template summary_data_int<ColMajor>(data_int.begin());
    bool use_int = internal::hdf5store_int(
        data_int.size(), data_int.data(), std::is_signed<size_type>());
    if (use_int) {
        Vector<int> data_int_small(data_int.size());
        std::copy(data_int.begin(), data_int.end(), data_int_small.begin());
        Vector<const int *> data_ptr_int(ncol_int);
        for (std::size_t j = 0; j != ncol_int; ++j)
            data_ptr_int[j] = data_int_small.data() + j * nrow;
        for (std::size_t j = 0; j != ncol_int; ++j)
            hdf5store_list_insert<int>(
                nrow, file_name, data_name, data_ptr_int[j], header_int[j]);
    } else {
        Vector<const size_type *> data_ptr_int(ncol_int);
        for (std::size_t j = 0; j != ncol_int; ++j)
            data_ptr_int[j] = data_int.data() + j * nrow;
        for (std::size_t j = 0; j != ncol_int; ++j)
            hdf5store_list_insert<size_type>(
                nrow, file_name, data_name, data_ptr_int[j], header_int[j]);
    }

    std::size_t ncol = sampler.summary_header_size();
    Vector<std::string> header(ncol);
    Vector<double> data(nrow * ncol);
    sampler.summary_header(header.begin());
    sampler.template summary_data<ColMajor>(data.begin());
    Vector<const double *> data_ptr(ncol);
    for (std::size_t j = 0; j != ncol; ++j)
        data_ptr[j] = data.data() + j * nrow;
    for (std::size_t j = 0; j != ncol; ++j)
        hdf5store_list_insert<double>(
            nrow, file_name, data_name, data_ptr[j], header[j]);
}

template <MatrixOrder Order, std::size_t Dim, typename T>
inline void hdf5store(const StateMatrix<Order, Dim, T> &state,
    const std::string &file_name, const std::string &data_name,
    bool append = false)
{
    hdf5store_matrix<Order, T>(
        state.size(), state.dim(), file_name, data_name, state.data(), append);
}

template <MatrixOrder Order, typename T, std::size_t StateSize,
    typename RealType, typename ID>
inline void hdf5store(const StateCL<StateSize, RealType, ID> &state,
    const std::string &file_name, const std::string &data_name,
    bool append = false)
{
    std::size_t nrow = state.size();
    std::size_t ncol = state.state_size() / sizeof(T);
    std::size_t N = nrow * ncol;
    Vector<T> data(N);
    state.manager().template read_buffer<T>(
        state.state_buffer().data(), N, data.data());
    hdf5store_matrix<Order, T>(
        nrow, ncol, file_name, data_name, data.data(), append);
}

template <typename T>
inline void hdf5store(const Particle<T> &particle,
    const std::string &file_name, const std::string &data_name,
    bool append = false)
{
    hdf5store_list_empty(file_name, data_name, append);
    hdf5store(particle.value(), file_name, data_name + "/value", true);
    hdf5store_matrix<ColMajor, double>(particle.size(), 1, file_name,
        data_name + "/weight", particle.weight().data(), true);
}

template <MatrixOrder Order, typename T, typename U>
inline void hdf5store(const Particle<U> &particle,
    const std::string &file_name, const std::string &data_name,
    bool append = false)
{
    hdf5store_list_empty(file_name, data_name, append);
    hdf5store<Order, T>(
        particle.value(), file_name, data_name + "/value", true);
    hdf5store_matrix<ColMajor, double>(particle.size(), 1, file_name,
        data_name + "/weight", particle.weight().data(), true);
}

} // namespace vsmc

#endif // VSMC_UTILITY_HDF5IO_HPP