diff --git a/triqs/arrays/h5.hpp b/triqs/arrays/h5.hpp index c22c5727..94feac14 100644 --- a/triqs/arrays/h5.hpp +++ b/triqs/arrays/h5.hpp @@ -18,7 +18,6 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_H5_H -#define TRIQS_ARRAYS_H5_H +#pragma once #include "./h5/simple_read_write.hpp" -#endif +#include "./h5/array_of_non_basic.hpp" diff --git a/triqs/arrays/h5/array_of_non_basic.hpp b/triqs/arrays/h5/array_of_non_basic.hpp new file mode 100644 index 00000000..e1a0cfdc --- /dev/null +++ b/triqs/arrays/h5/array_of_non_basic.hpp @@ -0,0 +1,111 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2011-2013 by O. Parcollet + * + * TRIQS 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. + * + * TRIQS 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 + * TRIQS. If not, see . + * + ******************************************************************************/ +#pragma once +#include +#include +#include +#include "../cache.hpp" + +namespace triqs { +namespace arrays { + + // details for generic save/read of arrays. + namespace h5_impl { + inline std::string _h5_name() { return ""; } + + template std::string _h5_name(T0 const& t0, Ts const&... ts) { + auto r = std::to_string(t0); + auto r1 = _h5_name(ts...); + if (r1 != "") r += "_" + r1; + return r; + } + +#ifndef __cpp_generic_lambdas + template struct _save_lambda { + ArrayType const& a; + h5::group g; + template void operator()(Is const&... is) const { h5_write(g, _h5_name(is...), a(is...)); } + }; + + template struct _load_lambda { + ArrayType& a; + h5::group g; + template void operator()(Is const&... is) { h5_read(g, _h5_name(is...), a(is...)); } + }; +#endif + } // details + + /* + * Write an array or a view into an hdf5 file when type is not fundamental + * ArrayType The type of the array/matrix/vector, etc.. + * g The h5 group + * name The name of the hdf5 array in the file/group where the stack will be stored + * A The array to be stored + * The HDF5 exceptions will be caught and rethrown as TRIQS_RUNTIME_ERROR (with a full stackstrace, cf triqs doc). + */ + template + std::c14::enable_if_t::value && !has_scalar_or_string_value_type::value> + h5_write(h5::group gr, std::string name, ArrayType const& a) { + if (a.is_empty()) TRIQS_RUNTIME_ERROR << " Can not save an empty array into hdf5"; + auto gr2 = gr.create_group(name); + gr2.write_triqs_hdf5_data_scheme(a); + // save the shape + array sha(ArrayType::rank); + for (int u = 0; u < ArrayType::rank; ++u) sha(u) = a.shape()[u]; + h5_write(gr2, "shape", sha); +#ifndef __cpp_generic_lambdas + foreach(a, h5_impl::_save_lambda{a, gr2}); +#else + foreach(a, [&](auto... is) { h5_write(gr2, h5_impl::_h5_name(is...), a(is...)); }); +#endif + } + + /* + * Read an array or a view from an hdf5 file when type is not fundamental + * ArrayType The type of the array/matrix/vector, etc.. + * g The h5 group + * name The name of the hdf5 array in the file/group where the stack will be stored + * A The array to be stored + * The HDF5 exceptions will be caught and rethrown as TRIQS_RUNTIME_ERROR (with a full stackstrace, cf triqs doc). + */ + template + std::c14::enable_if_t::value && !has_scalar_or_string_value_type::value> + h5_read(h5::group gr, std::string name, ArrayType& a) { + static_assert(!std::is_const::value, "Can not read in const object"); + auto gr2 = gr.open_group(name); + // TODO checking scheme... + // load the shape + auto sha2 = a.shape(); + array sha; + h5_read(gr2, "shape", sha); + if (first_dim(sha) != sha2.size()) + TRIQS_RUNTIME_ERROR << " array> load : rank mismatch. Expected " << sha2.size()<< " Got " << first_dim(sha); + for (int u = 0; u < sha2.size(); ++u) sha2[u] = sha(u); + if (a.shape() != sha2) a.resize(sha2); +#ifndef __cpp_generic_lambdas + foreach(a, h5_impl::_load_lambda{a, gr2}); +#else + foreach(a, [&](auto... is) { h5_read(gr2, h5_impl::_h5_name(is...), a(is...)); }); +#endif + } +} +} + diff --git a/triqs/arrays/h5/array_stack.hpp b/triqs/arrays/h5/array_stack.hpp index df49af21..28d0b47d 100644 --- a/triqs/arrays/h5/array_stack.hpp +++ b/triqs/arrays/h5/array_stack.hpp @@ -27,6 +27,18 @@ namespace triqs { namespace arrays { + // to be cleaned + namespace h5_impl { + template void* __get_array_data_ptr(A& x) { return h5::get_data_ptr(&(x.storage()[0])); } + + H5::DataSpace data_space_impl(array_stride_info info, bool is_complex); + + template H5::DataSpace data_space(ArrayType const& A) { + if (!A.indexmap().is_contiguous()) TRIQS_RUNTIME_ERROR << " h5 : internal error : array not contiguous"; + return data_space_impl(array_stride_info{A}, triqs::is_complex::value); + } + } + /// The implementation class template class array_stack_impl { static const size_t dim = R; diff --git a/triqs/arrays/h5/simple_read_write.cpp b/triqs/arrays/h5/simple_read_write.cpp new file mode 100644 index 00000000..e7ac0e2f --- /dev/null +++ b/triqs/arrays/h5/simple_read_write.cpp @@ -0,0 +1,150 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2011-2014 by O. Parcollet + * + * TRIQS 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. + * + * TRIQS 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 + * TRIQS. If not, see . + * + ******************************************************************************/ +#include "./simple_read_write.hpp" + +using dcomplex = std::complex; +namespace triqs { +namespace arrays { + namespace h5_impl { + + // the dataspace corresponding to the array. Contiguous data only... + H5::DataSpace data_space_impl(array_stride_info info, bool is_complex) { + hsize_t L[info.R], S[info.R]; + for (int u = 0; u < info.R; ++u) { + if (info.strides[u] <= 0) TRIQS_RUNTIME_ERROR << " negative strides not permitted in h5"; + S[u] = 1; + L[u] = info.lengths[u]; + } + return h5::dataspace_from_LS(info.R, is_complex, L, L, S); + } + + /// --------------------------- WRITE --------------------------------------------- + + template void write_array_impl(h5::group g, std::string const& name, const T* start, array_stride_info info) { + static_assert(!std::is_base_of::value, " Not implemented"); // 1d is below + bool is_complex = triqs::is_complex::value; + try { + H5::DataSet ds = g.create_dataset(name, h5::data_type_file(), data_space_impl(info, is_complex)); + ds.write(h5::get_data_ptr(start), h5::data_type_memory(), data_space_impl(info, is_complex)); + // if complex, to be python compatible, we add the __complex__ attribute + if (is_complex) h5::write_string_attribute(&ds, "__complex__", "1"); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + + template void write_array_impl(h5::group g, std::string const& name, const int* start, array_stride_info info); + template void write_array_impl(h5::group g, std::string const& name, const long* start, array_stride_info info); + template void write_array_impl(h5::group g, std::string const& name, const double* start, array_stride_info info); + template void write_array_impl(h5::group g, std::string const& name, const dcomplex* start, array_stride_info info); + + + // overload : special treatment for arrays of strings (one dimension only). + void write_array(h5::group g, std::string const& name, vector_const_view V) { + std::vector tmp(V.size()); + std::copy(begin(V), end(V), begin(tmp)); + h5_write(g, name, tmp); + } + + void write_array(h5::group g, std::string const& name, array_const_view V) { + std::vector tmp(first_dim(V)); + std::copy(begin(V), end(V), begin(tmp)); + h5_write(g, name, tmp); + } + + /// --------------------------- READ --------------------------------------------- + + +/* template void read_array(h5::group g, std::string const& name, ArrayType1&& A, bool C_reorder = true) { + typedef typename std::remove_reference::type ArrayType; + static_assert(!std::is_base_of::value, " Not implemented"); // 1d is below + try { + H5::DataSet ds = g.open_dataset(name); + H5::DataSpace dataspace = ds.getSpace(); + static const unsigned int Rank = ArrayType::rank + (triqs::is_complex::value ? 1 : 0); + int rank = dataspace.getSimpleExtentNdims(); + if (rank != Rank) + TRIQS_RUNTIME_ERROR << "triqs::array::h5::read. Rank mismatch : the array has rank = " << Rank + << " while the array stored in the hdf5 file has rank = " << rank; + mini_vector dims_out; + dataspace.getSimpleExtentDims(&dims_out[0], NULL); + mini_vector d2; + for (size_t u = 0; u < ArrayType::rank; ++u) d2[u] = dims_out[u]; + resize_or_check(A, d2); + if (C_reorder) { + read_array(g, name, make_cache(A).view(), false); + //read_array(g, name, cache(A).view(), false); + } else + ds.read(__get_array_data_ptr(A), h5::data_type_memory(), data_space(A), dataspace); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } +*/ + + std::vector get_array_lengths(int R, h5::group g, std::string const& name, bool is_complex) { + try { + H5::DataSet ds = g.open_dataset(name); + H5::DataSpace dataspace = ds.getSpace(); + int Rank = R + (is_complex ? 1 : 0); + int rank = dataspace.getSimpleExtentNdims(); + if (rank != Rank) + TRIQS_RUNTIME_ERROR << "triqs::array::h5::read. Rank mismatch : the array has rank = " << Rank + << " while the array stored in the hdf5 file has rank = " << rank; + std::vector d2(R); + hsize_t dims_out[rank]; + dataspace.getSimpleExtentDims(&dims_out[0], NULL); + for (size_t u = 0; u < R; ++u) d2[u] = dims_out[u]; + return d2; + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + + template void read_array_impl(h5::group g, std::string const& name, T* start, array_stride_info info) { + static_assert(!std::is_base_of::value, " Not implemented"); // 1d is below + bool is_complex = triqs::is_complex::value; + try { + H5::DataSet ds = g.open_dataset(name); + H5::DataSpace dataspace = ds.getSpace(); + ds.read(h5::get_data_ptr(start), h5::data_type_memory(), data_space_impl(info, is_complex), dataspace); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + + template void read_array_impl(h5::group g, std::string const& name, int* start, array_stride_info info); + template void read_array_impl(h5::group g, std::string const& name, long* start, array_stride_info info); + template void read_array_impl(h5::group g, std::string const& name, double* start, array_stride_info info); + template void read_array_impl(h5::group g, std::string const& name, dcomplex* start, array_stride_info info); + + void read_array(h5::group g, std::string const& name, arrays::vector& V) { + std::vector tmp; + h5_read(g, name, tmp); + V.resize(tmp.size()); + std::copy(begin(tmp), end(tmp), begin(V)); + } + + // I can not use the generic code, just because the resize of the array take a shape, not a size_t as std::vector and vector + void read_array(h5::group f, std::string const& name, arrays::array& V) { + arrays::vector res; + read_array(f, name, res); + V = res; + } + } +} +} diff --git a/triqs/arrays/h5/simple_read_write.hpp b/triqs/arrays/h5/simple_read_write.hpp index 7dbc74d0..74b295b7 100644 --- a/triqs/arrays/h5/simple_read_write.hpp +++ b/triqs/arrays/h5/simple_read_write.hpp @@ -18,8 +18,7 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_ARRAYS_H5_LOWLEVEL_H -#define TRIQS_ARRAYS_H5_LOWLEVEL_H +#pragma once #include #include #include @@ -29,24 +28,18 @@ namespace triqs { namespace arrays { namespace h5_impl { - template const void* __get_array_data_cptr(A const& a) { return h5::get_data_ptr(&(a.storage()[0])); } - template void* __get_array_data_ptr(A& x) { return h5::get_data_ptr(&(x.storage()[0])); } - - // the dataspace corresponding to the array. Contiguous data only... - template H5::DataSpace data_space(ArrayType const& A) { - if (!A.indexmap().is_contiguous()) TRIQS_RUNTIME_ERROR << " h5 : internal error : array not contiguous"; - static const unsigned int R = ArrayType::rank; - mini_vector S; - auto const& S1(A.indexmap().strides()); - for (int u = 0; u < R; ++u) { - if (S1[u] <= 0) TRIQS_RUNTIME_ERROR << " negative strides not permitted in h5"; - S[u] = 1; + struct array_stride_info { + int R; + size_t const* lengths; + long const* strides; + template explicit array_stride_info(A && a) { + R = std::c14::decay_t::rank; + lengths = a.indexmap().domain().lengths().ptr(); + strides = a.indexmap().strides().ptr(); } - static const bool is_complex = triqs::is_complex::value; - return h5::dataspace_from_LS(A.indexmap().domain().lengths(), A.indexmap().domain().lengths(), S); - } + }; - /******************** resize or check the size ****************************************************/ + /******************** resize or check the size ****************************************************/ template ENABLE_IF(is_amv_value_class) resize_or_check(A& a, mini_vector const& dimsf) { a.resize(indexmaps::cuboid::domain_t(dimsf)); @@ -60,92 +53,44 @@ namespace arrays { } /*********************************** WRITE array ****************************************************************/ - /* - * Write an array or a view into an hdf5 file - * - * f The h5 file or group of type H5::H5File or H5::Group - * name The name of the hdf5 array in the file/group where the stack will be stored - * A The array to be stored - * C_reorder bool If true [default] the data will be stored in C order in the hdf5, hence making a temporary - * cache of the data to reorder them in memory. - * If false, the array is stored as it [if you know what you are doing] - * The HDF5 exceptions will be caught and rethrown as TRIQS_RUNTIME_ERROR (with a full stackstrace, cf triqs doc). - */ - template - void write_array_impl(h5::group g, std::string const& name, array_const_view A, bool C_reorder) { - static_assert(!std::is_base_of::value, " Not implemented"); // 1d is below - if (C_reorder) { - write_array_impl(g, name, make_const_cache(A).view(), false); - return; - } - try { - H5::DataSet ds = g.create_dataset(name, h5::data_type_file(), data_space(A)); - ds.write(__get_array_data_cptr(A), h5::data_type_memory(), data_space(A)); - // if complex, to be python compatible, we add the __complex__ attribute - if (triqs::is_complex::value) h5::write_string_attribute(&ds, "__complex__", "1"); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } + + template void write_array_impl(h5::group g, std::string const& name, const T* start, array_stride_info info); template void write_array(h5::group g, std::string const& name, A const& a, bool C_reorder = true) { - write_array_impl(g, name, typename A::const_view_type{a}, C_reorder); - } - - /*********************************** READ array ****************************************************************/ - /* - * Read an array or a view from an hdf5 file - * ArrayType The type of the array/matrix/vector, etc.. - * f The h5 file or group of type H5::H5File or H5::Group - * name The name of the hdf5 array in the file/group where the stack will be stored - * A The array to be stored - * C_reorder bool If true [default] the data will be stored in C order in the hdf5, hence making a temporary - * cache of the data to reorder them in memory. If false, the array is stored as it [if you know what you are doing] - * The HDF5 exceptions will be caught and rethrown as TRIQS_RUNTIME_ERROR (with a full stackstrace, cf triqs doc). - */ - template void read_array(h5::group g, std::string const& name, ArrayType1&& A, bool C_reorder = true) { - typedef typename std::remove_reference::type ArrayType; - static_assert(!std::is_base_of::value, " Not implemented"); // 1d is below - try { - H5::DataSet ds = g.open_dataset(name); - H5::DataSpace dataspace = ds.getSpace(); - static const unsigned int Rank = ArrayType::rank + (triqs::is_complex::value ? 1 : 0); - int rank = dataspace.getSimpleExtentNdims(); - if (rank != Rank) - TRIQS_RUNTIME_ERROR << "triqs::array::h5::read. Rank mismatch : the array has rank = " << Rank - << " while the array stored in the hdf5 file has rank = " << rank; - mini_vector dims_out; - dataspace.getSimpleExtentDims(&dims_out[0], NULL); - mini_vector d2; - for (size_t u = 0; u < ArrayType::rank; ++u) d2[u] = dims_out[u]; - resize_or_check(A, d2); - if (C_reorder) { - read_array(g, name, make_cache(A).view(), false); - //read_array(g, name, cache(A).view(), false); - } else - ds.read(__get_array_data_ptr(A), h5::data_type_memory(), data_space(A), dataspace); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + if (C_reorder) { + auto b = make_const_cache(a).view(); + write_array_impl(g, name, b.data_start(), array_stride_info{b}); + } else + write_array_impl(g, name, a.data_start(), array_stride_info{a}); } // overload : special treatment for arrays of strings (one dimension only). - inline void write_array(h5::group f, std::string const& name, vector_const_view V) { - h5::detail::write_1darray_vector_of_string_impl(f, name, V); + void write_array(h5::group g, std::string const& name, vector_const_view V); + void write_array(h5::group g, std::string const& name, array_const_view V); + + /*********************************** READ array ****************************************************************/ + + std::vector get_array_lengths(int R, h5::group g, std::string const& name, bool is_complex); + template void read_array_impl(h5::group g, std::string const& name, T* start, array_stride_info info); + + template void read_array(h5::group g, std::string const& name, A&& a, bool C_reorder = true) { + resize_or_check(a, get_array_lengths(a.rank, g, name, triqs::is_complex::value_type>::value)); + if (C_reorder) { + { + // read_array(g, name, make_cache(a).view(), false); + // return ; + auto b = make_cache(a); + // auto b = make_cache(a).view(); + read_array_impl(g, name, b.view().data_start(), array_stride_info{b.view()}); + // read_array_impl(g, name, b.data_start(), array_stride_info{b}); + } + } else + read_array_impl(g, name, a.data_start(), array_stride_info{a}); } - inline void write_array(h5::group f, std::string const& name, array_const_view V) { - write_array(f, name, vector_const_view(V)); - } - - inline void read_array(h5::group f, std::string const& name, arrays::vector& V) { - h5::detail::read_1darray_vector_of_string_impl(f, name, V); - } - - // I can not use the generic code, just because the resize of the array take a shape, not a size_t as std::vector and vector - inline void read_array(h5::group f, std::string const& name, arrays::array& V) { - arrays::vector res; - read_array(f, name, res); - V = res; - } + // overload : special treatment for arrays of strings (one dimension only). + void read_array(h5::group g, std::string const& name, arrays::vector& V); + void read_array(h5::group f, std::string const& name, arrays::array& V); } // namespace h5_impl @@ -162,9 +107,7 @@ namespace arrays { template TYPE_ENABLE_IFC(std::string, is_amv_value_or_view_class::value) get_triqs_hdf5_data_scheme(ArrayType const&) { using triqs::get_triqs_hdf5_data_scheme; // for the basic types, not found by ADL - std::stringstream fs; - fs << "array<" << get_triqs_hdf5_data_scheme(typename ArrayType::value_type()) << "," << ArrayType::rank << ">"; - return fs.str(); + return "array<" + get_triqs_hdf5_data_scheme(typename ArrayType::value_type()) + "," + std::to_string(ArrayType::rank) + ">"; } /* @@ -196,86 +139,5 @@ namespace arrays { h5_impl::write_array(g, name, array_const_view(A)); } - // details for generic save/read of arrays. - namespace h5_impl { - inline std::string _h5_name() { return ""; } - - template std::string _h5_name(T0 const& t0, Ts const&... ts) { - auto r = std::to_string(t0); - auto r1 = _h5_name(ts...); - if (r1 != "") r += "_" + r1; - return r; - } - -#ifndef __cpp_generic_lambdas - template struct _save_lambda { - ArrayType const& a; - h5::group g; - template void operator()(Is const&... is) const { h5_write(g, _h5_name(is...), a(is...)); } - }; - - template struct _load_lambda { - ArrayType& a; - h5::group g; - template void operator()(Is const&... is) { h5_read(g, _h5_name(is...), a(is...)); } - }; -#endif - } // details - - /* - * Write an array or a view into an hdf5 file when type is not fundamental - * ArrayType The type of the array/matrix/vector, etc.. - * g The h5 group - * name The name of the hdf5 array in the file/group where the stack will be stored - * A The array to be stored - * The HDF5 exceptions will be caught and rethrown as TRIQS_RUNTIME_ERROR (with a full stackstrace, cf triqs doc). - */ - template - std::c14::enable_if_t::value && !has_scalar_or_string_value_type::value> - h5_write(h5::group gr, std::string name, ArrayType const& a) { - if (a.is_empty()) TRIQS_RUNTIME_ERROR << " Can not save an empty array into hdf5"; - auto gr2 = gr.create_group(name); - gr2.write_triqs_hdf5_data_scheme(a); - // save the shape - array sha(ArrayType::rank); - for (int u = 0; u < ArrayType::rank; ++u) sha(u) = a.shape()[u]; - h5_write(gr2, "shape", sha); -#ifndef __cpp_generic_lambdas - foreach(a, h5_impl::_save_lambda{a, gr2}); -#else - foreach(a, [&](auto... is) { h5_write(gr2, h5_impl::_h5_name(is...), a(is...)); }); -#endif - } - - /* - * Read an array or a view from an hdf5 file when type is not fundamental - * ArrayType The type of the array/matrix/vector, etc.. - * g The h5 group - * name The name of the hdf5 array in the file/group where the stack will be stored - * A The array to be stored - * The HDF5 exceptions will be caught and rethrown as TRIQS_RUNTIME_ERROR (with a full stackstrace, cf triqs doc). - */ - template - std::c14::enable_if_t::value && !has_scalar_or_string_value_type::value> - h5_read(h5::group gr, std::string name, ArrayType& a) { - static_assert(!std::is_const::value, "Can not read in const object"); - auto gr2 = gr.open_group(name); - // TODO checking scheme... - // load the shape - auto sha2 = a.shape(); - array sha; - h5_read(gr2, "shape", sha); - if (first_dim(sha) != sha2.size()) - TRIQS_RUNTIME_ERROR << " array> load : rank mismatch. Expected " << sha2.size()<< " Got " << first_dim(sha); - for (int u = 0; u < sha2.size(); ++u) sha2[u] = sha(u); - if (a.shape() != sha2) a.resize(sha2); -#ifndef __cpp_generic_lambdas - foreach(a, h5_impl::_load_lambda{a, gr2}); -#else - foreach(a, [&](auto... is) { h5_read(gr2, h5_impl::_h5_name(is...), a(is...)); }); -#endif - } -} -} -#endif +}} diff --git a/triqs/arrays/python/numpy_extractor.cpp b/triqs/arrays/python/numpy_extractor.cpp new file mode 100644 index 00000000..70e3c46e --- /dev/null +++ b/triqs/arrays/python/numpy_extractor.cpp @@ -0,0 +1,122 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2011 by O. Parcollet + * + * TRIQS 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. + * + * TRIQS 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 + * TRIQS. If not, see . + * + ******************************************************************************/ +#include "../impl/common.hpp" +#include "../impl/traits.hpp" +#include "../impl/flags.hpp" +#include "../../utility/mini_vector.hpp" +#include "./numpy_extractor.hpp" +#ifdef TRIQS_WITH_PYTHON_SUPPORT + +namespace triqs { namespace arrays { namespace numpy_interface { + + + PyObject * numpy_extractor_impl ( PyObject * X, bool allow_copy, std::string type_name, + int elementsType, int rank, size_t * lengths, std::ptrdiff_t * strides, size_t size_of_ValueType) { + + PyObject * numpy_obj; + + if (X==NULL) TRIQS_RUNTIME_ERROR<<"numpy interface : the python object is NULL !"; + if (_import_array()!=0) TRIQS_RUNTIME_ERROR <<"Internal Error in importing numpy"; + + static const char * error_msg = " A deep copy of the object would be necessary while views are supposed to guarantee to present a *view* of the python data.\n"; + + if (!allow_copy) { + if (!PyArray_Check(X)) throw copy_exception () << error_msg<<" Indeed the object was not even an array !\n"; + if ( elementsType != PyArray_TYPE((PyArrayObject*)X)) + throw copy_exception () << error_msg<<" The deep copy is caused by a type mismatch of the elements. Expected "<< type_name<< " and found XXX \n"; + PyArrayObject *arr = (PyArrayObject *)X; +#ifdef TRIQS_NUMPY_VERSION_LT_17 + if ( arr->nd != rank) throw copy_exception () << error_msg<<" Rank mismatch . numpy array is of rank "<< arr->nd << "while you ask for rank "<< rank<<". \n"; +#else + if ( PyArray_NDIM(arr) != rank) throw copy_exception () << error_msg<<" Rank mismatch . numpy array is of rank "<< PyArray_NDIM(arr) << "while you ask for rank "<< rank<<". \n"; +#endif + numpy_obj = X; Py_INCREF(X); + } + else { + // From X, we ask the numpy library to make a numpy, and of the correct type. + // This handles automatically the cases where : + // - we have list, or list of list/tuple + // - the numpy type is not the one we want. + // - adjust the dimension if needed + // If X is an array : + // - if Order is same, don't change it + // - else impose it (may provoque a copy). + // if X is not array : + // - Order = FortranOrder or SameOrder - > Fortran order otherwise C + + //bool ForceCast = false;// Unless FORCECAST is present in flags, this call will generate an error if the data type cannot be safely obtained from the object. + int flags = 0; //(ForceCast ? NPY_FORCECAST : 0) ;// do NOT force a copy | (make_copy ? NPY_ENSURECOPY : 0); + if (!(PyArray_Check(X) )) + //flags |= ( IndexMapType::traversal_order == indexmaps::mem_layout::c_order(rank) ? NPY_C_CONTIGUOUS : NPY_F_CONTIGUOUS); //impose mem order +#ifdef TRIQS_NUMPY_VERSION_LT_17 + flags |= (NPY_C_CONTIGUOUS); //impose mem order +#else + flags |= (NPY_ARRAY_C_CONTIGUOUS); //impose mem order +#endif + numpy_obj= PyArray_FromAny(X,PyArray_DescrFromType(elementsType), rank,rank, flags , NULL ); + + // do several checks + if (!numpy_obj) {// The convertion of X to a numpy has failed ! + if (PyErr_Occurred()) { + //PyErr_Print(); + PyErr_Clear(); + } + TRIQS_RUNTIME_ERROR<<"numpy interface : the python object is not convertible to a numpy. "; + } + assert (PyArray_Check(numpy_obj)); assert((numpy_obj->ob_refcnt==1) || ((numpy_obj ==X))); + + PyArrayObject *arr_obj; + arr_obj = (PyArrayObject *)numpy_obj; + try { +#ifdef TRIQS_NUMPY_VERSION_LT_17 + if (arr_obj->nd!=rank) TRIQS_RUNTIME_ERROR<<"numpy interface : internal error : dimensions do not match"; + if (arr_obj->descr->type_num != elementsType) + TRIQS_RUNTIME_ERROR<<"numpy interface : internal error : incorrect type of element :" <descr->type_num <<" vs "<type_num != elementsType) + TRIQS_RUNTIME_ERROR<<"numpy interface : internal error : incorrect type of element :" <type_num <<" vs "<nd; // we know that dim == rank + for (size_t i=0; i< dim ; ++i) { + lengths[i] = size_t(arr_obj-> dimensions[i]); + strides[i] = std::ptrdiff_t(arr_obj-> strides[i])/ size_of_ValueType; + } +#else + const size_t dim = PyArray_NDIM(arr_obj); // we know that dim == rank + for (size_t i=0; i< dim ; ++i) { + lengths[i] = size_t( PyArray_DIMS(arr_obj)[i]); + strides[i] = std::ptrdiff_t( PyArray_STRIDES(arr_obj)[i])/ size_of_ValueType; + } +#endif + + return numpy_obj; + } +}}} +#endif diff --git a/triqs/arrays/python/numpy_extractor.hpp b/triqs/arrays/python/numpy_extractor.hpp index c1d2925b..32e1c7ae 100644 --- a/triqs/arrays/python/numpy_extractor.hpp +++ b/triqs/arrays/python/numpy_extractor.hpp @@ -21,13 +21,15 @@ #ifndef TRIQS_ARRAYS_NUMPY_EXTRACTOR_H #define TRIQS_ARRAYS_NUMPY_EXTRACTOR_H -#ifdef TRIQS_WITH_PYTHON_SUPPORT #include "../storages/shared_block.hpp" #include "triqs/utility/exceptions.hpp" +#ifdef TRIQS_WITH_PYTHON_SUPPORT #include "numpy/arrayobject.h" namespace triqs { namespace arrays { namespace numpy_interface { + using utility::mini_vector; + inline std::string object_to_string (PyObject * p) { if (!PyString_Check(p)) TRIQS_RUNTIME_ERROR<<" Internal error, expected a python string ....."; return PyString_AsString(p); @@ -58,98 +60,9 @@ namespace triqs { namespace arrays { namespace numpy_interface { struct copy_exception : public triqs::runtime_error {}; // return a NEW (owned) reference - // - inline PyObject * numpy_extractor_impl ( PyObject * X, bool allow_copy, std::string type_name, - int elementsType, int rank, size_t * lengths, std::ptrdiff_t * strides, size_t size_of_ValueType) { - - PyObject * numpy_obj; - - if (X==NULL) TRIQS_RUNTIME_ERROR<<"numpy interface : the python object is NULL !"; - if (_import_array()!=0) TRIQS_RUNTIME_ERROR <<"Internal Error in importing numpy"; - - static const char * error_msg = " A deep copy of the object would be necessary while views are supposed to guarantee to present a *view* of the python data.\n"; - - if (!allow_copy) { - if (!PyArray_Check(X)) throw copy_exception () << error_msg<<" Indeed the object was not even an array !\n"; - if ( elementsType != PyArray_TYPE((PyArrayObject*)X)) - throw copy_exception () << error_msg<<" The deep copy is caused by a type mismatch of the elements. Expected "<< type_name<< " and found XXX \n"; - PyArrayObject *arr = (PyArrayObject *)X; -#ifdef TRIQS_NUMPY_VERSION_LT_17 - if ( arr->nd != rank) throw copy_exception () << error_msg<<" Rank mismatch . numpy array is of rank "<< arr->nd << "while you ask for rank "<< rank<<". \n"; -#else - if ( PyArray_NDIM(arr) != rank) throw copy_exception () << error_msg<<" Rank mismatch . numpy array is of rank "<< PyArray_NDIM(arr) << "while you ask for rank "<< rank<<". \n"; -#endif - numpy_obj = X; Py_INCREF(X); - } - else { - // From X, we ask the numpy library to make a numpy, and of the correct type. - // This handles automatically the cases where : - // - we have list, or list of list/tuple - // - the numpy type is not the one we want. - // - adjust the dimension if needed - // If X is an array : - // - if Order is same, don't change it - // - else impose it (may provoque a copy). - // if X is not array : - // - Order = FortranOrder or SameOrder - > Fortran order otherwise C - - //bool ForceCast = false;// Unless FORCECAST is present in flags, this call will generate an error if the data type cannot be safely obtained from the object. - int flags = 0; //(ForceCast ? NPY_FORCECAST : 0) ;// do NOT force a copy | (make_copy ? NPY_ENSURECOPY : 0); - if (!(PyArray_Check(X) )) - //flags |= ( IndexMapType::traversal_order == indexmaps::mem_layout::c_order(rank) ? NPY_C_CONTIGUOUS : NPY_F_CONTIGUOUS); //impose mem order -#ifdef TRIQS_NUMPY_VERSION_LT_17 - flags |= (NPY_C_CONTIGUOUS); //impose mem order -#else - flags |= (NPY_ARRAY_C_CONTIGUOUS); //impose mem order -#endif - numpy_obj= PyArray_FromAny(X,PyArray_DescrFromType(elementsType), rank,rank, flags , NULL ); - - // do several checks - if (!numpy_obj) {// The convertion of X to a numpy has failed ! - if (PyErr_Occurred()) { - //PyErr_Print(); - PyErr_Clear(); - } - TRIQS_RUNTIME_ERROR<<"numpy interface : the python object is not convertible to a numpy. "; - } - assert (PyArray_Check(numpy_obj)); assert((numpy_obj->ob_refcnt==1) || ((numpy_obj ==X))); - - PyArrayObject *arr_obj; - arr_obj = (PyArrayObject *)numpy_obj; - try { -#ifdef TRIQS_NUMPY_VERSION_LT_17 - if (arr_obj->nd!=rank) TRIQS_RUNTIME_ERROR<<"numpy interface : internal error : dimensions do not match"; - if (arr_obj->descr->type_num != elementsType) - TRIQS_RUNTIME_ERROR<<"numpy interface : internal error : incorrect type of element :" <descr->type_num <<" vs "<type_num != elementsType) - TRIQS_RUNTIME_ERROR<<"numpy interface : internal error : incorrect type of element :" <type_num <<" vs "<nd; // we know that dim == rank - for (size_t i=0; i< dim ; ++i) { - lengths[i] = size_t(arr_obj-> dimensions[i]); - strides[i] = std::ptrdiff_t(arr_obj-> strides[i])/ size_of_ValueType; - } -#else - const size_t dim = PyArray_NDIM(arr_obj); // we know that dim == rank - for (size_t i=0; i< dim ; ++i) { - lengths[i] = size_t( PyArray_DIMS(arr_obj)[i]); - strides[i] = std::ptrdiff_t( PyArray_STRIDES(arr_obj)[i])/ size_of_ValueType; - } -#endif - - return numpy_obj; - } - + PyObject * numpy_extractor_impl ( PyObject * X, bool allow_copy, std::string type_name, + int elementsType, int rank, size_t * lengths, std::ptrdiff_t * strides, size_t size_of_ValueType); + // a little template class template struct numpy_extractor { diff --git a/triqs/gfs/meshes/matsubara_freq.hpp b/triqs/gfs/meshes/matsubara_freq.hpp index b10137f0..aecdf41e 100644 --- a/triqs/gfs/meshes/matsubara_freq.hpp +++ b/triqs/gfs/meshes/matsubara_freq.hpp @@ -129,7 +129,7 @@ namespace gfs { h5_write(gr, "max", ((Fermion ? 1 : 0) + 2 * m.size()) * M_PI / beta); h5_write(gr, "kind", 2); } else { // A strange way : to preserve backward compatibility for old archive. - h5_write(gr, "start_at_0", m._positive_only); + h5_write(gr, "start_at_0", (m._positive_only?1:0)); } } @@ -138,11 +138,11 @@ namespace gfs { h5::group gr = fg.open_group(subgroup_name); typename matsubara_freq_mesh::domain_t dom; int L; - bool s = true; + int s = 1; h5_read(gr, "domain", dom); h5_read(gr, "size", L); if (gr.has_key("start_at_0")) h5_read(gr, "start_at_0", s); - m = matsubara_freq_mesh{std::move(dom), L, s}; + m = matsubara_freq_mesh{std::move(dom), L, (s==1)}; } friend class boost::serialization::access; diff --git a/triqs/h5/base.cpp b/triqs/h5/base.cpp new file mode 100644 index 00000000..dda485b3 --- /dev/null +++ b/triqs/h5/base.cpp @@ -0,0 +1,96 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2011-2014 by O. Parcollet + * + * TRIQS 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. + * + * TRIQS 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 + * TRIQS. If not, see . + * + ******************************************************************************/ +#include "./base.hpp" +#include + +namespace triqs { + +namespace h5 { + + // dataspace from lengths and strides. Correct for the complex. strides must be >0 + H5::DataSpace dataspace_from_LS(int R, bool is_complex, hsize_t const *Ltot, hsize_t const *L, hsize_t const *S, + hsize_t const *offset) { + int rank = R + (is_complex ? 1 : 0); + hsize_t totdimsf[rank], dimsf[rank], stridesf[rank], offsetf[rank]; // dataset dimensions + for (size_t u = 0; u < R; ++u) { + offsetf[u] = (offset ? offset[u] : 0); + dimsf[u] = L[u]; + totdimsf[u] = Ltot[u]; + stridesf[u] = S[u]; + } + if (is_complex) { + offsetf[rank - 1] = 0; + dimsf[rank - 1] = 2; + totdimsf[rank - 1] = 2; + stridesf[rank - 1] = 1; + } + H5::DataSpace ds(rank, totdimsf); + ds.selectHyperslab(H5S_SELECT_SET, dimsf, offsetf, stridesf); + return ds; + } + + /****************** Write string attribute *********************************************/ + + void write_string_attribute(H5::H5Object const *obj, std::string name, std::string value) { + H5::DataSpace attr_dataspace = H5::DataSpace(H5S_SCALAR); + // Create new string datatype for attribute + H5::StrType strdatatype(H5::PredType::C_S1, value.size()); + // Set up write buffer for attribute + // const H5std_string strwritebuf (value); + // Create attribute and write to it + H5::Attribute myatt_in = obj->createAttribute(name.c_str(), strdatatype, attr_dataspace); + // myatt_in.write(strdatatype, strwritebuf); + myatt_in.write(strdatatype, (void *)(value.c_str())); + } + + + /****************** Read string attribute *********************************************/ + + /// Return the attribute name of obj, and "" if the attribute does not exist. + std::string read_string_attribute(H5::H5Object const *obj, std::string name) { + std::string value = ""; + H5::Attribute attr; + if (H5LTfind_attribute(obj->getId(), name.c_str()) == 0) return value; // not present + // can not find how to get the size with hl. Using full interface + // herr_t err2 = H5LTget_attribute_string(gr.getId(), x.c_str(), name.c_str() , &(buf.front()) ) ; + // value.append( &(buf.front()) ); + try { + attr = obj->openAttribute(name.c_str()); + } + catch (H5::AttributeIException) { + return value; + } + try { + H5::DataSpace dataspace = attr.getSpace(); + int rank = dataspace.getSimpleExtentNdims(); + if (rank != 0) TRIQS_RUNTIME_ERROR << "Reading a string attribute and got rank !=0"; + size_t size = attr.getStorageSize(); + H5::StrType strdatatype(H5::PredType::C_S1, size + 1); + std::vector buf(size + 1, 0x00); + attr.read(strdatatype, (void *)(&buf[0])); + value.append(&(buf.front())); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + return value; + } +} +} + diff --git a/triqs/h5/base.hpp b/triqs/h5/base.hpp index fbceff5a..3c12ae1d 100644 --- a/triqs/h5/base.hpp +++ b/triqs/h5/base.hpp @@ -2,7 +2,7 @@ * * TRIQS: a Toolbox for Research in Interacting Quantum Systems * - * Copyright (C) 2011-2013 by O. Parcollet + * Copyright (C) 2011-2014 by O. Parcollet * * TRIQS 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 @@ -18,12 +18,10 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_H5_BASE_H -#define TRIQS_H5_BASE_H +#pragma once #include #include #include -#include "hdf5_hl.h" #include namespace triqs { @@ -108,21 +106,9 @@ namespace triqs { get_data_ptr (S * p) {return p;} // dataspace from lengths and strides. Correct for the complex. strides must be >0 - template - H5::DataSpace dataspace_from_LS ( - mini_vector const & Ltot, - mini_vector const & L, - mini_vector const & S, - mini_vector const & offset = mini_vector() ) - { - static const unsigned int rank = R + (IsComplex ? 1 : 0); - hsize_t totdimsf[rank], dimsf [rank], stridesf[rank], offsetf[rank]; // dataset dimensions - for (size_t u=0; ucreateAttribute(name.c_str(), strdatatype, attr_dataspace); - //myatt_in.write(strdatatype, strwritebuf); - myatt_in.write(strdatatype, (void *)(value.c_str())); - } + void write_string_attribute(H5::H5Object const *obj, std::string name, std::string value); - - /****************** Read string attribute *********************************************/ - - /// Return the attribute name of obj, and "" if the attribute does not exist. - inline std::string read_string_attribute (H5::H5Object const * obj, std::string name ) { - std::string value =""; - H5::Attribute attr; - if (H5LTfind_attribute(obj -> getId(), name.c_str() )==0) return value;// not present - // can not find how to get the size with hl. Using full interface - //herr_t err2 = H5LTget_attribute_string(gr.getId(), x.c_str(), name.c_str() , &(buf.front()) ) ; - //value.append( &(buf.front()) ); - try { attr= obj->openAttribute(name.c_str());} - catch (H5::AttributeIException) { return value;} - try { - H5::DataSpace dataspace = attr.getSpace(); - int rank = dataspace.getSimpleExtentNdims(); - if (rank != 0) TRIQS_RUNTIME_ERROR << "Reading a string attribute and got rank !=0"; - size_t size = attr.getStorageSize(); - H5::StrType strdatatype(H5::PredType::C_S1, size+1); - std::vector buf(size+1, 0x00); - attr.read(strdatatype, (void *)(&buf[0])); - value.append( &(buf.front()) ); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - return value; - } + /// Returns the attribute name of obj, and "" if the attribute does not exist. + std::string read_string_attribute(H5::H5Object const *obj, std::string name); }} -#endif diff --git a/triqs/h5/group.cpp b/triqs/h5/group.cpp index 8f86ce00..93ba624b 100644 --- a/triqs/h5/group.cpp +++ b/triqs/h5/group.cpp @@ -1,36 +1,144 @@ #include "./group.hpp" +#include -namespace triqs { namespace h5 { +namespace triqs { +namespace h5 { + + group::group(H5::Group g, H5::Group parent, std::string name_in_parent) + : _g(g), _parent(parent), _name_in_parent(name_in_parent) {} + + //group::group(H5::Group g) : _g(g) {} + + group::group(H5::H5File f) : _g(f.openGroup("/")) {} // can not fail, right ? + + group::group(std::string const &filename, int flag) { + H5::H5File f(filename, H5F_ACC_TRUNC); + *this = group(f); + } + + group::group(hid_t id_, bool is_group) { + if (is_group) { + _g.setId(id_); + } else { + H5::H5File f; + f.setId(id_); + *this = group(f); + } + } + + void group::_write_triqs_hdf5_data_scheme(const char *a) { h5::write_string_attribute(&_g, "TRIQS_HDF5_data_scheme", a); } + + std::string group::read_triqs_hdf5_data_scheme() const { return read_string_attribute(&_g, "TRIQS_HDF5_data_scheme"); } + + bool group::has_key(std::string const &key) const { return (H5Lexists(_g.getId(), key.c_str(), H5P_DEFAULT)); } + + void group::unlink_key_if_exists(std::string const &key) const { + if (this->has_key(key)) _g.unlink(key.c_str()); + } + + group group::open_group(std::string const &key) const { + if (!has_key(key)) TRIQS_RUNTIME_ERROR << "no subgroup " << key << " in the group"; + group res; + try { + res = group(_g.openGroup(key.c_str()), _g, key); + } // has a parent + catch (H5::GroupIException const &e) { + TRIQS_RUNTIME_ERROR << "Error in opening the subgroup " << key << "\n H5 error message : \n " << e.getCDetailMsg(); + } + return res; + } + + + /// Open an existing DataSet. Throw it if does not exists + H5::DataSet group::open_dataset(std::string const &key) const { + if (!has_key(key)) TRIQS_RUNTIME_ERROR << "no dataset " << key << " in the group"; + H5::DataSet res; + try { + res = _g.openDataSet(key.c_str()); + } + catch (H5::GroupIException const &e) { + TRIQS_RUNTIME_ERROR << "Error in opening the dataset " << key << "\n H5 error message : \n " << e.getCDetailMsg(); + } + return res; + } + /** + * \brief Create a subgroup. + * \param key The name of the subgroup + * \param delete_if_exists Unlink the group if it exists + */ + group group::create_group(std::string const &key, bool delete_if_exists) const { + unlink_key_if_exists(key); + return group(_g.createGroup(key.c_str()), _g, key); + } + /** + * \brief Create a dataset. + * \param key The name of the subgroup + * \param args Other parameters are forwarded to H5::Group + * + * NB : It unlinks the dataset if it exists. + */ + + H5::DataSet group::create_dataset(std::string const &key, const H5::DataType &data_type, const H5::DataSpace &data_space, + const H5::DSetCreatPropList &create_plist) const { + unlink_key_if_exists(key); + H5::DataSet res; + try { + res = _g.createDataSet(key.c_str(), data_type, data_space, create_plist); + } + catch (H5::GroupIException const &e) { + TRIQS_RUNTIME_ERROR << "Error in creating the dataset " << key << "\n H5 error message : \n " << e.getCDetailMsg(); + } + return res; + } + + void group::write_string_attribute (std::string const & obj_name, std::string const & attr_name, std::string const & value){ + herr_t err = H5LTset_attribute_string(_g.getId(),obj_name.c_str(),attr_name.c_str(), value.c_str() ) ; + if (err<0) TRIQS_RUNTIME_ERROR << "Error in setting attribute of "<< obj_name<<" named "<< attr_name << " to " << value; + } + + //----------------------------------------------------------------------- extern "C" { - herr_t get_group_elements_name_ds(hid_t loc_id, const char *name, void *opdata) { - H5O_info_t object_info; herr_t err =H5Oget_info_by_name(loc_id,name,&object_info,H5P_DEFAULT ); - if (err<0) TRIQS_RUNTIME_ERROR << "get_group_elements_name_ds internal"; - if (object_info.type == H5O_TYPE_DATASET) static_cast *>(opdata)->push_back(name); - return 0; - } - herr_t get_group_elements_name_grp(hid_t loc_id, const char *name, void *opdata) { - H5O_info_t object_info; herr_t err =H5Oget_info_by_name(loc_id,name,&object_info,H5P_DEFAULT ); - if (err<0) TRIQS_RUNTIME_ERROR << "get_group_elements_name_grp internal"; - if (object_info.type == H5O_TYPE_GROUP) static_cast *>(opdata)->push_back(name); - return 0; - } + herr_t get_group_elements_name_ds(hid_t loc_id, const char *name, void *opdata) { + H5O_info_t object_info; + herr_t err = H5Oget_info_by_name(loc_id, name, &object_info, H5P_DEFAULT); + if (err < 0) TRIQS_RUNTIME_ERROR << "get_group_elements_name_ds internal"; + if (object_info.type == H5O_TYPE_DATASET) static_cast *>(opdata)->push_back(name); + return 0; + } + herr_t get_group_elements_name_grp(hid_t loc_id, const char *name, void *opdata) { + H5O_info_t object_info; + herr_t err = H5Oget_info_by_name(loc_id, name, &object_info, H5P_DEFAULT); + if (err < 0) TRIQS_RUNTIME_ERROR << "get_group_elements_name_grp internal"; + if (object_info.type == H5O_TYPE_GROUP) static_cast *>(opdata)->push_back(name); + return 0; + } } //----------------------------------------------------------------------- - std::vector group::get_all_subgroup_names(std::string const & key) const { + std::vector group::get_all_subgroup_names(std::string const &key) const { std::vector grp_name; H5::Group F(_g); - F.iterateElems(key.c_str(), NULL, get_group_elements_name_grp, static_cast(&grp_name)); + F.iterateElems(key.c_str(), NULL, get_group_elements_name_grp, static_cast(&grp_name)); return grp_name; } - std::vector group::get_all_dataset_names(std::string const & key) const { + std::vector group::get_all_dataset_names(std::string const &key) const { std::vector ds_name; H5::Group F(_g); - F.iterateElems(key.c_str(), NULL, get_group_elements_name_ds, static_cast(&ds_name)); + F.iterateElems(key.c_str(), NULL, get_group_elements_name_ds, static_cast(&ds_name)); return ds_name; } -}} + std::vector group::get_all_subgroup_names() const { + if (!has_parent()) TRIQS_RUNTIME_ERROR << "Group hdf5 : parent not found"; + return group(_parent).get_all_subgroup_names(_name_in_parent); + } + + std::vector group::get_all_dataset_names() const { + if (!has_parent()) TRIQS_RUNTIME_ERROR << "Group hdf5 : parent not found"; + return group(_parent).get_all_dataset_names(_name_in_parent); + } +} +} diff --git a/triqs/h5/group.hpp b/triqs/h5/group.hpp index d26a63f6..95c5cb5d 100644 --- a/triqs/h5/group.hpp +++ b/triqs/h5/group.hpp @@ -18,107 +18,84 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_H5_GROUP_H -#define TRIQS_H5_GROUP_H +#pragma once #include "./base.hpp" -namespace triqs { namespace h5 { + +namespace triqs { +namespace h5 { + + // using hid_t = int; + /** * \brief A local derivative of Group. - * Rational : use ADL for h5_read/h5_write, catch and rethrow exception, add some policy for opening/creating + * Rational : use ADL for h5_read/h5_write, catch and rethrow exception, add some policy for opening/creating */ class group { - H5::Group _g, _parent; - std::string _name_in_parent; - group(H5::Group g, H5::Group parent, std::string name_in_parent) : _g(g), _parent(parent), _name_in_parent(name_in_parent) {} + // hid_t _g_id, _parent_id; + H5::Group _g, _parent; + std::string _name_in_parent; + group(H5::Group g, H5::Group parent, std::string name_in_parent); + void _write_triqs_hdf5_data_scheme(const char *a); + public: - group() = default; - group(group const &) = default; - group(H5::Group g) : _g(g) {} - - /// Takes the "/" group at the top of the file. - group (H5::H5File f) : _g(f.openGroup("/")) {} // can not fail, right ? + group() = default; + group(group const &) = default; + group(H5::Group g) : _g(g) {} + group(H5::H5File f); /// Takes the "/" group at the top of the file + group(std::string const &filename, int flag); + group(hid_t id_, bool is_group); - group(std::string const & filename, int flag) { H5::H5File f(filename, H5F_ACC_TRUNC); *this = group(f);} + bool has_parent() const { return _name_in_parent != ""; } - group(hid_t id_, bool is_group) { - if (is_group) { _g.setId(id_); } - else { H5::H5File f; f.setId(id_); *this = group(f); } - } + /// Write the triqs tag of the group if it is an object. + template void write_triqs_hdf5_data_scheme(T const &obj) { + _write_triqs_hdf5_data_scheme(get_triqs_hdf5_data_scheme(obj).c_str()); + } - bool has_parent() const { return _name_in_parent != "";} + /// Read the triqs tag of the group if it is an object. Returns"" if attribute is not present + std::string read_triqs_hdf5_data_scheme() const; - /// Write the triqs tag of the group if it is an object. - template void write_triqs_hdf5_data_scheme (T const & obj) { - h5::write_string_attribute( &_g, "TRIQS_HDF5_data_scheme" , get_triqs_hdf5_data_scheme(obj).c_str() ) ; - } - /// Read the triqs tag of the group if it is an object. Returns"" if attribute is not present - std::string read_triqs_hdf5_data_scheme() const { return read_string_attribute(&_g, "TRIQS_HDF5_data_scheme") ; } - /// - bool has_key (std::string const & key) const { return (H5Lexists(_g.getId(), key.c_str(), H5P_DEFAULT)); } - /// - void unlink_key_if_exists(std::string const & key) const { if (this->has_key(key)) _g.unlink(key.c_str()); } - /// Open a subgroup. Throw it if does not exists - group open_group(std::string const & key) const { - if (!has_key(key)) TRIQS_RUNTIME_ERROR << "no subgroup "< - H5::DataSet create_dataset(std::string const & key, Args && ... args) const { - unlink_key_if_exists(key); - H5::DataSet res; - try{ res = _g.createDataSet(key.c_str(), std::forward(args)...);} - catch (H5::GroupIException const & e){ TRIQS_RUNTIME_ERROR << "Error in creating the dataset "<< key <<"\n H5 error message : \n "<< e.getCDetailMsg(); } - return res; - } + /// + bool has_key(std::string const &key) const; + /// - /// Returns all names of subgroup of key in G - std::vector get_all_subgroup_names(std::string const & key) const; + void unlink_key_if_exists(std::string const &key) const; - std::vector get_all_subgroup_names() const { - if (!has_parent()) TRIQS_RUNTIME_ERROR << "Group hdf5 : parent not found"; - return group(_parent).get_all_subgroup_names(_name_in_parent); - } + /// Open a subgroup. Throw it if does not exists + group open_group(std::string const &key) const; - /// Returns all names of dataset of key in G - std::vector get_all_dataset_names(std::string const & key) const; + /// Open an existing DataSet. Throw it if does not exists + H5::DataSet open_dataset(std::string const &key) const; - std::vector get_all_dataset_names() const { - if (!has_parent()) TRIQS_RUNTIME_ERROR << "Group hdf5 : parent not found"; - return group(_parent).get_all_dataset_names(_name_in_parent); - } + /** + * \brief Create a subgroup. + * \param key The name of the subgroup + * \param delete_if_exists Unlink the group if it exists + */ + group create_group(std::string const &key, bool delete_if_exists = true) const; + + /** + * \brief Create a dataset. + * \param key The name of the subgroup + * + * NB : It unlinks the dataset if it exists. + */ + H5::DataSet create_dataset(std::string const &key, const H5::DataType &data_type, const H5::DataSpace &data_space, + const H5::DSetCreatPropList &create_plist= H5::DSetCreatPropList::DEFAULT) const; - void write_string_attribute (std::string const & obj_name, std::string const & attr_name, std::string const & value){ - herr_t err = H5LTset_attribute_string(_g.getId(),obj_name.c_str(),attr_name.c_str(), value.c_str() ) ; - if (err<0) TRIQS_RUNTIME_ERROR << "Error in setting attribute of "<< obj_name<<" named "<< attr_name << " to " << value; - } + /// Returns all names of subgroup of key in G + std::vector get_all_subgroup_names(std::string const &key) const; + /// Returns all names of subgroup of G + std::vector get_all_subgroup_names() const; + + /// Returns all names of dataset of key in G + std::vector get_all_dataset_names(std::string const &key) const; + + /// Returns all names of dataset of G + std::vector get_all_dataset_names() const; + + void write_string_attribute (std::string const & obj_name, std::string const & attr_name, std::string const & value); }; -}} -#endif +} +} diff --git a/triqs/h5/scalar.cpp b/triqs/h5/scalar.cpp new file mode 100644 index 00000000..736153be --- /dev/null +++ b/triqs/h5/scalar.cpp @@ -0,0 +1,74 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2011-2014 by O. Parcollet + * + * TRIQS 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. + * + * TRIQS 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 + * TRIQS. If not, see . + * + ******************************************************************************/ +#include "./scalar.hpp" +namespace triqs { namespace h5 { + + namespace { + + // TODO : Fix complex number .... + + // S must be scalar + template + void h5_write_impl (group g, std::string const & name, S const & A) { + try { + g.unlink_key_if_exists(name); + H5::DataSet ds = g.create_dataset( name, data_type_file(), H5::DataSpace() ); + ds.write( (void *)(&A), data_type_memory(), H5::DataSpace() ); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + + template + void h5_read_impl (group g, std::string const & name, S & A) { + try { + H5::DataSet ds = g.open_dataset(name); + H5::DataSpace dataspace = ds.getSpace(); + int rank = dataspace.getSimpleExtentNdims(); + if (rank != 0) TRIQS_RUNTIME_ERROR << "triqs::array::h5::read. Rank mismatch : expecting a scalar (rank =0)" + <<" while the array stored in the hdf5 file has rank = "<(), H5::DataSpace() , H5::DataSpace() ); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + } + + void h5_write(group g, std::string const &name, int const &x) { h5_write_impl(g, name, x); } + void h5_write(group g, std::string const &name, long const &x) { h5_write_impl(g, name, x); } + void h5_write(group g, std::string const &name, size_t const &x) { h5_write_impl(g, name, x); } + void h5_write(group g, std::string const &name, char const &x) { h5_write_impl(g, name, x); } + void h5_write(group g, std::string const &name, bool const &x) { h5_write_impl(g, name, x); } + void h5_write(group g, std::string const &name, double const &x) { h5_write_impl(g, name, x); } + void h5_write(group g, std::string const &name, std::complex const &x) { h5_write_impl(g, name, x); } + + + void h5_read(group g, std::string const &name, int &x) { h5_read_impl(g, name, x); } + void h5_read(group g, std::string const &name, long &x) { h5_read_impl(g, name, x); } + void h5_read(group g, std::string const &name, size_t &x) { h5_read_impl(g, name, x); } + void h5_read(group g, std::string const &name, char &x) { h5_read_impl(g, name, x); } + void h5_read(group g, std::string const &name, bool &x) { h5_read_impl(g, name, x); } + void h5_read(group g, std::string const &name, double &x) { h5_read_impl(g, name, x); } + void h5_read(group g, std::string const &name, std::complex &x) { h5_read_impl(g, name, x); } + + +}} + + + diff --git a/triqs/h5/scalar.hpp b/triqs/h5/scalar.hpp index c83cff99..1b4c7bf6 100644 --- a/triqs/h5/scalar.hpp +++ b/triqs/h5/scalar.hpp @@ -18,39 +18,29 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_H5_SCALAR_H -#define TRIQS_H5_SCALAR_H +#pragma once #include "./group.hpp" namespace triqs { namespace h5 { - template - typename std::enable_if::value>::type - h5_write (group g, std::string const & name, S const & A) { - try { - g.unlink_key_if_exists(name); - H5::DataSet ds = g.create_dataset( name, data_type_file(), H5::DataSpace() ); - ds.write( (void *)(&A), data_type_memory(), H5::DataSpace() ); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } - - template - typename std::enable_if::value>::type - h5_read (group g, std::string const & name, S & A) { - try { - H5::DataSet ds = g.open_dataset(name); - H5::DataSpace dataspace = ds.getSpace(); - int rank = dataspace.getSimpleExtentNdims(); - if (rank != 0) TRIQS_RUNTIME_ERROR << "triqs::array::h5::read. Rank mismatch : expecting a scalar (rank =0)" - <<" while the array stored in the hdf5 file has rank = "<(), H5::DataSpace() , H5::DataSpace() ); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } - - // the complex number is missing here... -}} -#endif + // Issue several types are *implicitly* convertible to bool + // it could be confusing. Better to use int in hdf5 files. + void h5_write(group g, std::string const &name, int const &x); + void h5_write(group g, std::string const &name, long const &x); + void h5_write(group g, std::string const &name, size_t const &x); + void h5_write(group g, std::string const &name, bool const &x); + void h5_write(group g, std::string const &name, char const &x); + void h5_write(group g, std::string const &name, double const &x); + void h5_write(group g, std::string const &name, std::complex const &x); + void h5_read(group g, std::string const &name, int &x); + void h5_read(group g, std::string const &name, long &x); + void h5_read(group g, std::string const &name, size_t &x); + void h5_read(group g, std::string const &name, bool &x); + void h5_read(group g, std::string const &name, char &x); + void h5_read(group g, std::string const &name, double &x); + void h5_read(group g, std::string const &name, std::complex &x); + // the implementation complex number is missing ... +} +} diff --git a/triqs/h5/string.cpp b/triqs/h5/string.cpp new file mode 100644 index 00000000..dc0a03cb --- /dev/null +++ b/triqs/h5/string.cpp @@ -0,0 +1,51 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2011-2014 by O. Parcollet + * + * TRIQS 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. + * + * TRIQS 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 + * TRIQS. If not, see . + * + ******************************************************************************/ +#include "./string.hpp" +namespace triqs { + + namespace h5 { + + void h5_write (group g, std::string const & name, std::string const & value) { + try { + H5::StrType strdatatype(H5::PredType::C_S1, value.size()); + H5::DataSet ds = g.create_dataset(name, strdatatype, H5::DataSpace()); + ds.write((void*)(value.c_str()), strdatatype ); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + + void h5_read (group g, std::string const & name, std::string & value) { + try { + H5::DataSet ds = g.open_dataset(name); + H5::DataSpace dataspace = ds.getSpace(); + int rank = dataspace.getSimpleExtentNdims(); + if (rank != 0) TRIQS_RUNTIME_ERROR << "Reading a string and got rank !=0"; + size_t size = ds.getStorageSize(); + H5::StrType strdatatype(H5::PredType::C_S1, size); + std::vector buf(size+1, 0x00); + ds.read( (void *)(&buf[0]), strdatatype, H5::DataSpace(), H5::DataSpace() ); + value = ""; value.append( &(buf.front()) ); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + +}} + diff --git a/triqs/h5/string.hpp b/triqs/h5/string.hpp index 8e91ca7d..10007349 100644 --- a/triqs/h5/string.hpp +++ b/triqs/h5/string.hpp @@ -18,10 +18,10 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_H5_STRING_H -#define TRIQS_H5_STRING_H +#pragma once #include "./group.hpp" #include + namespace triqs { inline std::string get_triqs_hdf5_data_scheme(std::string const & ) { return "string";} @@ -35,14 +35,9 @@ namespace triqs { * \param value The string * \exception The HDF5 exceptions will be caught and rethrown as TRIQS_RUNTIME_ERROR (with a full stackstrace, cf triqs doc). */ - inline void h5_write (group g, std::string const & name, std::string const & value) { - try { - H5::StrType strdatatype(H5::PredType::C_S1, value.size()); - H5::DataSet ds = g.create_dataset(name, strdatatype, H5::DataSpace()); - ds.write((void*)(value.c_str()), strdatatype ); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } + void h5_write (group g, std::string const & name, std::string const & value); + + inline void h5_write (group g, std::string const & name, const char * s) { h5_write(g,name,std::string{s});} /** * \brief Read a string from an hdf5 file @@ -51,74 +46,9 @@ namespace triqs { * \param value The string to fill * \exception The HDF5 exceptions will be caught and rethrown as TRIQS_RUNTIME_ERROR (with a full stackstrace, cf triqs doc). */ - inline void h5_read (group g, std::string const & name, std::string & value) { - try { - H5::DataSet ds = g.open_dataset(name); - H5::DataSpace dataspace = ds.getSpace(); - int rank = dataspace.getSimpleExtentNdims(); - if (rank != 0) TRIQS_RUNTIME_ERROR << "Reading a string and got rank !=0"; - size_t size = ds.getStorageSize(); - H5::StrType strdatatype(H5::PredType::C_S1, size); - std::vector buf(size+1, 0x00); - ds.read( (void *)(&buf[0]), strdatatype, H5::DataSpace(), H5::DataSpace() ); - value = ""; value.append( &(buf.front()) ); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } - - - /*********************************** 1d array and std::vector of string *********************************/ - namespace detail { - - template - void write_1darray_vector_of_string_impl (group g, std::string const & name, ArrayVectorOfStringType const & V) { - size_t s=0; for (auto & x : V) s = std::max(s,x.size()); - try { - H5::DataSet ds; - H5::StrType strdatatype(H5::PredType::C_S1, s); - const size_t n = V.size(); - std::vector buf(n*(s+1), 0x00); - size_t i=0; for (auto & x : V) {strcpy( &buf[i*s], x.c_str()); ++i;} - - mini_vector L; L[0]=V.size(); - mini_vector S; S[0]=1; - auto d_space = dataspace_from_LS<1,false > (L,L,S); - - ds = g.create_dataset(name, strdatatype, d_space ); - ds.write( (void *)(&buf[0]),strdatatype, d_space ); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } - - //------------------------------------------------------------ - - template - void read_1darray_vector_of_string_impl (group g, std::string const & name, ArrayVectorOfStringType & V) { - try { - H5::DataSet ds = g.open_dataset(name); - H5::DataSpace dataspace = ds.getSpace(); - mini_vector dims_out; - int ndims = dataspace.getSimpleExtentDims( &dims_out[0], NULL); - if (ndims !=1) TRIQS_RUNTIME_ERROR << "triqs::h5 : Trying to read 1d array/vector . Rank mismatch : the array stored in the hdf5 file has rank = "< buf(Len*(size+1), 0x00); - - mini_vector L; L[0]=V.size(); - mini_vector S; S[0]=1; - auto d_space = dataspace_from_LS<1,false > (L,L,S); - - ds.read( (void *)(&buf[0]),strdatatype, d_space ); - size_t i=0; for (auto & x : V) { x = ""; x.append(&buf[i*(size)]); ++i;} - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } - }// detail + void h5_read (group g, std::string const & name, std::string & value); + inline void h5_read (group g, std::string const & name, char * s) = delete; + }} -#endif diff --git a/triqs/h5/vector.cpp b/triqs/h5/vector.cpp new file mode 100644 index 00000000..0a34979b --- /dev/null +++ b/triqs/h5/vector.cpp @@ -0,0 +1,120 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 2011-2014 by O. Parcollet + * + * TRIQS 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. + * + * TRIQS 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 + * TRIQS. If not, see . + * + ******************************************************************************/ +#include "./vector.hpp" + +namespace triqs { + namespace h5 { + + void h5_write (group g, std::string const & name, std::vector const & V) { + size_t s=0; for (auto & x : V) s = std::max(s,x.size()); + try { + H5::DataSet ds; + H5::StrType strdatatype(H5::PredType::C_S1, s); + const size_t n = V.size(); + std::vector buf(n*(s+1), 0x00); + size_t i=0; for (auto & x : V) {strcpy( &buf[i*s], x.c_str()); ++i;} + + hsize_t L[1], S[1]; + L[0] = V.size(); + S[0] = 1; + auto d_space = dataspace_from_LS(1,false,L,L,S); + + ds = g.create_dataset(name, strdatatype, d_space ); + ds.write( (void *)(&buf[0]),strdatatype, d_space ); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + + void h5_read (group g, std::string const & name, std::vector & V) { + try { + H5::DataSet ds = g.open_dataset(name); + H5::DataSpace dataspace = ds.getSpace(); + mini_vector dims_out; + int ndims = dataspace.getSimpleExtentDims( &dims_out[0], NULL); + if (ndims !=1) TRIQS_RUNTIME_ERROR << "triqs::h5 : Trying to read 1d array/vector . Rank mismatch : the array stored in the hdf5 file has rank = "< buf(Len*(size+1), 0x00); + + hsize_t L[1], S[1]; + L[0] = V.size(); + S[0] = 1; + auto d_space = dataspace_from_LS(1,false, L,L,S); + + ds.read( (void *)(&buf[0]),strdatatype, d_space ); + size_t i=0; for (auto & x : V) { x = ""; x.append(&buf[i*(size)]); ++i;} + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + + // implementation for vector of double and complex + namespace { + + // the dataspace corresponding to the array. Contiguous data only... + template H5::DataSpace data_space_for_vector(std::vector const &V) { + hsize_t L[1], S[1]; + S[0] = 1; + L[0] = V.size(); + return h5::dataspace_from_LS(1, triqs::is_complex::value, L, L, S); + } + + template + inline void h5_write_vector_impl (group g, std::string const & name, std::vector const & V) { + try { + H5::DataSet ds = g.create_dataset(name, h5::data_type_file(), data_space_for_vector(V) ); + ds.write( &V[0], h5::data_type_memory(), data_space_for_vector(V) ); + // if complex, to be python compatible, we add the __complex__ attribute + if (triqs::is_complex::value) h5::write_string_attribute(&ds,"__complex__","1"); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + + template + inline void h5_read_impl (group g, std::string const & name, std::vector & V) { + try { + H5::DataSet ds = g.open_dataset(name); + H5::DataSpace dataspace = ds.getSpace(); + static const unsigned int Rank = 1 + (triqs::is_complex::value ? 1 : 0); + int rank = dataspace.getSimpleExtentNdims(); + if (rank != Rank) TRIQS_RUNTIME_ERROR << "triqs : h5 : read vector. Rank mismatch : the array stored in the hdf5 file has rank = "< dims_out; + dataspace.getSimpleExtentDims( &dims_out[0], NULL); + V.resize(dims_out[0]); + ds.read( &V[0], h5::data_type_memory(), data_space_for_vector(V) , dataspace ); + } + TRIQS_ARRAYS_H5_CATCH_EXCEPTION; + } + }// impl namespace + + + void h5_write (group f, std::string const & name, std::vector const & V) { h5_write_vector_impl(f,name,V);} + void h5_write (group f, std::string const & name, std::vector> const & V) { h5_write_vector_impl(f,name,V);} + void h5_read (group f, std::string const & name, std::vector & V) { h5_read_impl(f,name,V);} + void h5_read (group f, std::string const & name, std::vector> & V) { h5_read_impl(f,name,V);} + + } +} + + diff --git a/triqs/h5/vector.hpp b/triqs/h5/vector.hpp index 6b951714..224ab1c3 100644 --- a/triqs/h5/vector.hpp +++ b/triqs/h5/vector.hpp @@ -18,8 +18,7 @@ * TRIQS. If not, see . * ******************************************************************************/ -#ifndef TRIQS_H5_VECTOR_H -#define TRIQS_H5_VECTOR_H +#pragma once #include "./group.hpp" #include "./string.hpp" #include @@ -28,72 +27,23 @@ namespace triqs { inline std::string get_triqs_hdf5_data_scheme(std::vector const & ) { return "vector";} - template - std::string get_triqs_hdf5_data_scheme(std::vector const&) { - using triqs::get_triqs_hdf5_data_scheme;// for the basic types, not found by ADL - std::stringstream fs; - fs<<"std::vector<"<"; - return fs.str(); - } + template std::string get_triqs_hdf5_data_scheme(std::vector const&) { + using triqs::get_triqs_hdf5_data_scheme; // for the basic types, not found by ADL + return "std::vector<" + get_triqs_hdf5_data_scheme(T()) + ">"; + } namespace h5 { - // special case of vector of string - inline void h5_write (group f, std::string const & name, std::vector const & V) { - detail::write_1darray_vector_of_string_impl(f,name,V); - } + void h5_write (group f, std::string const & name, std::vector const & V); + void h5_read (group f, std::string const & name, std::vector & V); - inline void h5_read (group f, std::string const & name, std::vector & V) { - detail::read_1darray_vector_of_string_impl(f,name,V); - } + void h5_write (group f, std::string const & name, std::vector const & V); + void h5_write (group f, std::string const & name, std::vector> const & V); - // implementation for vector of double and complex - namespace vector_impl { - - // the dataspace corresponding to the array. Contiguous data only... - template - H5::DataSpace data_space_for_vector (std::vector const & V) { - mini_vector L,S; S[0]=1;L[0] = V.size(); - return h5::dataspace_from_LS<1, triqs::is_complex::value > (L,L,S); - } - - template - inline void h5_write_vector_impl (group g, std::string const & name, std::vector const & V) { - try { - H5::DataSet ds = g.create_dataset(name, h5::data_type_file(), data_space_for_vector(V) ); - ds.write( &V[0], h5::data_type_memory(), data_space_for_vector(V) ); - // if complex, to be python compatible, we add the __complex__ attribute - if (triqs::is_complex::value) h5::write_string_attribute(&ds,"__complex__","1"); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } - - template - inline void h5_read_impl (group g, std::string const & name, std::vector & V) { - try { - H5::DataSet ds = g.open_dataset(name); - H5::DataSpace dataspace = ds.getSpace(); - static const unsigned int Rank = 1 + (triqs::is_complex::value ? 1 : 0); - int rank = dataspace.getSimpleExtentNdims(); - if (rank != Rank) TRIQS_RUNTIME_ERROR << "triqs : h5 : read vector. Rank mismatch : the array stored in the hdf5 file has rank = "< dims_out; - dataspace.getSimpleExtentDims( &dims_out[0], NULL); - V.resize(dims_out[0]); - ds.read( &V[0], h5::data_type_memory(), data_space_for_vector(V) , dataspace ); - } - TRIQS_ARRAYS_H5_CATCH_EXCEPTION; - } - }// impl namespace - - - inline void h5_write (group f, std::string const & name, std::vector const & V) { vector_impl::h5_write_vector_impl(f,name,V);} - inline void h5_write (group f, std::string const & name, std::vector> const & V) { vector_impl::h5_write_vector_impl(f,name,V);} - - inline void h5_read (group f, std::string const & name, std::vector & V) { vector_impl::h5_read_impl(f,name,V);} - inline void h5_read (group f, std::string const & name, std::vector> & V) { vector_impl::h5_read_impl(f,name,V);} + void h5_read (group f, std::string const & name, std::vector & V); + void h5_read (group f, std::string const & name, std::vector> & V); } } -#endif