gf product and curry. First draft

- works on simple case, to be reread.
- curry in general. To be reread ..
- added some tests.

test/triqs/gf/gf_re_im_freq_time.cpp

@@ -3,6 +3,8 @@
 #include <triqs/gfs/re_im_freq.hpp> 
 #include <triqs/gfs/re_im_time.hpp> 
 #include <triqs/gfs/refreq_imtime.hpp> 
+#include <triqs/gfs/product.hpp> 
+#include <triqs/gfs/curry.hpp> 
 
 #include <triqs/gfs/local/fourier_real.hpp> 
 #include <triqs/arrays.hpp>
@@ -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<re_im_time, scalar_valued>(   tmin, tmax, n_re_time, beta, Fermion, n_im_time);
+ //auto G_t_tau= make_gf<re_im_time, scalar_valued>(   tmin, tmax, n_re_time, beta, Fermion, n_im_time);
  auto G_w_wn = make_gf<re_im_freq, scalar_valued>(   wmin, wmax, n_re_freq, beta, Fermion, n_im_freq);
  auto G_w_tau= make_gf<refreq_imtime, scalar_valued>(wmin, wmax, n_re_freq, beta, Fermion, n_im_time);
  auto G_w= make_gf<refreq, scalar_valued>(wmin, wmax, n_re_freq);
