diff --git a/test/triqs/gf/gf_re_im_freq_time.cpp b/test/triqs/gf/gf_re_im_freq_time.cpp index 55f264b1..d2601806 100644 --- a/test/triqs/gf/gf_re_im_freq_time.cpp +++ b/test/triqs/gf/gf_re_im_freq_time.cpp @@ -3,6 +3,8 @@ #include #include #include +#include +#include #include #include @@ -11,7 +13,9 @@ namespace tql= triqs::clef; using namespace triqs::gfs; int main() { - + +try { + double precision=10e-9; double beta =1.; @@ -25,16 +29,35 @@ int main() { int n_re_freq=100; int n_im_freq=100; - auto G_t_tau= make_gf( tmin, tmax, n_re_time, beta, Fermion, n_im_time); + //auto G_t_tau= make_gf( tmin, tmax, n_re_time, beta, Fermion, n_im_time); auto G_w_wn = make_gf( wmin, wmax, n_re_freq, beta, Fermion, n_im_freq); auto G_w_tau= make_gf(wmin, wmax, n_re_freq, beta, Fermion, n_im_time); auto G_w= make_gf(wmin, wmax, n_re_freq); - + + auto G_t_tau= make_gf, scalar_valued>( make_gf_mesh(tmin, tmax, n_re_time), make_gf_mesh(beta, Fermion, n_im_time)); + //auto G_t_tau_N= make_gf, scalar_valued>( {tmin, tmax, n_re_time}, {beta, Fermion, n_im_time}); + + auto G_w_wn2 = make_gf, scalar_valued>( make_gf_mesh(wmin, wmax, n_re_freq), make_gf_mesh(beta, Fermion, n_im_freq)); + auto G_w_tau2 = make_gf, scalar_valued>( make_gf_mesh(wmin, wmax, n_re_freq), make_gf_mesh(beta, Fermion, n_im_time,full_bins)); + + + //auto g_tau = slice_mesh1(G_w_tau(),1); + //auto g_wn = G_w_wn_curry0[1]; + + //std::cout << G_t_tau_N (0.1,0.2) << std::endl; + + auto G_w_wn2_view = G_w_wn2(); + auto G_w_wn_sl0_a = partial_eval<0>(G_w_wn2(), std::make_tuple(8)); + //auto G_w_wn_curry0_a = curry0(G_w_wn2); + //auto G_w_wn_sl0_a = slice_mesh0(G_w_wn2(), 8); + triqs::clef::placeholder<0> w_; triqs::clef::placeholder<1> wn_; triqs::clef::placeholder<2> tau_; G_w_wn(w_,wn_)<<1/(wn_-1)/( pow(w_,3) ); + G_w_wn2(w_,wn_)<<1/(wn_-1)/( pow(w_,3) ); G_w_tau(w_,tau_)<< exp( -2*tau_ ) / (w_*w_ + 1 ); + G_w_tau2(w_,tau_)<< exp( -2*tau_ ) / (w_*w_ + 1 ); int index = n_re_freq/3; double tau = std::get<1>(G_w_tau.mesh().components())[index]; @@ -45,21 +68,38 @@ int main() { G_w.singularity()(0)=triqs::arrays::matrix{{0}}; G_w.singularity()(1)=triqs::arrays::matrix{{0}}; G_w.singularity()(2)=triqs::arrays::matrix{{0}}; + //auto G_w2 = slice_mesh1(G_w_tau(), index); auto G_w2 = slice_mesh_imtime(G_w_tau, index); for(auto& w:G_w.mesh()) if ( std::abs(G_w[w]-G_w2[w]) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice error : w="<< w <<" ,G_w="<< G_w[w]<<" ,G_w2="<< G_w2[w] <<"\n"; //test of the interpolation std::cout << G_t_tau(0.789,0.123) << std::endl; - std::cout << G_w_wn( 0.789,0.123) << std::endl; - std::cout << G_w_tau(0.789,0.123) << std::endl; - + std::cout << "G_w_wn( 0.789,0.123) "<< G_w_wn( 0.789,0.123) << std::endl; + std::cout << "G_w_wn( 0.789,0.123) "<(G_w_wn2); + auto G_w_wn_curry1 = curry<1>(G_w_wn2); + auto G_w_wn2_view2 = G_w_wn2(); + + std::cout << " curry "< precision) TRIQS_RUNTIME_ERROR<<" fourier_slice_re_time error : t="<< t <<" ,G_t="<< G_t[t] <<" ,G_t2="<< G_t2[t] <<"\n"; - +*/ +} +catch(std::exception const & e ) { std::cout << "error "<< e.what()<< std::endl;} } diff --git a/triqs/gfs/curry.hpp b/triqs/gfs/curry.hpp new file mode 100644 index 00000000..05dcf36c --- /dev/null +++ b/triqs/gfs/curry.hpp @@ -0,0 +1,160 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 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 . + * + ******************************************************************************/ +#ifndef TRIQS_GF_CURRY_H +#define TRIQS_GF_CURRY_H +#include "./product.hpp" + +namespace triqs { namespace gfs { + + template struct lambda_valued {}; + + namespace gfs_implementation { + + /// --------------------------- data access --------------------------------- + + template struct data_proxy,Opt> : data_proxy_lambda {}; + + /// --------------------------- Factories --------------------------------- + + template + struct factories, lambda_valued, Opt> {}; + + // detail + template struct cartesian_product_add_front; + template struct cartesian_product_add_front>{ typedef cartesian_product type; }; + + // ------------------------------------------------- + // Partial evaluation of the gf + // ------------------------------------------------- + // + // Given a cartesian_product of meshes (CP), and a compile time list of int (I) + // - metacompute the list of Ms without position those at position 0,2 (type) + // - provide 2 runtimes functions : + // - sl : given empty tuple () and a tuple (it) of indices + // return the tuple of indices and range, where range are at the position defined by I, + // and the indices in other places, in order. + // - m : returns from a CP object the corresponding tuple of meshes of the remaining meshes + // after partial eval (of the type computed by "type"). + // - auxiliary data : + // pos : position in the CP tuple ( CP::size-1 ->0) + // ip : position in the tuple of indices (for sl) + // MP : accumulation of the final type metacomputed. + // + template struct pv_impl; + template struct pv_ : pv_impl,I...>{}; + + template struct pv_impl<-1, ip, CP, MP, I... > { + // the final type is a cartesian_product<...> if there is more than one mess + // and otherwise the only mesh remaining... (avoiding cartesian_product e.g. which makes little sense). + typedef typename std::conditional::type, MP>::type type; + template static T sl(T t, IT const & it) {return t;} + template static T m (T t, MT const & mt) {return t;} + }; + + template struct pv_impl { + typedef pv_impl B; + typedef typename B::type type; + template static auto sl (T t, IT const & it) DECL_AND_RETURN( B::sl(triqs::tuple::push_front(t,arrays::range()),it)); + template static auto m (T t, MT const & mt) DECL_AND_RETURN( B::m(t,mt)); + }; + + template struct pv_impl { + typedef typename cartesian_product_add_front::type, MP>::type MP2; + typedef pv_impl B; + typedef typename B::type type; + template static auto sl (T t, IT const & it) DECL_AND_RETURN( B::sl(triqs::tuple::push_front(t,std::get(it)),it)); + template static auto m (T t, MT const & mt) DECL_AND_RETURN( B::m (triqs::tuple::push_front(t,std::get(mt)),mt)); + }; + + // partial_eval<0> (g, 1) : returns : x -> g(1,x) + // partial_eval<1> (g, 3) : returns : x -> g(x,3) + // + template + gf_view,pos...>::type ,Target, Opt> + partial_eval(gf_impl< cartesian_product, Target,Opt,B> const & g, IT index) { + auto arr = reinterpret_linear_array(g.mesh(),g.data()); + typedef pv_,pos...> pv_t; + typedef gf_view< typename pv_t::type,Target, Opt> r_t; + auto comp = pv_t::m(std::make_tuple(),g.mesh().components()); + auto arr_args = pv_t::sl(std::make_tuple(),index); + // generalize this get<0> ---> flatten the tuple (construct from a tuple of args...) + return r_t{ std::get<0>(comp), triqs::tuple::apply(arr, arr_args), typename r_t::singularity_non_view_t{}, typename r_t::symmetry_t{} }; + } + + // to adapt the partial_eval as a polymorphic lambda (replace by a lambda in c++14) + template struct curry_polymorphic_lambda { + Gview g; + template + auto operator()(I ... i) const DECL_AND_RETURN( partial_eval(g,std::make_tuple(i...))); + }; + + // curry<0>(g) returns : x-> y... -> g(x,y...) + // curry<1>(g) returns : x-> y,z... -> g(y,x,z...) + // and so on + template + gf_view< typename pv_,pos...>::type, + lambda_valued, Target,Opt>,pos...>>, + Opt> + curry (gf_impl, Target,Opt,B> const & g) { + auto comp = pv_,pos...>::m(std::make_tuple(),g.mesh().components()); + typedef typename mesh< typename pv_,pos...>::type,Opt>::type m_t; + return {triqs::tuple::apply_construct(comp),curry_polymorphic_lambda, Target,Opt>, pos ...>{g}, nothing(), nothing()}; + }; + + } // gf_implementation + + using gfs_implementation::partial_eval; + using gfs_implementation::curry; + + /// ----- first implementation + /* + // slicing on first arg + template + gf_view slice_mesh0 (gf_view< cartesian_product, scalar_valued,Opt> g, size_t index) { + auto arr = reinterpret_linear_array(g.mesh(),g.data()); + typedef gf_view r_t; + return { std::get<1>(g.mesh().components()), arr(index,arrays::range()), typename r_t::singularity_non_view_t{}, typename r_t::symmetry_t{} }; + } + + // slicing on first arg + template + gf_view slice_mesh1 (gf_view< cartesian_product, scalar_valued,Opt> g, size_t index) { + auto arr = reinterpret_linear_array(g.mesh(),g.data()); + typedef gf_view r_t; + return { std::get<0>(g.mesh().components()), arr(arrays::range(), index), typename r_t::singularity_non_view_t{}, typename r_t::symmetry_t{} }; + } + + + template struct curry_lambda0 { + Gview g; + auto operator()(size_t i) const DECL_AND_RETURN( slice_mesh0(g,i)); + }; + + template + gf_view, scalar_valued,Opt>>>, Opt> + curry0 (gf_impl, scalar_valued,Opt,B> const & g) { + return {std::get<0>(g.mesh().components()),curry_lambda0, scalar_valued,Opt>>{g}, nothing(), nothing()}; + }; + +*/ +}} +#endif + diff --git a/triqs/gfs/data_proxies.hpp b/triqs/gfs/data_proxies.hpp index bfe2e52e..328914d9 100644 --- a/triqs/gfs/data_proxies.hpp +++ b/triqs/gfs/data_proxies.hpp @@ -42,9 +42,9 @@ namespace triqs { namespace gfs { arrays::matrix_view_proxy operator()(storage_view_t & data, size_t i) const { return arrays::matrix_view_proxy(data,i); } arrays::const_matrix_view_proxy operator()(storage_view_t const & data, size_t i) const { return arrays::const_matrix_view_proxy(data,i); } - template static void assign_no_resize (S & data, RHS && rhs) { data() = rhs;} - template static void assign_to_scalar (S & data, RHS && rhs) { data() = rhs;} - template static void assign_with_resize (storage_t & data, RHS && rhs) { data = rhs;} + template static void assign_no_resize (S & data, RHS && rhs) { data() = std::forward(rhs);} + template static void assign_to_scalar (S & data, RHS && rhs) { data() = std::forward(rhs);} + template static void assign_with_resize (storage_t & data, RHS && rhs) { data = std::forward(rhs);} template static void rebind (storage_view_t & data, RHS && rhs) { data.rebind(rhs.data()); } }; @@ -61,9 +61,9 @@ namespace triqs { namespace gfs { auto operator()(storage_view_t & data,size_t i) const -> decltype(data(i)) { return data(i);} auto operator()(storage_view_t const & data,size_t i) const -> decltype(data(i)) { return data(i);} - template static void assign_no_resize (S & data, RHS && rhs) { data() = rhs;} - template static void assign_to_scalar (S & data, RHS && rhs) { data() = rhs;} - template static void assign_with_resize (storage_t & data, RHS && rhs) { data = rhs;} + template static void assign_no_resize (S & data, RHS && rhs) { data() = std::forward(rhs);} + template static void assign_to_scalar (S & data, RHS && rhs) { data() = std::forward(rhs);} + template static void assign_with_resize (storage_t & data, RHS && rhs) { data = std::forward(rhs);} template static void rebind (storage_view_t & data, RHS && rhs) { data.rebind(rhs.data()); } }; @@ -92,6 +92,25 @@ namespace triqs { namespace gfs { template static void rebind (storage_view_t & data, RHS && rhs) { data.clear(); for (auto & x : rhs.data()) data.push_back(x);} }; + //---------------------------- vector ---------------------------------- + + template struct data_proxy_lambda { + + /// The storage + typedef F storage_t; + typedef F storage_view_t; + + /// The data access + auto operator()(storage_t & data, size_t i) DECL_AND_RETURN( data(i)); + auto operator()(storage_t const & data, size_t i) const DECL_AND_RETURN( data(i)); + + template static void assign_no_resize (S & data, RHS && rhs) { data() = std::forward(rhs);} + template static void assign_with_resize (S & data, RHS && rhs) = delete; + template static void assign_to_scalar (S & data, RHS && rhs) = delete; + template static void rebind (storage_view_t & data, RHS && rhs) = delete;// { data = std::forward(rhs);} + }; + + }} #endif diff --git a/triqs/gfs/gf.hpp b/triqs/gfs/gf.hpp index 820d1840..96f140aa 100644 --- a/triqs/gfs/gf.hpp +++ b/triqs/gfs/gf.hpp @@ -102,6 +102,9 @@ namespace triqs { namespace gfs { typedef gf_desc descriptor_t; + typedef Variable variable_t; + typedef Opt option_t; + typedef typename gfs_implementation::mesh::type mesh_t; typedef typename mesh_t::domain_t domain_t; typedef typename mesh_t::mesh_point_t mesh_point_t; diff --git a/triqs/gfs/meshes/product.hpp b/triqs/gfs/meshes/product.hpp index aa5d9fe1..5d77a885 100644 --- a/triqs/gfs/meshes/product.hpp +++ b/triqs/gfs/meshes/product.hpp @@ -48,7 +48,7 @@ namespace triqs { namespace gfs { /// Conversions point <-> index <-> linear_index struct _aux1 { template void operator()(P & p, M const & m, I const& i) {p = m.index_to_point(i);}}; - typename domain_t::point_t index_to_point(index_t const & ind) const { domain_pt_t res; triqs::tuple::apply_on_zip3(_aux1(), res,m_tuple,ind); return res;} + typename domain_t::point_t index_to_point(index_t const & ind) const { domain_pt_t res; triqs::tuple::apply_on_zip(_aux1(), res,m_tuple,ind); return res;} // index[0] + component[0].size * (index[1] + component[1].size* (index[2] + ....)) struct _aux2 { template size_t operator()(M const & m, I const & i,size_t R) {return m.index_to_linear(i) + R * m.size();}}; diff --git a/triqs/gfs/product.hpp b/triqs/gfs/product.hpp new file mode 100644 index 00000000..0601256f --- /dev/null +++ b/triqs/gfs/product.hpp @@ -0,0 +1,173 @@ +/******************************************************************************* + * + * TRIQS: a Toolbox for Research in Interacting Quantum Systems + * + * Copyright (C) 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 . + * + ******************************************************************************/ +#ifndef TRIQS_GF_PRODUCT_H +#define TRIQS_GF_PRODUCT_H +#include "./tools.hpp" +#include "./gf.hpp" +#include "./meshes/product.hpp" + +namespace triqs { namespace gfs { + + template struct cartesian_product{ + typedef std::tuple type; + static constexpr size_t size = sizeof...(Ms); + }; + + namespace gfs_implementation { + // the mesh is simply a cartesian product + template struct mesh,Opt> { + typedef mesh_product< typename mesh::type ... > type; + typedef std::tuple mesh_name_t; + static type make (typename mesh::type ... ms) { return type{std::move(ms)...};} + }; + + // h5 name : name1_x_name2_..... + template struct h5_name,matrix_valued,Opt> { + static std::string invoke(){ + return triqs::tuple::fold( + [](std::string a, std::string b) { return a + std::string(b.empty()?"" : "_x_") + b;}, + std::make_tuple(h5_name::invoke()...), + std::string()); + } + }; + + /// --------------------------- data access --------------------------------- + + template struct data_proxy,matrix_valued,Opt> : data_proxy_array,3> {}; + template struct data_proxy,scalar_valued,Opt> : data_proxy_array,1> {}; + + /// --------------------------- evaluator --------------------------------- + + struct evaluator_grid_simple { + size_t n; + evaluator_grid_simple() = default; + + template + evaluator_grid_simple (MeshType const & m, PointType const & p) { n=p; } + template auto operator()(F const & f) const DECL_AND_RETURN(f (n)); + }; + + struct evaluator_grid_linear_interpolation { + double w1, w2; size_t n1, n2; + + evaluator_grid_linear_interpolation() = default; + + template + evaluator_grid_linear_interpolation (MeshType const & m, PointType const & p, double prefactor=1) { // delegate ! + bool in; double w; + std::tie(in, n1, w) = windowing(m,p); + //std::cout << in << " "<< n1 << " "<< w << " " << p << std::endl; + if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds"; + w1 = prefactor * w; w2 = prefactor *(1-w); n2 = n1 +1; + } + + template auto operator()(F const & f) const DECL_AND_RETURN(w1 * f(n1) + w2 * f (n2)); + }; + + template struct evaluator_fnt_on_mesh; + + // can not use inherited constructors, too recent... +#define TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(NEWCLASS,CLASS) : CLASS { template NEWCLASS(T &&... t) : CLASS(std::forward(t)...){};}; + + template<> struct evaluator_fnt_on_mesh TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_simple); + template<> struct evaluator_fnt_on_mesh TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_linear_interpolation); + template<> struct evaluator_fnt_on_mesh TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_linear_interpolation); + template<> struct evaluator_fnt_on_mesh TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_linear_interpolation); + + + /** + * This the multi-dimensional evaluator. + * It combine the evaluator of each components, as long as they are a linear form + * eval(g, x) = \sum_i w_i g( n_i(x)) , with w some weight and n_i some points on the grid. + * Mathematically, it is written as (example of evaluating g(x1,x2,x3,x4)). + * Notation : eval(X) : g -> g(X) + * eval(x1,x2,x3,x4) (g) = eval (x1) ( binder ( g, (), (x2,x3,x4)) ) + * binder( g, (), (x2,x3,x4)) (p1) = eval(x2)(binder (g,(p1),(x3,x4))) + * binder( g, (p1), (x3,x4)) (p2) = eval(x3)(binder (g,(p1,p2),(x4))) + * binder( g, (p1,p2), (x4)) (p3) = eval(x4)(binder (g,(p1,p2,p3),())) + * binder( g, (p1,p2,p3),()) (p4) = g[p1,p2,p3,p4] + * + * p_i are points on the grids, x_i points in the domain. + * + * Unrolling the formula gives (for 2 variables, with 2 points interpolation) + * eval(xa,xb) (g) = eval (xa) ( binder ( g, (), (xb)) ) = w_1(xa) binder ( g, (), (xb))( n_1(xa)) + w_2(xa) binder ( g, (), (xb))( n_2(xa)) + * = w_1(xa) ( eval(xb)( binder ( g, (n_1(xa) ), ()))) + 1 <-> 2 + * = w_1(xa) ( W_1(xb) * binder ( g, (n_1(xa) ), ())(N_1(xb)) + 1<->2 ) + 1 <-> 2 + * = w_1(xa) ( W_1(xb) * g[n_1(xa), N_1(xb)] + 1<->2 ) + 1 <-> 2 + * = w_1(xa) ( W_1(xb) * g[n_1(xa), N_1(xb)] + W_2(xb) * g[n_1(xa), N_2(xb)] ) + 1 <-> 2 + * which is the expected formula + */ + // implementation : G = gf, Tn : tuple of n points, Ev : tuple of evaluators (the evals functions), pos = counter from #args-1 =>0 + // NB : the tuple is build in reverse with respect to the previous comment. + template struct binder; + + template + binder make_binder(G const * g, Tn tn, Ev const & ev) { return binder{g, std::move(tn), ev}; } + + template struct binder { + G const * g; Tn tn; Ev const & evals; + auto operator()(size_t p) const DECL_AND_RETURN( std::get(evals) ( make_binder(g, triqs::tuple::push_front(tn,p), evals) )); + }; + + template struct binder { + G const * g; Tn tn; Ev const & evals; + auto operator()(size_t p) const DECL_AND_RETURN( triqs::tuple::apply(on_mesh(*g), triqs::tuple::push_front(tn,p))); + }; + + // now the multi d evaluator itself. + template + struct evaluator,Target,Opt> { + static constexpr int arity = sizeof...(Ms); + mutable std::tuple< evaluator_fnt_on_mesh ... > evals; + + struct _poly_lambda {// replace by a polymorphic lambda in C++14 + template void operator()(A & a, B const & b, C const & c) const { a = A{b,c};} + }; + + template + std::complex operator() (G const * g, Args && ... args) const { + static constexpr int R = sizeof...(Args); + // build the evaluators, as a tuple of ( evaluator ( mesh_component, args)) + triqs::tuple::call_on_zip(_poly_lambda(), evals, g->mesh().components(), std::make_tuple(args...)); + return std::get(evals) (make_binder (g, std::make_tuple(), evals) ); + } + }; + + // ------------------------------- Factories -------------------------------------------------- + + template + struct factories, scalar_valued,Opt> { + typedef gf, scalar_valued,Opt> gf_t; + + template + static gf_t make_gf(Meshes && ... meshes) { + auto m = make_gf_mesh,Opt>(meshes...); + typename gf_t::data_non_view_t A(m.size()); + A() =0; + return gf_t (m, std::move(A), nothing(), nothing()); + } + }; + + } // gf_implementation + +}} +#endif + diff --git a/triqs/gfs/refreq.hpp b/triqs/gfs/refreq.hpp index b94f47bd..197f127a 100644 --- a/triqs/gfs/refreq.hpp +++ b/triqs/gfs/refreq.hpp @@ -54,6 +54,7 @@ namespace triqs { namespace gfs { typedef typename std::conditional < std::is_same::value, arrays::matrix >, std::complex>::type rtype; template rtype operator() (G const * g,double w0) const { + //auto operator() (G const * g,double w0) const -> typename decltype ((*g)[0])::non_view_type { size_t n; double w; bool in; std::tie(in, n, w) = windowing(g->mesh(),w0); if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds"; diff --git a/triqs/gfs/refreq_imtime.hpp b/triqs/gfs/refreq_imtime.hpp index ea7010d5..796baebd 100644 --- a/triqs/gfs/refreq_imtime.hpp +++ b/triqs/gfs/refreq_imtime.hpp @@ -69,6 +69,9 @@ namespace triqs { namespace gfs { if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds"; auto gg = on_mesh(*g); //[g]( size_t n1, size_t n2) {return g->on_mesh(n1,n2);}; auto res = w1 *( w2*gg(n1,n2) + (1-w2)*gg(n1,n2+1)) + (1-w1) * ( w2*gg(n1+1,n2) + (1-w2)*gg(n1+1,n2+1)); + + //std::cout << "eval reref imtim"<< n1 << " "<< n2 << " "<< w1 << " " << w2 << " "<< omega << " "<< tau<< std::endl; + return ((std::get<1>(g->mesh().components()).domain().statistic == Fermion) && (p%2==1) ? -res : res); } };