diff --git a/triqs/arrays/indexmaps/range.hpp b/triqs/arrays/indexmaps/range.hpp index 43ef283f..c02cb98d 100644 --- a/triqs/arrays/indexmaps/range.hpp +++ b/triqs/arrays/indexmaps/range.hpp @@ -40,6 +40,8 @@ namespace arrays { range(std::ptrdiff_t first__, std::ptrdiff_t last__, std::ptrdiff_t step__ = 1) : first_(first__), last_(last__), step_(step__) {} + range(std::ptrdiff_t i) : range (i,i+1,1){} + std::ptrdiff_t first() const { return first_; } std::ptrdiff_t last() const { return last_; } std::ptrdiff_t step() const { return step_; } diff --git a/triqs/gfs/block.hpp b/triqs/gfs/block.hpp index ce0d6598..cd4ab913 100644 --- a/triqs/gfs/block.hpp +++ b/triqs/gfs/block.hpp @@ -97,13 +97,13 @@ namespace gfs { block_gf make_block_gf(int n, gf const &g) { auto V = std::vector>{}; for (int i = 0; i < n; ++i) V.push_back(g); - return {{n}, std::move(V), nothing{}, nothing{}}; + return {{n}, std::move(V), nothing{}, nothing{}, nothing{}}; } // from a vector of gf (moving directly) template block_gf make_block_gf(std::vector> V) { - return {{int(V.size())}, std::move(V), nothing{}, nothing{}}; + return {{int(V.size())}, std::move(V), nothing{}, nothing{}, nothing{}}; } // from a vector of gf : generalized to have a different type of gf in the vector (e.g. views...) @@ -111,13 +111,13 @@ namespace gfs { block_gf make_block_gf(std::vector const &V) { auto V2 = std::vector>{}; for (auto const &g : V) V2.push_back(g); - return {{int(V.size())}, std::move(V2), nothing{}, nothing{}}; + return {{int(V.size())}, std::move(V2), nothing{}, nothing{}, nothing{}}; } // from a init list of GF with the correct type template block_gf make_block_gf(std::initializer_list> const &V) { - return {{int(V.size())}, V, nothing{}, nothing{}}; + return {{int(V.size())}, V, nothing{}, nothing{}, nothing{}}; } // from vector and a gf to be copied @@ -125,7 +125,7 @@ namespace gfs { block_gf make_block_gf(std::vector block_names, gf const &g) { auto V = std::vector>{}; for (int i = 0; i < block_names.size(); ++i) V.push_back(g); - return {{block_names}, std::move(V), nothing{}, nothing{}}; + return {{block_names}, std::move(V), nothing{}, nothing{}, nothing{}}; } // from vector, vector @@ -133,7 +133,7 @@ namespace gfs { block_gf make_block_gf(std::vector block_names, std::vector> V) { if (block_names.size() != V.size()) TRIQS_RUNTIME_ERROR << "make_block_gf(vector, vector) : the two vectors do not have the same size !"; - return {{block_names}, std::move(V), nothing{}, nothing{}}; + return {{block_names}, std::move(V), nothing{}, nothing{}, nothing{}}; } // from vector, init_list @@ -141,7 +141,7 @@ namespace gfs { block_gf make_block_gf(std::vector block_names, std::initializer_list> const &V) { if (block_names.size() != V.size()) TRIQS_RUNTIME_ERROR << "make_block_gf(vector, init_list) : size mismatch !"; - return {{block_names}, V, nothing{}, nothing{}}; + return {{block_names}, V, nothing{}, nothing{}, nothing{}}; } // ------------------------------- Free Factories for view type -------------------------------------------------- @@ -149,17 +149,17 @@ namespace gfs { template gf_view::type::view_type> make_block_gf_view(G0 &&g0, G &&... g) { auto V = std::vector::type::view_type>{std::forward(g0), std::forward(g)...}; - return {{int(V.size())}, std::move(V), nothing{}, nothing{}}; + return {{int(V.size())}, std::move(V), nothing{}, nothing{}, nothing{}}; // return { gf_mesh {int(V.size())}, std::move(V), nothing{}, nothing{} } ; } template gf_view make_block_gf_view_from_vector(std::vector V) { - return {{int(V.size())}, std::move(V), nothing{}, nothing{}}; + return {{int(V.size())}, std::move(V), nothing{}, nothing{}, nothing{}}; } // for cython proxy only. do not document. template gf_view make_block_gf_view_from_vector_of_cython_proxy(std::vector V) { - return {{int(V.size())}, std::move(V), nothing{}, nothing{}}; + return {{int(V.size())}, std::move(V), nothing{}, nothing{}, nothing{}}; } // ------------------------------- Extend reinterpret_scalar_valued_gf_as_matrix_valued for block gf ------ diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp index 554b5a9c..e134086b 100644 --- a/triqs/gfs/gf.hpp +++ b/triqs/gfs/gf.hpp @@ -72,6 +72,15 @@ namespace gfs { using type = nothing; }; + // indices + template struct indices { + using type = nothing; + }; + + template struct indices { + using type = indices_2; + }; + // data_proxy contains function to manipulate the data array, but no data itself. // this is used to specialize this part of the code to array of dim 3 (matrix gf), dim 1 (scalar gf) and vector (e.g. block gf, // ...) @@ -119,6 +128,7 @@ namespace gfs { using mesh_point_t = typename mesh_t::mesh_point_t; using mesh_index_t = typename mesh_t::index_t; using symmetry_t = typename gfs_implementation::symmetry::type; + using indices_t = typename gfs_implementation::indices::type; using evaluator_t = gfs_implementation::evaluator; using data_proxy_t = gfs_implementation::data_proxy; @@ -138,6 +148,7 @@ namespace gfs { singularity_t &singularity() { return _singularity; } singularity_t const &singularity() const { return _singularity; } symmetry_t const &symmetry() const { return _symmetry; } + indices_t const &indices() const { return _indices; } evaluator_t const &get_evaluator() const { return _evaluator; } auto get_data_shape() const DECL_AND_RETURN(get_shape(this -> data())); @@ -147,6 +158,10 @@ namespace gfs { data_t _data; singularity_t _singularity; symmetry_t _symmetry; + indices_t _indices; + public : + std::string name; + protected : evaluator_t _evaluator; data_proxy_t _data_proxy; @@ -160,6 +175,8 @@ namespace gfs { _data(factory(x.data())), _singularity(factory(x.singularity())), _symmetry(x.symmetry()), + _indices(x.indices()), + name(x.name), _evaluator(x._evaluator) {} gf_impl(gf_impl &&) = default; @@ -171,15 +188,21 @@ namespace gfs { _data(factory(x.data())), _singularity(factory(x.singularity())), _symmetry(x.symmetry()), + _indices(x.indices()), + name(x.name), _evaluator(x.get_evaluator()) {} template - gf_impl(M &&m, D &&dat, S &&sing, SY &&sy, EV &&ev) - : _mesh(std::forward(m)), - _data(std::forward(dat)), - _singularity(std::forward(sing)), - _symmetry(std::forward(sy)), - _evaluator(std::forward(ev)) {} + gf_impl(M &&m, D &&dat, S &&sing, SY &&sy, indices_t ind, std::string name, EV &&ev) + : _mesh(std::forward(m)) + , _data(std::forward(dat)) + , _singularity(std::forward(sing)) + , _symmetry(std::forward(sy)) + , _indices(std::move(ind)) + , name(name) + , _evaluator(std::forward(ev)) { + if (!_indices.check_size(this)) TRIQS_RUNTIME_ERROR << "Size of indices mismatch with data size"; + } void operator=(gf_impl const &rhs) = delete; // done in derived class. @@ -189,6 +212,8 @@ namespace gfs { swap(this->_data, b._data); swap(this->_singularity, b._singularity); swap(this->_symmetry, b._symmetry); + swap(this->_indices, b._indices); + swap(this->name, b.name); swap(this->_evaluator, b._evaluator); } @@ -320,6 +345,8 @@ namespace gfs { ar &TRIQS_MAKE_NVP("singularity", _singularity); ar &TRIQS_MAKE_NVP("mesh", _mesh); ar &TRIQS_MAKE_NVP("symmetry", _symmetry); + ar &TRIQS_MAKE_NVP("indices", _indices); + ar &TRIQS_MAKE_NVP("name", name); } /// print @@ -391,17 +418,19 @@ namespace gfs { *this = x; } - gf(typename B::mesh_t m, typename B::data_t dat, typename B::singularity_view_t const &si, typename B::symmetry_t const &s) - : B(std::move(m), std::move(dat), si, s, typename B::evaluator_t{}) {} + gf(typename B::mesh_t m, typename B::data_t dat, typename B::singularity_view_t const &si, typename B::symmetry_t const &s, + typename B::indices_t const &ind, std::string name = "") + : B(std::move(m), std::move(dat), si, s, ind, name, typename B::evaluator_t{}) {} using target_shape_t = typename factory::target_shape_t; - gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{}) - : B(std::move(m), factory::make_data(m, shape), factory::make_singularity(m, shape), typename B::symmetry_t{}, + gf(typename B::mesh_t m, target_shape_t shape = target_shape_t{}, typename B::indices_t const &ind = typename B::indices_t{}) + : B(std::move(m), factory::make_data(m, shape), factory::make_singularity(m, shape), typename B::symmetry_t{}, ind, {}, // clean unncessary types typename B::evaluator_t{}) {} friend void swap(gf &a, gf &b) noexcept { a.swap_impl(b); } + // CLEAN THIS : one op only gf rhs, + move gf &operator=(gf const &rhs) { *this = gf(rhs); return *this; @@ -424,6 +453,7 @@ namespace gfs { this->_singularity = rhs.singularity(); // to be implemented : there is none in the gf_expr in particular.... // this->_symmetry = rhs.symmetry(); + // indices and name are not affected by it ??? } }; @@ -454,14 +484,17 @@ namespace gfs { gf_view(gf_impl &&g) noexcept : B(std::move(g), bool {}) {} // from a gf && template - gf_view(typename B::mesh_t const &m, D const &dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s) - : B(m, factory(dat), t, s, typename B::evaluator_t{}) {} + gf_view(typename B::mesh_t const &m, D const &dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s, + typename B::indices_t const &ind, std::string name = "") + : B(m, factory(dat), t, s, ind, name, typename B::evaluator_t{}) {} void rebind(gf_view const &X) noexcept { this->_mesh = X._mesh; this->_symmetry = X._symmetry; this->_data_proxy.rebind(this->_data, X); this->_singularity.rebind(X._singularity); + this->_indices = X._indices; + this->name = X.name; } void rebind(gf_view const &X) noexcept { @@ -489,8 +522,9 @@ namespace gfs { gf_view(gf_impl &&g) noexcept : B(std::move(g), bool {}) {} // from a gf && template - gf_view(typename B::mesh_t const &m, D const &dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s) - : B(m, factory(dat), t, s, typename B::evaluator_t{}) {} + gf_view(typename B::mesh_t const &m, D const &dat, typename B::singularity_view_t const &t, typename B::symmetry_t const &s, + typename B::indices_t const &ind = typename B::indices_t{}, std::string name = "") + : B(m, factory(dat), t, s, ind, name, typename B::evaluator_t{}) {} friend void swap(gf_view &a, gf_view &b) noexcept { a.swap_impl(b); } @@ -499,6 +533,8 @@ namespace gfs { this->_symmetry = X._symmetry; this->_data_proxy.rebind(this->_data, X); this->_singularity.rebind(X._singularity); + this->_indices = X._indices; + this->name = X.name; } gf_view &operator=(gf_view const &rhs) { @@ -540,7 +576,7 @@ namespace gfs { static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); using arrays::range; return {g.mesh(), g.data()(range(), std::forward(args)...), - slice_target(g.singularity(), std::forward(args)...), g.symmetry()}; + slice_target(g.singularity(), std::forward(args)...), g.symmetry(), slice(g.indices(),std::forward(args)...), g.name }; } template @@ -560,7 +596,7 @@ namespace gfs { static_assert(std::is_same::value, "slice_target only for matrix_valued GF's"); using arrays::range; return {g.mesh(), g.data()(range(), std::forward(args)...), - slice_target(g.singularity(), range(args, args + 1)...), g.symmetry()}; + slice_target(g.singularity(), range(args, args + 1)...), g.symmetry(), {}, g.name}; } template @@ -573,6 +609,15 @@ namespace gfs { return slice_target_to_scalar(g(), std::forward(args)...); } + /*// slice with string arguments, for python mainly. + template + gf_view slice_target_string_indices(gf_view g, std::string const & a1, std::string const & a2) { + int n1 = g.indices().convert_index(a1,0); + int n2 = g.indices().convert_index(a2,1); + return slice_target (g, n1, n2); // the slice will have no indices. Ok ? + } +*/ + // ---------------------------------- target reinterpretation ------------------------------------ // a scalar_valued gf can be viewed as a 1x1 matrix @@ -581,7 +626,7 @@ namespace gfs { reinterpret_scalar_valued_gf_as_matrix_valued(gf_view g) { using a_t = typename gf_view::data_view_t; auto a = a_t{typename a_t::indexmap_type(join(g.data().shape(), make_shape(1, 1))), g.data().storage()}; - return {g.mesh(), a, g.singularity(), g.symmetry()}; + return {g.mesh(), a, g.singularity(), g.symmetry(), {} , g.name}; } template @@ -677,6 +722,8 @@ namespace gfs { h5_write(gr, "singularity", g._singularity); h5_write(gr, "mesh", g._mesh); h5_write(gr, "symmetry", g._symmetry); + h5_write(gr, "indices", g._indices); + //h5_write(gr, "name", g.name); } template static void read(h5::group gr, gf_impl &g) { @@ -684,6 +731,8 @@ namespace gfs { h5_read(gr, "singularity", g._singularity); h5_read(gr, "mesh", g._mesh); h5_read(gr, "symmetry", g._symmetry); + h5_read(gr, "indices", g._indices); + //h5_read(gr, "name", g.name); } }; } // gfs_implementation diff --git a/triqs/gfs/local/no_tail.hpp b/triqs/gfs/local/no_tail.hpp index 80be390d..6b38cc05 100644 --- a/triqs/gfs/local/no_tail.hpp +++ b/triqs/gfs/local/no_tail.hpp @@ -30,7 +30,7 @@ namespace gfs { template gf_view make_gf_view_without_tail(gf_impl const &g) { - return {g.mesh(), g.data(), {}, g.symmetry()}; + return {g.mesh(), g.data(), {}, g.symmetry(), g.indices(), g.name}; } namespace details { // dispatch the test for scalar_valued and matrix_valued diff --git a/triqs/gfs/tools.hpp b/triqs/gfs/tools.hpp index 3b7294ac..9a0276b9 100644 --- a/triqs/gfs/tools.hpp +++ b/triqs/gfs/tools.hpp @@ -74,6 +74,73 @@ namespace gfs { template closest_pt_wrap closest_mesh_pt(T &&... x) { return closest_pt_wrap{std::forward(x)...}; } + //------------------------------------------------------ + + // A simple indice struct + // Replace empty state by optional ? + struct indices_2 { + + std::vector> ind; + + indices_2() = default; + indices_2(std::vector> _ind) : ind(std::move(_ind)) {} + + bool is_empty() const { return ind.size() == 0; } + + template bool check_size(G *g) const { + return (is_empty() || + ((ind.size() == 2) && (ind[0].size() == get_target_shape(*g)[0]) && (ind[1].size() == get_target_shape(*g)[1]))); + } + + std::vector operator[](int i) const { + if (is_empty()) + return {}; + else + return ind[i]; + } + + arrays::range convert_index(std::string const &s, int l_r) const { + if (!is_empty()) { + auto b = ind[l_r].begin(), e = ind[l_r].end(); + auto it = std::find(b, e, s); + if (it != e) return it - b; + } + TRIQS_RUNTIME_ERROR << "There are no string indices for this Green function"; + } + + friend class boost::serialization::access; + template void serialize(Archive &ar, const unsigned int version) { ar &TRIQS_MAKE_NVP("ind", ind); } + + friend void h5_write(h5::group fg, std::string subgroup_name, indices_2 const &g) { + if (g.is_empty()) return; + auto gr = fg.create_group(subgroup_name); + h5_write(gr, "left", g.ind[0]); + h5_write(gr, "right", g.ind[1]); + } + + friend void h5_read(h5::group fg, std::string subgroup_name, indices_2 &g) { + h5::group gr; + try { + gr = fg.open_group(subgroup_name); + } + catch (...) { + g = indices_2{}; // empty, no file + return; + } + h5_read(gr, "left", g.ind[0]); + h5_read(gr, "right", g.ind[1]); + } + }; + + // inline indices_2 slice(indices_2 const & ind, arrays::range rl, arrays::range rr) { + // return {}; + // } + + inline indices_2 slice(indices_2 const &ind, arrays::range rl, arrays::range rr) { + if (ind.is_empty()) return {}; + return {}; // MODIFY : slice the indices !!! + TRIQS_RUNTIME_ERROR << "Not implemented : slice of string indices"; + } //------------------------------------------------------ @@ -87,10 +154,12 @@ namespace gfs { template void operator=(RHS &&) {} friend void h5_write(h5::group, std::string subgroup_name, nothing) {} friend void h5_read(h5::group, std::string subgroup_name, nothing) {} + template friend nothing slice(nothing, A...) { return nothing(); } friend class boost::serialization::access; template void serialize(Archive &ar, const unsigned int version) {} friend nothing operator+(nothing, nothing) { return nothing(); } template friend void assign_from_expression(nothing &, RHS) {} + template bool check_size(A) {return true;} }; template nothing slice_target(nothing, T...) { return nothing(); }