gf/clef. Bug fix, and added scalar->matrix reinterpret for blocks

- clef : fix a little bug in storage when evaluating
  (was using the wrong trait to deduce storage type).
- gf : block :
     - added reinterpret_scalar_valued_gf_as_matrix_valued
       for block function
     - cleaned make_block_gf_view_from_vector
     - added make_block_gf_view_from_vectormake_block_gf_view_from_vector_of_cython_proxy
       and changed the cython accordingly because it requires a slightly different syntax.
     - updated tests
- gf : cleaned some template.

pytriqs/gf/local/imfreq.pxd

@@ -48,7 +48,7 @@ cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
         gf_imfreq & operator [](int)
         discrete_mesh & mesh()
 
-    cdef gf_block_imfreq  make_gf_block_imfreq "triqs::gfs::make_block_gf_view_from_vector<triqs::gfs::gf<triqs::gfs::imfreq>>" (vector[gf_imfreq] &) 
+    cdef gf_block_imfreq  make_gf_block_imfreq "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::imfreq>>" (vector[gf_imfreq] &) 
     #cdef gf_block_imfreq  make_gf_block_imfreq "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::imfreq>>" (  vector[gf_imfreq] &) 
 
 cdef gf_block_imfreq  as_gf_block_imfreq (G) except +

pytriqs/gf/local/imtime.pxd

@@ -49,7 +49,7 @@ cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
         gf_imtime & operator [](int)
         discrete_mesh & mesh()
 
-    cdef gf_block_imtime  make_gf_block_imtime "triqs::gfs::make_block_gf_view_from_vector<triqs::gfs::gf<triqs::gfs::imtime>>" (vector[gf_imtime] &) 
+    cdef gf_block_imtime  make_gf_block_imtime "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::imtime>>" (vector[gf_imtime] &) 
     #cdef gf_block_imtime  make_gf_block_imtime "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::imtime>>" (  vector[gf_imtime] &) 
 
 cdef gf_block_imtime  as_gf_block_imtime (G) except +

pytriqs/gf/local/legendre.pxd

@@ -46,7 +46,7 @@ cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
         gf_legendre & operator [](int)
         discrete_mesh & mesh()
 
-    cdef gf_block_legendre  make_gf_block_legendre "triqs::gfs::make_block_gf_view_from_vector<triqs::gfs::gf<triqs::gfs::legendre>>" (vector[gf_legendre] &) 
+    cdef gf_block_legendre  make_gf_block_legendre "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::legendre>>" (vector[gf_legendre] &) 
     #cdef gf_block_legendre  make_gf_block_legendre "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::legendre>>" (  vector[gf_legendre] &) 
 
 cdef gf_block_legendre  as_gf_block_legendre (G) except +

pytriqs/gf/local/refreq.pxd

@@ -48,7 +48,7 @@ cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
         gf_refreq & operator [](int)
         discrete_mesh & mesh()
 
-    cdef gf_block_refreq  make_gf_block_refreq "triqs::gfs::make_block_gf_view_from_vector<triqs::gfs::gf<triqs::gfs::refreq>>" (vector[gf_refreq] &) 
+    cdef gf_block_refreq  make_gf_block_refreq "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::refreq>>" (vector[gf_refreq] &) 
     #cdef gf_block_refreq  make_gf_block_refreq "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::refreq>>" (  vector[gf_refreq] &) 
 
 cdef gf_block_refreq  as_gf_block_refreq (G) except +

pytriqs/gf/local/retime.pxd

@@ -48,7 +48,7 @@ cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
         gf_retime & operator [](int)
         discrete_mesh & mesh()
 
-    cdef gf_block_retime  make_gf_block_retime "triqs::gfs::make_block_gf_view_from_vector<triqs::gfs::gf<triqs::gfs::retime>>" (vector[gf_retime] &) 
+    cdef gf_block_retime  make_gf_block_retime "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::retime>>" (vector[gf_retime] &) 
     #cdef gf_block_retime  make_gf_block_retime "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::retime>>" (  vector[gf_retime] &) 
 
 cdef gf_block_retime  as_gf_block_retime (G) except +

pytriqs/gf/local/two_real_times.pxd

@@ -51,7 +51,7 @@ cdef extern from "triqs/gfs/block.hpp" namespace "triqs::gfs" :
         gf_two_real_times & operator [](int)
         discrete_mesh & mesh()
 
