From 1c9d6dacfa5f16d86c1623ac59e78d4fa62cedad Mon Sep 17 00:00:00 2001 From: Olivier Parcollet Date: Thu, 13 Feb 2014 17:28:55 +0100 Subject: [PATCH] arrays: h5 read/write for arrays of complex types - array of complex type (not fundamental) can now be saved/loaded to h5 - with a test with array> --- test/triqs/gfs/array_gf.cpp | 40 +++++++++++++ test/triqs/parameters/param.cpp | 3 +- triqs/arrays/h5/simple_read_write.hpp | 83 ++++++++++++++++++++++++++- triqs/parameters/opaque_object_h5.cpp | 3 +- triqs/utility/mini_vector.hpp | 4 +- 5 files changed, 129 insertions(+), 4 deletions(-) create mode 100644 test/triqs/gfs/array_gf.cpp diff --git a/test/triqs/gfs/array_gf.cpp b/test/triqs/gfs/array_gf.cpp new file mode 100644 index 00000000..ddff66d9 --- /dev/null +++ b/test/triqs/gfs/array_gf.cpp @@ -0,0 +1,40 @@ +#include +#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl; + +using namespace triqs::gfs; +using namespace triqs::arrays; + +int main(int argc, char* argv[]) { + + try { + + triqs::clef::placeholder<0> w_; + auto agf = array, 2>{2, 3}; + auto bgf = array, 2>{2, 3}; + agf() = gf{{10.0, Fermion}, {1, 1}}; + agf(0, 0)(w_) << 1 / (w_ + 2); + + array A(2,2); A()=0;A(0,0) = 1.3; A(1,1) = -8.2; + array, 1> aa(2), bb; + aa(0) = A; + aa(1) = 2 * A; + bb = aa; + + { + H5::H5File file("ess_array_gf.h5", H5F_ACC_TRUNC); + h5_write(file, "Agf", agf); + h5_write(file, "aa", aa); + } + { + H5::H5File file("ess_array_gf.h5", H5F_ACC_RDONLY); + h5_read(file, "Agf", bgf); + h5_read(file, "aa", bb); + } + { + H5::H5File file("ess_array_gf2.h5", H5F_ACC_TRUNC); + h5_write(file, "Agf", bgf); + h5_write(file, "aa", bb); + } + } + TRIQS_CATCH_AND_ABORT; +} diff --git a/test/triqs/parameters/param.cpp b/test/triqs/parameters/param.cpp index 1fe7ff8c..32b81904 100644 --- a/test/triqs/parameters/param.cpp +++ b/test/triqs/parameters/param.cpp @@ -92,7 +92,7 @@ int main() { C = extract(P["B"]); std::cout << "C" << C << std::endl; - array, 1> aa(3); + array, 1> aa(2); aa(0) = A; aa(1) = 2*A; P["aa"] = aa; @@ -143,6 +143,7 @@ int main() { parameters P4; parameters::register_type>(); + //parameters::register_type>(); std::cout << "P4 before : "<< P4<< std::endl ; { H5::H5File file( "ess.h5", H5F_ACC_RDONLY ); diff --git a/triqs/arrays/h5/simple_read_write.hpp b/triqs/arrays/h5/simple_read_write.hpp index f55ef991..d897c85c 100644 --- a/triqs/arrays/h5/simple_read_write.hpp +++ b/triqs/arrays/h5/simple_read_write.hpp @@ -147,7 +147,7 @@ namespace arrays { V = res; } - } // namespace h5impl + } // namespace h5_impl // a trait to detect if A::value_type exists and is a scalar or a string // used to exclude array> @@ -192,8 +192,89 @@ namespace arrays { template ENABLE_IFC(is_amv_value_or_view_class::value&& has_scalar_or_string_value_type::value) h5_write(h5::group g, std::string const& name, ArrayType const& A) { + if (A.is_empty()) TRIQS_RUNTIME_ERROR << " Can not save an empty array into hdf5"; 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/parameters/opaque_object_h5.cpp b/triqs/parameters/opaque_object_h5.cpp index 4f7dd339..2a68d42e 100644 --- a/triqs/parameters/opaque_object_h5.cpp +++ b/triqs/parameters/opaque_object_h5.cpp @@ -48,7 +48,8 @@ namespace triqs { namespace utility { size_t type_hash_element = type_hash; auto it2= _object::code_element_rank_to_code_array.find(std::make_pair(type_hash_element,rank)); if (it2 == _object::code_element_rank_to_code_array.end()) - TRIQS_RUNTIME_ERROR << " code_element_rank_to_code_array : type not found" << rank; + TRIQS_RUNTIME_ERROR << " code_element_rank_to_code_array : type not found : " << name << " " << type_hash_element << " " + << rank; type_hash = it2->second; } } diff --git a/triqs/utility/mini_vector.hpp b/triqs/utility/mini_vector.hpp index 5ba82ba9..de735c95 100644 --- a/triqs/utility/mini_vector.hpp +++ b/triqs/utility/mini_vector.hpp @@ -41,7 +41,7 @@ namespace triqs { namespace utility { public : typedef T value_type; - + mini_vector(){init();} #define AUX(z,p,unused) _data[p] = x_##p; @@ -69,6 +69,8 @@ namespace triqs { namespace utility { friend void swap(mini_vector & a, mini_vector & b) { std::swap(a._data, b._data);} + int size() const { return Rank;} + template mini_vector & operator=(const mini_vector & x){ for (int i=0;i