- 
+
+ auto G_t_tau= make_gf<cartesian_product<retime,imtime>, scalar_valued>(  make_gf_mesh<retime>(tmin, tmax, n_re_time), make_gf_mesh<imtime>(beta, Fermion, n_im_time));
+ //auto G_t_tau_N= make_gf<cartesian_product<retime,imtime>, scalar_valued>(  {tmin, tmax, n_re_time}, {beta, Fermion, n_im_time});
+
+ auto G_w_wn2 = make_gf<cartesian_product<refreq,imfreq>, scalar_valued>( make_gf_mesh<refreq>(wmin, wmax, n_re_freq), make_gf_mesh<imfreq>(beta, Fermion, n_im_freq));
+ auto G_w_tau2 = make_gf<cartesian_product<refreq,imtime>, scalar_valued>( make_gf_mesh<refreq>(wmin, wmax, n_re_freq), make_gf_mesh<imtime>(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<double>{{0}};
  G_w.singularity()(1)=triqs::arrays::matrix<double>{{0}};
  G_w.singularity()(2)=triqs::arrays::matrix<double>{{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_wn( 0.789,0.123) "<< 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_wn2( 0.789,0.123) << std::endl; 
- 
+ std::cout << "G_w_tau(0.789,0.123)" << G_w_tau(0.789,0.123) << std::endl; 
+ std::cout << "G_w_tau(0.789,0.123)" << G_w_tau2(0.789,0.123) << std::endl; 
+
+ // test curry
+ std::cout  << "curry no"<< G_w_wn.on_mesh(8,3) << std::endl ;
+
+ auto G_w_wn_curry0 = curry<0>(G_w_wn2);
+ auto G_w_wn_curry1 = curry<1>(G_w_wn2);
+ auto G_w_wn2_view2 = G_w_wn2();
+
+ std::cout  << " curry "<<G_w_wn_curry0[8] << G_w_wn_curry0[8][3] << G_w_wn2_view2.on_mesh(8,3) << G_w_wn2_view2.on_mesh(8,3) <<std::endl ;
+ std::cout  << " curry "<<G_w_wn_curry1[3][8] << std::endl; 
+ std::cout  << "G_w_wn_sl0_a [3]"<<G_w_wn_sl0_a[3] << std::endl ;
+
+
  // test hdf5 
  H5::H5File file("gf_re_im_freq_time.h5", H5F_ACC_TRUNC );
  h5_write(file, "g_t_tau", G_t_tau);
  h5_write(file, "g_w_wn",  G_w_wn);
+ h5_write(file, "g_w_wn2", G_w_wn2);
  h5_write(file, "g_w_tau", G_w_tau);
 
+ /*
  // try to slice it
  auto gt = slice_mesh_imtime(G_t_tau, 1);
  h5_write(file, "gt0", gt);
@@ -71,5 +111,7 @@ int main() {
  auto G_t2 = inverse_fourier(slice_mesh_imtime(G_w_tau, index) );
  for(auto& t:G_t.mesh())
   if ( std::abs(G_t[t]-G_t2[t]) > 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;}
 }

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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef TRIQS_GF_CURRY_H
+#define TRIQS_GF_CURRY_H
+#include "./product.hpp"
+
+namespace triqs { namespace gfs { 
+
+ template<typename F> struct lambda_valued {};
+
+ namespace gfs_implementation { 
+
+  /// ---------------------------  data access  ---------------------------------
+
+  template<typename Opt, typename F, typename M> struct data_proxy<M,lambda_valued<F>,Opt> : data_proxy_lambda<F> {};
+
+  /// ---------------------------  Factories ---------------------------------
+
+  template<typename F, typename Opt, typename ... Ms>
+   struct factories<cartesian_product<Ms...>, lambda_valued<F>, Opt> {};
+
+  // detail
+  template<typename M0, typename CP> struct cartesian_product_add_front;
+  template<typename M0, typename ... Ms> struct cartesian_product_add_front<M0,cartesian_product<Ms...>>{ typedef cartesian_product<M0,Ms...> 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<int pos, int ip, typename CP, typename MP, int ... I> struct pv_impl; 
+  template<typename CP, int ... I> struct pv_ : pv_impl<CP::size-1,0,CP,cartesian_product<>,I...>{}; 
+ 
+  template<typename CP, int ip, typename MP, int ... I>  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<imfreq> e.g. which makes little sense).
+   typedef typename std::conditional<MP::size==1, typename std::tuple_element<0,typename MP::type>::type, MP>::type type; 
+   template<typename T, typename IT> static T sl(T t, IT const & it) {return t;}
+   template<typename T, typename MT> static T m (T t, MT const & mt) {return t;}
+  }; 
+
+  template<int pos, int ip, typename CP, typename MP, int ... I> struct pv_impl {
+   typedef pv_impl<pos-1, ip, CP,MP, I...> B;
+   typedef typename B::type type;
+   template<typename T, typename IT> static auto sl (T t, IT const & it) DECL_AND_RETURN( B::sl(triqs::tuple::push_front(t,arrays::range()),it));
+   template<typename T, typename MT> static auto m  (T t, MT const & mt) DECL_AND_RETURN( B::m(t,mt));
+  };
+
+  template<int pos, int ip, typename CP, typename MP, int ... I> struct pv_impl<pos, ip, CP, MP, pos ,I...> {
+   typedef typename cartesian_product_add_front<typename  std::tuple_element<pos,typename CP::type>::type, MP>::type MP2; 
+   typedef pv_impl<pos-1,ip+1, CP,MP2,I...> B;
+   typedef typename B::type type;
+   template<typename T, typename IT> static auto sl (T t, IT const & it) DECL_AND_RETURN( B::sl(triqs::tuple::push_front(t,std::get<ip >(it)),it));
+   template<typename T, typename MT> static auto m  (T t, MT const & mt) DECL_AND_RETURN( B::m (triqs::tuple::push_front(t,std::get<pos>(mt)),mt));
+  };
+
+  // partial_eval<0> (g, 1) : returns :  x -> g(1,x)
+  // partial_eval<1> (g, 3) : returns :  x -> g(x,3)
+  //
+  template<int ... pos, typename Opt, typename Target, bool B, typename IT, typename ... Ms>
+   gf_view<typename pv_<cartesian_product<Ms...>,pos...>::type ,Target, Opt> 
+   partial_eval(gf_impl< cartesian_product<Ms...>, Target,Opt,B> const & g, IT index) { 
+    auto arr = reinterpret_linear_array(g.mesh(),g.data());
+    typedef pv_<cartesian_product<Ms...>,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<typename Gview, int ... pos> struct curry_polymorphic_lambda { 
+   Gview g;
+   template<typename ...I>
+    auto operator()(I ... i) const DECL_AND_RETURN( partial_eval<pos...>(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<int ... pos, typename Target, typename Opt, bool B, typename ... Ms>
+   gf_view< typename pv_<cartesian_product<Ms...>,pos...>::type,
+	   lambda_valued<curry_polymorphic_lambda<gf_view<cartesian_product<Ms...>, Target,Opt>,pos...>>, 
+	   Opt>
+    curry (gf_impl<cartesian_product<Ms...>, Target,Opt,B> const & g) { 
+    auto comp =  pv_<cartesian_product<Ms...>,pos...>::m(std::make_tuple(),g.mesh().components());
+    typedef typename mesh< typename pv_<cartesian_product<Ms...>,pos...>::type,Opt>::type m_t;
+    return {triqs::tuple::apply_construct<m_t>(comp),curry_polymorphic_lambda<gf_view<cartesian_product<Ms...>, Target,Opt>, pos ...>{g}, nothing(), nothing()}; 
+   };
+
+ } // gf_implementation
+
+ using gfs_implementation::partial_eval;
+ using gfs_implementation::curry;
+
+ /// ----- first implementation
+ /*
+ // slicing on first arg
+ template<typename Opt, typename M0,  typename M1>
+ gf_view<M1,scalar_valued, Opt> slice_mesh0 (gf_view< cartesian_product<M0,M1>, scalar_valued,Opt> g, size_t index) { 
+ auto arr = reinterpret_linear_array(g.mesh(),g.data());
+ typedef gf_view<M1,scalar_valued, Opt> 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<typename Opt, typename M0, typename M1>
+ gf_view<M0,scalar_valued, Opt> slice_mesh1 (gf_view< cartesian_product<M0,M1>, scalar_valued,Opt> g, size_t index) { 
+ auto arr = reinterpret_linear_array(g.mesh(),g.data()); 
+ typedef gf_view<M0,scalar_valued, Opt> 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<typename Gview> struct curry_lambda0 { 
+ Gview g;
+ auto operator()(size_t i) const DECL_AND_RETURN( slice_mesh0(g,i));
+ };
+
+ template<typename Opt, bool B, typename M0, typename M1>
+ gf_view<M0,lambda_valued<curry_lambda0<gf_view<cartesian_product<M0,M1>, scalar_valued,Opt>>>, Opt> 
+ curry0 (gf_impl<cartesian_product<M0,M1>, scalar_valued,Opt,B> const & g) { 
+ return {std::get<0>(g.mesh().components()),curry_lambda0<gf_view<cartesian_product<M0,M1>, scalar_valued,Opt>>{g}, nothing(), nothing()}; 
+ };
+
+*/
+}}
+#endif
+

triqs/gfs/data_proxies.hpp

@@ -42,9 +42,9 @@ namespace triqs { namespace gfs {
   arrays::matrix_view_proxy<storage_view_t,0>       operator()(storage_view_t       & data, size_t i) const { return arrays::matrix_view_proxy<storage_view_t,0>(data,i); } 
   arrays::const_matrix_view_proxy<storage_view_t,0> operator()(storage_view_t const & data, size_t i) const { return arrays::const_matrix_view_proxy<storage_view_t,0>(data,i); } 
 
-  template<typename S, typename RHS> static void assign_no_resize (S & data, RHS && rhs)           { data() = rhs;}
+  template<typename S, typename RHS> static void assign_no_resize (S & data, RHS && rhs)           { data() = std::forward<RHS>(rhs);}
-  template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs)           { data() = rhs;}
+  template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs)           { data() = std::forward<RHS>(rhs);}
-  template<typename RHS>             static void assign_with_resize (storage_t & data, RHS && rhs) { data = rhs;}
+  template<typename RHS>             static void assign_with_resize (storage_t & data, RHS && rhs) { data = std::forward<RHS>(rhs);}
   template<typename RHS>             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<typename S, typename RHS> static void assign_no_resize (S & data, RHS && rhs)           { data() = rhs;}
+  template<typename S, typename RHS> static void assign_no_resize (S & data, RHS && rhs)           { data() = std::forward<RHS>(rhs);}
-  template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs)           { data() = rhs;}
+  template<typename S, typename RHS> static void assign_to_scalar (S & data, RHS && rhs)           { data() = std::forward<RHS>(rhs);}
-  template<typename RHS>             static void assign_with_resize (storage_t & data, RHS && rhs) { data = rhs;}
+  template<typename RHS>             static void assign_with_resize (storage_t & data, RHS && rhs) { data = std::forward<RHS>(rhs);}
   template<typename RHS>             static void rebind (storage_view_t & data, RHS && rhs)        { data.rebind(rhs.data()); }
  };
 