-    cdef gf_block_two_real_times  make_gf_block_two_real_times "triqs::gfs::make_block_gf_view_from_vector<triqs::gfs::gf<triqs::gfs::two_real_times>>" (  vector[gf_two_real_times] &) 
+    cdef gf_block_two_real_times  make_gf_block_two_real_times "triqs::gfs::make_block_gf_view_from_vector_of_cython_proxy<triqs::gfs::gf<triqs::gfs::two_real_times>>" (  vector[gf_two_real_times] &) 
     #cdef gf_block_two_real_times  make_gf_block_two_real_times "triqs::gfs::make_gf_view<triqs::gfs::block_index,triqs::gfs::gf<triqs::gfs::two_real_times>>" (  vector[gf_two_real_times] &) 
 
 cdef gf_block_two_real_times  as_gf_block_two_real_times (G) except +

test/triqs/gfs/block.cpp

@@ -1,75 +1,88 @@
 #define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
-#include <triqs/gfs.hpp> 
+#include <triqs/gfs.hpp>
 using namespace triqs::gfs;
 using namespace triqs;
-#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl;
+#define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> " << (X) << std::endl << std::endl;
 
 int main() {
  try {
-  double beta =1;
+  double beta = 1;
-  auto G1 = gf<imfreq> ({beta, Fermion}, {2,2});
+  auto G1 = gf<imfreq>({beta, Fermion}, {2, 2});
   auto G2 = G1;
   auto G3 = G2;
 
   // construct some block functions
-  auto B0 = block_gf<imfreq> (3);
+  auto B0 = block_gf<imfreq>(3);
 
-  auto B1 = make_block_gf<imfreq> (3, G1); 
+  auto B1 = make_block_gf<imfreq>(3, G1);
-  auto B2 = make_block_gf<imfreq> ({G1,G1,G1});
+  auto B2 = make_block_gf<imfreq>({G1, G1, G1});
-  auto B3 = make_block_gf<imfreq> ({"a","b","c"}, {G1,G1,G1}); 
+  auto B3 = make_block_gf<imfreq>({"a", "b", "c"}, {G1, G1, G1});
-  auto B4 = block_gf<imfreq> (1);
+  auto B4 = block_gf<imfreq>(1);
 
-  // test hdf5 
+  // test hdf5
   {
-    H5::H5File file("ess_gf.h5", H5F_ACC_TRUNC );
+   H5::H5File file("ess_gf.h5", H5F_ACC_TRUNC);
-    h5_write(file, "B3", B3);
+   h5_write(file, "B3", B3);
   }
 
   {
-    H5::H5File file("ess_gf.h5", H5F_ACC_RDONLY);
+   H5::H5File file("ess_gf.h5", H5F_ACC_RDONLY);
-    std::cout  << "B4 mesh" << B4.mesh().size()<<std::endl;
+   std::cout << "B4 mesh" << B4.mesh().size() << std::endl;
-    h5_read(file, "B3", B4);
+   h5_read(file, "B3", B4);
-    std::cout  << "B4 mesh" << B4.mesh().size()<<std::endl;
+   std::cout << "B4 mesh" << B4.mesh().size() << std::endl;
   }
 
   B1[0][0] = 98;
-  //not implemented yet
+  // not implemented yet
-  //B3["a"][0] = 98;
+  // B3["a"][0] = 98;
 
-  auto View =  make_block_gf_view(G1,G2,G3);
+  auto View = make_block_gf_view(G1, G2, G3);
 
-  std::cout  << "Number of blocks " << View.mesh().size()<<std::endl ;
+  std::cout << "Number of blocks " << View.mesh().size() << std::endl;
   auto g0 = View[0];
   auto g0v = View[0]();
 
   auto Gv = g0();
 
   Gv[0] = 20;
-  TEST( G1( 0) ) ;
+  TEST(G1(0));
   Gv[0] = 0;
 
   g0v[0] = 3.2;
-  TEST( G1( 0) ) ;
+  TEST(G1(0));
 
   // Operation
   g0[0] = 3.2;
-  TEST( View[0](0) ) ;
+  TEST(View[0](0));
-  View = View/2; 
+  View = View / 2;
-  TEST( View[0](0) ) ;
+  TEST(View[0](0));
 
   // try the loop over the block.
-  for (auto & g : View) { g[0] = 20;}
+  for (auto& g : View) {
-  for (auto & g : B1) { g[0] = 20;}
+   g[0] = 20;
+  }
+  for (auto& g : B1) {
+   g[0] = 20;
+  }
 
- // check chaining of clef
+  // check chaining of clef
   clef::placeholder<0> b_;
   clef::placeholder<1> om_;
-  B1[b_][om_] << b_ / ( om_ + 2);
+  B1[b_][om_] << b_ / (om_ + 2);
 
- // does not work
+  auto B11 = B1;
- //std::vector<gf<imfreq>> green_v;
- //green_v.emplace_back({beta, Fermion}, {2,2} ); 
 
+  B1[b_](om_) << B11[b_](om_) * B1[b_](om_) * B11[b_](om_);
+
+  // does not work
+  // std::vector<gf<imfreq>> green_v;
+  // green_v.emplace_back({beta, Fermion}, {2,2} );
+
+  // test reinterpretation
+  // compile only, add more test here
+  auto gs1 = gf<imfreq, scalar_valued>({beta, Fermion});
+  auto bgs = make_block_gf<imfreq>(3, gs1);
+  auto bg = reinterpret_scalar_valued_gf_as_matrix_valued(bgs);
  }
-TRIQS_CATCH_AND_ABORT;
+ TRIQS_CATCH_AND_ABORT;
 }

