3
0
mirror of https://github.com/triqs/dft_tools synced 2024-12-24 13:23:37 +01:00

Fix g(k,om) for tests

- evaluator
- G(k,tau) is real
- partial_eval for matrix_valued functions
- details : simplifying traits (using decay_t)
This commit is contained in:
Olivier Parcollet 2014-02-26 15:40:14 +01:00
parent 300b016b75
commit e1c113b745
7 changed files with 63 additions and 35 deletions

View File

@ -2,16 +2,20 @@
#include <triqs/gfs.hpp>
#include <triqs/gfs/bz.hpp>
// FAIRE make_value !!
//
using namespace triqs::gfs;
using namespace triqs::clef;
using namespace triqs::arrays;
using namespace triqs::lattice;
template <typename Function, typename Mesh>
// requires ( is_function_on_mesh<Function,Mesh>())
auto sum(Function const &f, Mesh const &m)->decltype(f(typename Mesh::mesh_point_t{})) {
auto sum(Function const &f, Mesh const &m) ->decltype(make_matrix(0*f(*(m.begin())))) {
//auto sum(Function const &f, Mesh const &m) {
auto res = decltype(f(typename Mesh::mesh_point_t{})) {};
for (auto const &x : m) res = res + f(x);
//auto res = typename triqs::regular_type_if_exists_else_type<decltype(f(typename Mesh::mesh_point_t{}))>::type (f(m.begin()));
auto res = make_matrix(0*f(*(m.begin())));
for (auto const &x : m) res = res + f(x);
return res;
}
namespace triqs {
@ -21,7 +25,7 @@ namespace triqs {
}
}
#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
#define TEST(...) std::cout << BOOST_PP_STRINGIZE((__VA_ARGS__)) << " ---> " << (__VA_ARGS__) << std::endl << std::endl;
int main() {
try {
@ -41,16 +45,25 @@ int main() {
G_k_iom(k_, w_) << 1 / (w_ - eps(k_));
auto G_loc = gf<imfreq>{{beta,Fermion, 100}};
//G_loc(w_) << sum(k_ >> G_k_iom(k_,w_), g_eps.mesh());//does not compile
auto G_loc = gf<imfreq, matrix_valued, no_tail>{{beta, Fermion, 100}, {1, 1}};
auto r = G_k_iom(k_t{0, 0}, matsubara_freq{0, beta, Fermion});
auto r5 = sum(k_ >> G_k_iom(k_,0), g_eps.mesh());
G_loc(w_) << sum(k_ >> G_k_iom(k_,w_), g_eps.mesh());
TEST(G_loc(0));
auto G_k_tau = gf<cartesian_product<bz, imtime>>{{{bz_, 20}, {beta, Fermion, 100}}, {1, 1}};
//curry<1>(G_k_tau) [k_] << inverse_fourier(curry<1>(G_k_iom)[k_]); //does not compile
//auto r3 = partial_eval<0>(G_k_iom,0);
//auto r4 = partial_eval<0>(G_k_tau,0);
//auto gt = curry<0>(G_k_tau) [0];
//auto gw = curry<0>(G_k_iom)[0];
//curry<0>(G_k_tau) [k_] << inverse_fourier(curry<0>(G_k_iom)[k_]);
//TEST(G_k_tau(0,0));
//TEST(G_k_tau[{0,0}]);
// hdf5
//H5::H5File file("ess_g_k_om.h5", H5F_ACC_TRUNC );

View File

@ -170,8 +170,7 @@ namespace triqs { namespace clef {
}; \
} \
template <typename L, typename R> \
typename std::enable_if<is_any_lazy<L, R>::value, expr<tags::TAG, expr_storage_t<L>, expr_storage_t<R>>>::type operator OP( \
L&& l, R&& r) { \
std14::enable_if_t<is_any_lazy<L, R>::value, expr<tags::TAG, expr_storage_t<L>, expr_storage_t<R>>> operator OP(L&& l, R&& r) { \
return {tags::TAG(), std::forward<L>(l), std::forward<R>(r)}; \
} \
template <> struct operation<tags::TAG> { \
@ -198,7 +197,7 @@ namespace triqs { namespace clef {
}; \
} \
template <typename L> \
typename std::enable_if<is_any_lazy<L>::value, expr<tags::TAG, expr_storage_t<L>>>::type operator OP(L&& l) { \
std14::enable_if_t<is_any_lazy<L>::value, expr<tags::TAG, expr_storage_t<L>>> operator OP(L&& l) { \
return {tags::TAG(), std::forward<L>(l)}; \
} \
template <> struct operation<tags::TAG> { \
@ -223,11 +222,7 @@ namespace triqs { namespace clef {
/* ---------------------------------------------------------------------------------------------------
* Evaluation of the expression tree.
* --------------------------------------------------------------------------------------------------- */
//template<typename ... T> struct _or;
//template<typename T0, typename ... T> struct _or<T0,T...> : std::integral_constant<bool,T0::value || _or<T...>::value>{};
//template<> struct _or<> : std::false_type{};
// just to try
constexpr bool __or() { return false;}
template<typename ... B> constexpr bool __or(bool b, B... bs) { return b || __or(bs...); }
@ -283,8 +278,8 @@ namespace triqs { namespace clef {
static constexpr bool is_lazy = __or(evaluator<Childs, Pairs...>::is_lazy...);
auto operator()(expr<Tag, Childs...> const& ex, Pairs const&... pairs) const {
auto eval_in_context = [&pairs](auto const& _child) { return eval(_child, pairs); };
return triqs::tuple::apply_compose(op_dispatch<Tag, is_lazy>{}, eval_in_context, ex.childs);
auto eval_in_context = [&pairs...](auto const& _child) { return eval(_child, pairs...); };
return tuple::apply_compose(op_dispatch<Tag, is_lazy>{}, eval_in_context, ex.childs);
}
};
#else
@ -299,8 +294,10 @@ namespace triqs { namespace clef {
template <typename Tag, typename... Childs, typename... Pairs> struct evaluator_node_gal<2, expr<Tag, Childs...>, Pairs...> {
static constexpr bool is_lazy = __or(evaluator<Childs, Pairs...>::is_lazy...);
auto operator()(expr<Tag, Childs...> const& ex, Pairs const&... pairs) const DECL_AND_RETURN(op_dispatch<Tag, is_lazy> {
}(eval(std::get<0>(ex.childs), pairs...), eval(std::get<1>(ex.childs), pairs...)));
decltype(auto) operator()(expr<Tag, Childs...> const &ex,
Pairs const &... pairs) const {
return op_dispatch<Tag, is_lazy> {}(eval(std::get<0>(ex.childs), pairs...), eval(std::get<1>(ex.childs), pairs...));
}
};
// the general case for more than 2 nodes. I put 1 and 2 nodes apart, just because it is the most frequent

View File

@ -50,10 +50,13 @@ namespace gfs {
// simple evaluation : take the point on the grid...
template <> struct evaluator_fnt_on_mesh<bz> {
lattice::k_t k;
size_t n;
evaluator_fnt_on_mesh() = default;
template <typename MeshType> evaluator_fnt_on_mesh(MeshType const &m, lattice::k_t const &k_) { k = k_; }
template <typename F> auto operator()(F const &f) const DECL_AND_RETURN(f(k));
template <typename MeshType> evaluator_fnt_on_mesh(MeshType const &m, lattice::k_t const &k) {
n = m.locate_neighbours(k); // TO BE IMPROVED
}
//template <typename F> auto operator()(F const &f) const DECL_AND_RETURN(f(k));
template <typename F> decltype(auto) operator()(F const &f) const { return f(n); }
};
// ------------- evaluator -------------------

View File

@ -44,7 +44,7 @@ namespace triqs { namespace gfs {
template<typename ... Ms> struct cart_prod_impl;
template<typename ... Ms> using cart_prod = typename cart_prod_impl<Ms...>::type;
template<typename ... Ms> struct cart_prod_impl<std::tuple<Ms...>> { using type = cartesian_product<Ms...>;};
template<typename M> struct cart_prod_impl<std::tuple<M>> { typedef M type;}; //using type = M;};
template<typename M> struct cart_prod_impl<std::tuple<M>> { using type = M;};
template<typename M0, typename M1, typename ...M> auto rm_tuple_of_size_one(std::tuple<M0,M1,M...> const & t) DECL_AND_RETURN(t);
template<typename M> auto rm_tuple_of_size_one(std::tuple<M> const & t) DECL_AND_RETURN(std::get<0>(t));
@ -55,9 +55,9 @@ namespace triqs { namespace gfs {
template<typename ... T> std::tuple<T...> const & as_tuple(std::tuple<T...> const & x) { return x;}
template<typename ... T> std::tuple<T...> & as_tuple(std::tuple<T...> & x) { return x;}
template<int ... pos, typename Opt, typename Target, bool B, bool IsConst, typename IT, typename ... Ms>
gf_view< typename cart_prod_impl< triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>::type ,Target, Opt,IsConst>
partial_eval(gf_impl< cartesian_product<Ms...>, Target,Opt,B,IsConst> const & g, IT index) {
template <int... pos, typename Opt, typename Target, bool IsConst, typename IT, typename... Ms>
gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, IsConst>
partial_eval(gf_view<cartesian_product<Ms...>, Target, Opt, IsConst> g, IT index) {
// meshes of the returned gf_view : just drop the mesh of the evaluated variables
auto meshes_tuple_partial = triqs::tuple::filter_out<pos...>(g.mesh().components());
// a view of the array of g, with the dimension sizeof...(Ms)
@ -65,12 +65,24 @@ namespace triqs { namespace gfs {
// now rebuild a tuple of the size sizeof...(Ms), containing the indices and range at the position of evaluated variables.
auto arr_args = triqs::tuple::inverse_filter<sizeof...(Ms),pos...>(as_tuple(index), arrays::range());
// from it, we make a slice of the array of g, corresponding to the data of the returned gf_view
auto arr2 = triqs::tuple::apply(arr, arr_args);
auto arr2 = triqs::tuple::apply(arr, std::tuple_cat(arr_args, std::make_tuple(arrays::ellipsis{})));
// finally, we build the view on this data.
using r_t = gf_view< cart_prod< triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>> ,Target, Opt,IsConst>;
return r_t{ rm_tuple_of_size_one(meshes_tuple_partial), arr2, typename r_t::singularity_non_view_t{}, typename r_t::symmetry_t{} };
}
template <int... pos, typename Opt, typename Target, typename IT, typename... Ms>
gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, false>
partial_eval(gf<cartesian_product<Ms...>, Target, Opt> & g, IT index) {
return partial_eval<pos...>(g(),index);
}
template <int... pos, typename Opt, typename Target, typename IT, typename... Ms>
gf_view<cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>, Target, Opt, true>
partial_eval(gf<cartesian_product<Ms...>, Target, Opt> const& g, IT index) {
return partial_eval<pos...>(g(),index);
}
/// --------------------------- curry ---------------------------------
// curry<0>(g) returns : x-> y... -> g(x,y...)
// curry<1>(g) returns : y-> x,z... -> g(x,y,z...)

View File

@ -92,7 +92,10 @@ namespace gfs {
struct data_proxy<cartesian_product<Ms...>, matrix_valued, Opt> : data_proxy_array<std::complex<double>, 3> {};
template <int R, typename Opt, typename... Ms>
struct data_proxy<cartesian_product<Ms...>, tensor_valued<R>, Opt> : data_proxy_array<std::complex<double>, R + 1> {};
template <typename Opt, typename M0>
struct data_proxy<cartesian_product<M0,imtime>, matrix_valued, Opt> : data_proxy_array<double, 3> {};
/// --------------------------- evaluator ---------------------------------
/**

View File

@ -45,7 +45,7 @@ namespace lattice {
}
pos++;
}
return pos;
return pos_min;
}
/// Write into HDF5

View File

@ -35,15 +35,15 @@ namespace std {
template<typename ... T> _triqs_reversed_tuple<std::tuple<T...>const &> reverse(std::tuple<T...> const & x) { return {x};}
template<int pos, typename TU> auto get(_triqs_reversed_tuple<TU> const & t)
DECL_AND_RETURN(std::get<std::tuple_size<typename std::remove_const<typename std::remove_reference<TU>::type>::type>::value-1-pos>(t._x));
DECL_AND_RETURN(std::get<std::tuple_size<std14::decay_t<TU>>::value-1-pos>(t._x));
template<int pos, typename TU> auto get(_triqs_reversed_tuple<TU> & t)
DECL_AND_RETURN(std::get<std::tuple_size<typename std::remove_const<typename std::remove_reference<TU>::type>::type>::value-1-pos>(t._x));
DECL_AND_RETURN(std::get<std::tuple_size<std14::decay_t<TU>>::value-1-pos>(t._x));
template<int pos, typename TU> auto get(_triqs_reversed_tuple<TU> && t)
DECL_AND_RETURN(std::get<std::tuple_size<typename std::remove_const<typename std::remove_reference<TU>::type>::type>::value-1-pos>(move(t)._x));
DECL_AND_RETURN(std::get<std::tuple_size<std14::decay_t<TU>>::value-1-pos>(move(t)._x));
template<typename TU> class tuple_size<_triqs_reversed_tuple<TU>> : public tuple_size<typename std::remove_const<typename std::remove_reference<TU>::type>::type>{};
template<typename TU> class tuple_size<_triqs_reversed_tuple<TU>> : public tuple_size<std14::decay_t<TU>>{};
}
@ -177,7 +177,7 @@ namespace triqs { namespace tuple {
};
template<typename F, typename T1>
auto apply_on_tuple (F && f,T1 && t1) DECL_AND_RETURN( apply_on_tuple_impl<std::tuple_size<typename std::remove_const<typename std::remove_reference<T1>::type>::type>::value-1>()(std::forward<F>(f),std::forward<T1>(t1)));
auto apply_on_tuple (F && f,T1 && t1) DECL_AND_RETURN( apply_on_tuple_impl<std::tuple_size<std14::decay_t<T1>>::value-1>()(std::forward<F>(f),std::forward<T1>(t1)));
/**
* apply_on_zip(f, t1,t2)
@ -197,7 +197,7 @@ namespace triqs { namespace tuple {
};
template<typename F, typename T1, typename T2>
auto apply_on_zip (F && f,T1 && t1, T2 && t2) DECL_AND_RETURN( apply_on_zip_impl<std::tuple_size<typename std::remove_const<typename std::remove_reference<T1>::type>::type>::value-1>()(std::forward<F>(f),std::forward<T1>(t1),std::forward<T2>(t2)));
auto apply_on_zip (F && f,T1 && t1, T2 && t2) DECL_AND_RETURN( apply_on_zip_impl<std::tuple_size<std14::decay_t<T1>>::value-1>()(std::forward<F>(f),std::forward<T1>(t1),std::forward<T2>(t2)));
/**
* apply_on_zip(f, t1,t2,t3)