@@ -92,6 +92,25 @@ namespace triqs { namespace gfs {
   template<typename RHS> static void rebind (storage_view_t & data, RHS && rhs) { data.clear(); for (auto & x : rhs.data()) data.push_back(x);}
  };
 
+ //---------------------------- vector ----------------------------------
+
+ template<typename F> 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<typename S, typename RHS> static void assign_no_resize (S & data, RHS && rhs) { data() = std::forward<RHS>(rhs);}
+  template<typename S, typename RHS> static void assign_with_resize (S & data, RHS && rhs) = delete;
+  template<typename S, typename RHS> static void assign_to_scalar   (S & data, RHS && rhs) = delete;
+  template<typename RHS> static void rebind (storage_view_t & data, RHS && rhs) = delete;// { data = std::forward<RHS>(rhs);}
+ };
+
+
 }}
 #endif
 

triqs/gfs/gf.hpp

@@ -102,6 +102,9 @@ namespace triqs { namespace gfs {
 
     typedef gf_desc<Variable,Target,Opt>    descriptor_t;
 
+    typedef Variable variable_t;
+    typedef Opt option_t;
+
     typedef typename gfs_implementation::mesh<Variable,Opt>::type                                mesh_t;
     typedef typename mesh_t::domain_t                                                           domain_t;
     typedef typename mesh_t::mesh_point_t                                                       mesh_point_t;

triqs/gfs/meshes/product.hpp

@@ -48,7 +48,7 @@ namespace triqs { namespace gfs {
 
    /// Conversions point <-> index <-> linear_index
    struct _aux1 { template<typename P, typename M, typename I> 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<typename I, typename M> size_t operator()(M const & m, I const & i,size_t R) {return m.index_to_linear(i) + R * m.size();}};

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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#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<typename ... Ms> struct cartesian_product{ 
+  typedef std::tuple<Ms...> type; 
+  static constexpr size_t size = sizeof...(Ms); 
+ };
+
+ namespace gfs_implementation { 
+  // the mesh is simply a cartesian product
+  template<typename Opt, typename ... Ms> struct mesh<cartesian_product<Ms...>,Opt>  { 
+   typedef mesh_product< typename mesh<Ms,Opt>::type ... > type;
+   typedef std::tuple<Ms...> mesh_name_t;
+   static type make (typename mesh<Ms,Opt>::type ... ms) { return type{std::move(ms)...};} 
+  };
+
+  // h5 name : name1_x_name2_.....
+  template<typename Opt, typename ... Ms> struct h5_name<cartesian_product<Ms...>,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<Ms,matrix_valued,Opt>::invoke()...), 
+      std::string());
+   }
+  };
+
+  /// ---------------------------  data access  ---------------------------------
+
+  template<typename Opt, typename ... Ms> struct data_proxy<cartesian_product<Ms...>,matrix_valued,Opt> : data_proxy_array<std::complex<double>,3> {};
+  template<typename Opt, typename ... Ms> struct data_proxy<cartesian_product<Ms...>,scalar_valued,Opt> : data_proxy_array<std::complex<double>,1> {};
+
+  /// ---------------------------  evaluator ---------------------------------
+
+  struct evaluator_grid_simple { 
+   size_t n; 
+   evaluator_grid_simple() = default;
+
+   template<typename MeshType, typename PointType> 
+    evaluator_grid_simple (MeshType const & m, PointType const & p) { n=p; }
+   template<typename F> 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<typename MeshType, typename PointType> 
+    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<typename F> auto operator()(F const & f) const DECL_AND_RETURN(w1 * f(n1) + w2 * f (n2));
+  };
+
+  template<typename MeshType> struct evaluator_fnt_on_mesh;
+
+  // can not use inherited constructors, too recent...
+#define TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(NEWCLASS,CLASS) : CLASS { template<typename ...T> NEWCLASS(T &&... t) : CLASS(std::forward<T>(t)...){};}; 
+
+  template<> struct evaluator_fnt_on_mesh<imfreq> TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_simple);
+  template<> struct evaluator_fnt_on_mesh<imtime> TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_linear_interpolation);
+  template<> struct evaluator_fnt_on_mesh<retime> TRIQS_INHERIT_AND_FORWARD_CONSTRUCTOR(evaluator_fnt_on_mesh, evaluator_grid_linear_interpolation);
+  template<> struct evaluator_fnt_on_mesh<refreq> 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<typename G, typename Tn, typename Ev, int pos> struct binder;
+
+  template<int pos, typename G, typename Tn, typename Ev> 
+   binder<G,Tn,Ev,pos> make_binder(G const * g, Tn tn, Ev const & ev) { return binder<G,Tn,Ev,pos>{g, std::move(tn), ev}; }
+
+  template<typename G, typename Tn, typename Ev, int pos> struct binder {
+   G const * g; Tn tn; Ev const & evals;
+   auto operator()(size_t p) const DECL_AND_RETURN( std::get<pos>(evals) ( make_binder<pos-1>(g, triqs::tuple::push_front(tn,p), evals) ));
+  };
+
+  template<typename G, typename Tn, typename Ev> struct binder<G,Tn,Ev,-1> {
+   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<typename Target, typename Opt, typename ... Ms>
+   struct evaluator<cartesian_product<Ms...>,Target,Opt> {
+    static constexpr int arity = sizeof...(Ms);
+    mutable std::tuple< evaluator_fnt_on_mesh<Ms> ... > evals;
+
+    struct _poly_lambda {// replace by a polymorphic lambda in C++14 
+     template<typename A, typename B, typename C> void operator()(A & a, B const & b, C const & c) const { a = A{b,c};}
+    };
+
+    template<typename G, typename ... Args>
+     std::complex<double> operator() (G const * g, Args && ... args)  const {
+      static constexpr int R = sizeof...(Args);
+      // build the evaluators, as a tuple of ( evaluator<Ms> ( mesh_component, args))
+      triqs::tuple::call_on_zip(_poly_lambda(), evals, g->mesh().components(), std::make_tuple(args...));
+      return std::get<R-1>(evals) (make_binder<R-2> (g, std::make_tuple(), evals) );
+     } 
+   };
+
+  // -------------------------------   Factories  --------------------------------------------------
+
+  template<typename Opt, typename ... Ms>
+   struct factories<cartesian_product<Ms...>, scalar_valued,Opt> {
+    typedef gf<cartesian_product<Ms...>, scalar_valued,Opt> gf_t;
+
+    template<typename ... Meshes>
+     static gf_t make_gf(Meshes && ... meshes) { 
+      auto m =  make_gf_mesh<cartesian_product<Ms...>,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
+

triqs/gfs/refreq.hpp

@@ -54,6 +54,7 @@ namespace triqs { namespace gfs {
     typedef typename std::conditional < std::is_same<Target, matrix_valued>::value, arrays::matrix<std::complex<double> >, std::complex<double>>::type rtype; 
     template<typename G>
      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";

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);
       } 
     };