test/triqs/gfs/block.output

@@ -1,11 +1,11 @@
 B4 mesh1
 B4 mesh3
 Number of blocks 3
-(G1( 0)) ---> 
+(G1(0)) ---> 
 [[(20,0),(0,0)]
  [(0,0),(20,0)]]
 
-(G1( 0)) ---> 
+(G1(0)) ---> 
 [[(3.2,0),(0,0)]
  [(0,0),(3.2,0)]]
 

triqs/clef/clef.hpp

@@ -139,19 +139,21 @@ namespace triqs { namespace clef {
 
  // a little function to clean the reference_wrapper
  template<typename U> U  _cl(U && x) { return std::forward<U>(x);}
- template<typename U> U & _cl(std::reference_wrapper<U> x) { return x.get();}
+ template<typename U> auto _cl(std::reference_wrapper<U> x) DECL_AND_RETURN(x.get());
 
  /// Terminal 
  template<> struct operation<tags::terminal> { template<typename L> L operator()(L&& l) const { return std::forward<L>(l);} };
 
  /// Function call 
- template<> struct operation<tags::function> { 
+ template<> struct operation<tags::function> {
-  template<typename F, typename... Args> auto operator()(F const & f, Args const & ... args) const DECL_AND_RETURN(_cl(f)(_cl(args)...));
+  template <typename F, typename... Args>
+  auto operator()(F&& f, Args&&... args) const DECL_AND_RETURN(_cl(std::forward<F>(f))(_cl(std::forward<Args>(args))...));
  };
 
  /// [ ] Call
  template<> struct operation<tags::subscript> { 
-  template<typename F, typename Args> auto operator()(F const & f, Args const & args) const DECL_AND_RETURN(_cl(f)[_cl(args)]);
+  template <typename F, typename Args>
+  auto operator()(F&& f, Args&& args) const DECL_AND_RETURN(_cl(std::forward<F>(f))[_cl(std::forward<Args>(args))]);
  };
 
  // all binary operators....
@@ -228,7 +230,7 @@ namespace triqs { namespace clef {
  template<typename Tag, bool IsLazy> struct operation2;
  template<typename Tag> struct operation2<Tag, true> { 
   template<typename... Args> 
-   expr<Tag,typename remove_cv_ref<Args>::type ...> operator()(Args && ...  args) const  {
+   expr<Tag,typename expr_storage_t<Args>::type ...> operator()(Args && ...  args) const  {
     return {Tag(), std::forward<Args>(args)...};
    } 
  };

triqs/gfs/block.hpp

@@ -153,10 +153,37 @@ namespace gfs {
   // return { gf_mesh<block_index, Opt> {int(V.size())}, std::move(V), nothing{}, nothing{} } ;
  }
 
- template <typename GF, typename GF2> gf_view<block_index, GF> make_block_gf_view_from_vector(std::vector<GF2> V) {
+ template <typename GF> gf_view<block_index, typename GF::regular_type> make_block_gf_view_from_vector(std::vector<GF> V) {
   return {{int(V.size())}, std::move(V), nothing{}, nothing{}};
  }
 
+ // for cython proxy only. do not document.
+ template <typename GF, typename GF2> gf_view<block_index, GF> make_block_gf_view_from_vector_of_cython_proxy(std::vector<GF2> V) {
+  return {{int(V.size())}, std::move(V), nothing{}, nothing{}};
+ }
+
+ // ------------------------------- Extend reinterpret_scalar_valued_gf_as_matrix_valued for block gf   ------
+
+ template <typename Variable, typename Opt, typename Opt2, bool IsConst>
+ gf_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2, IsConst>
+ reinterpret_scalar_valued_gf_as_matrix_valued(gf_view<block_index, gf<Variable, scalar_valued, Opt>, Opt2, IsConst> bg) {
+  std::vector<gf_view<Variable, matrix_valued, Opt>> V;
+  for (auto &g : bg) V.push_back(reinterpret_scalar_valued_gf_as_matrix_valued(g));
+  return make_block_gf_view_from_vector(std::move(V));
+ }
+
+ template <typename Variable, typename Opt, typename Opt2>
+ gf_const_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2>
+ reinterpret_scalar_valued_gf_as_matrix_valued(gf<block_index, gf<Variable, scalar_valued, Opt>, Opt2> const &bg) {
+  return reinterpret_scalar_valued_gf_as_matrix_valued(bg());
+ }
+
+ template <typename Variable, typename Opt, typename Opt2>
+ gf_view<block_index, gf<Variable, matrix_valued, Opt>, Opt2>
+ reinterpret_scalar_valued_gf_as_matrix_valued(gf<block_index, gf<Variable, scalar_valued, Opt>, Opt2> &bg) {
+  return reinterpret_scalar_valued_gf_as_matrix_valued(bg());
+ }
+
  // -------------------------------   Free functions   --------------------------------------------------
 
  // a simple function to get the number of blocks
@@ -165,17 +192,19 @@ namespace gfs {
 
  // -------------------------------   an iterator over the blocks --------------------------------------------------
 
+ template<typename T> using __get_target = typename std::remove_reference<decltype(std::declval<T>()[0])>::type;
+
  // iterator
- template <typename Target, typename Opt, bool B, bool C>
+ template <typename G>
  class block_gf_iterator
-     : public boost::iterator_facade<block_gf_iterator<Target, Opt, B, C>, Target, boost::forward_traversal_tag, Target &> {
+     : public boost::iterator_facade<block_gf_iterator<G>, __get_target<G>, boost::forward_traversal_tag, __get_target<G> &> {
   friend class boost::iterator_core_access;
-  typedef gf_impl<block_index, Target, Opt, B, C> big_gf_t;
+  typedef typename std::remove_reference<G>::type big_gf_t;
   big_gf_t &big_gf;
   typedef typename big_gf_t::mesh_t::const_iterator mesh_iterator_t;
   mesh_iterator_t mesh_it;
 
-  Target &dereference() const { return big_gf[*mesh_it]; }
+  __get_target<G> &dereference() const { return big_gf[*mesh_it]; }
   bool equal(block_gf_iterator const &other) const { return ((mesh_it == other.mesh_it)); }
 
   public:
@@ -186,27 +215,27 @@ namespace gfs {
 
  //------------
  template <typename Target, typename Opt, bool B, bool C>
- block_gf_iterator<Target, Opt, B, C> begin(gf_impl<block_index, Target, Opt, B, C> &bgf) {
+ block_gf_iterator<gf_impl<block_index, Target, Opt, B, C>> begin(gf_impl<block_index, Target, Opt, B, C> &bgf) {
   return {bgf, false};
  }
 
  //------------
  template <typename Target, typename Opt, bool B, bool C>
- block_gf_iterator<Target, Opt, B, C> end(gf_impl<block_index, Target, Opt, B, C> &bgf) {
+ block_gf_iterator<gf_impl<block_index, Target, Opt, B, C>> end(gf_impl<block_index, Target, Opt, B, C> &bgf) {
   return {bgf, true};
  }
 
  //----- const iterator
- template <typename Target, typename Opt, bool B, bool C>
+ template <typename G>
- class block_gf_const_iterator : public boost::iterator_facade<block_gf_const_iterator<Target, Opt, B, C>, Target,
+ class block_gf_const_iterator
-                                                               boost::forward_traversal_tag, Target const &> {
+     : public boost::iterator_facade<block_gf_const_iterator<G>, __get_target<G>, boost::forward_traversal_tag, __get_target<G> const &> {
   friend class boost::iterator_core_access;
-  typedef gf_impl<block_index, Target, Opt, B, C> big_gf_t;
+  typedef typename std::remove_reference<G>::type big_gf_t;
   big_gf_t const &big_gf;
   typedef typename big_gf_t::mesh_t::const_iterator mesh_iterator_t;
   mesh_iterator_t mesh_it;
 
-  Target const &dereference() const { return big_gf[*mesh_it]; }
+  __get_target<G> const &dereference() const { return big_gf[*mesh_it]; }
   bool equal(block_gf_const_iterator const &other) const { return ((mesh_it == other.mesh_it)); }
 
   public:
@@ -216,22 +245,22 @@ namespace gfs {
  };
 
  template <typename Target, typename Opt, bool B, bool C>
- block_gf_const_iterator<Target, Opt, B, C> begin(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
+ block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> begin(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
   return {bgf, false};
  }
 
  template <typename Target, typename Opt, bool B, bool C>
- block_gf_const_iterator<Target, Opt, B, C> end(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
+ block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> end(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
   return {bgf, true};
  }
 
  template <typename Target, typename Opt, bool B, bool C>
- block_gf_const_iterator<Target, Opt, B, C> cbegin(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
+ block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> cbegin(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
   return {bgf, false};
  }
 
  template <typename Target, typename Opt, bool B, bool C>
- block_gf_const_iterator<Target, Opt, B, C> cend(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
+ block_gf_const_iterator<gf_impl<block_index, Target, Opt, B, C>> cend(gf_impl<block_index, Target, Opt, B, C> const &bgf) {
   return {bgf, true};
  }
 }

triqs/gfs/gf.hpp

@@ -546,7 +546,7 @@ namespace triqs { namespace gfs {
   }
 
   // a scalar_valued gf can be viewed as a 1x1 matrix
-  template <typename Variable, typename Opt, bool IsConst, typename... Args>
+  template <typename Variable, typename Opt, bool IsConst>
   gf_view<Variable, matrix_valued, Opt, IsConst>
   reinterpret_scalar_valued_gf_as_matrix_valued(gf_view<Variable, scalar_valued, Opt, IsConst> g) {
    typedef typename gf_view<Variable, matrix_valued, Opt, IsConst>::data_view_t a_t;
@@ -554,12 +554,12 @@ namespace triqs { namespace gfs {
    return {g.mesh(), a, g.singularity(), g.symmetry()};
   }
 
-  template <typename Variable, typename Opt, typename... Args>
+  template <typename Variable, typename Opt>
   gf_view<Variable, matrix_valued, Opt> reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Opt> &g) {
    return reinterpret_scalar_valued_gf_as_matrix_valued(g());
   }
 
-  template <typename Variable, typename Opt, typename... Args>
+  template <typename Variable, typename Opt>
   gf_const_view<Variable, matrix_valued, Opt>
   reinterpret_scalar_valued_gf_as_matrix_valued(gf<Variable, scalar_valued, Opt> const &g) {
    return reinterpret_scalar_valued_gf_as_matrix_valued(g());
@@ -574,6 +574,8 @@ namespace triqs { namespace gfs {
  namespace gfs_implementation { // implement some default traits
 
   // -------------------------  default factories ---------------------
+
+  // ----- tensor_valued
   template <int R, typename Var, typename Opt> struct factories<Var, tensor_valued<R>, Opt> {
    typedef gf<Var, tensor_valued<R>, Opt> gf_t;
    typedef tqa::mini_vector<size_t, R> target_shape_t;
@@ -590,6 +592,7 @@ namespace triqs { namespace gfs {
    }
   };
 
+  // ----- matrix_valued
   template <typename Var, typename Opt> struct factories<Var, matrix_valued, Opt> {
    typedef gf<Var, matrix_valued, Opt> gf_t;
    typedef tqa::mini_vector<size_t, 2> target_shape_t;
@@ -606,6 +609,7 @@ namespace triqs { namespace gfs {
    }
   };
 
+  // ----- scalar_valued
   template <typename Var, typename Opt> struct factories<Var, scalar_valued, Opt> {
    typedef gf<Var, scalar_valued, Opt> gf_t;
    struct target_shape_t {