From c06e183ae52eaf36f928fdce7c41f0528bfc164d Mon Sep 17 00:00:00 2001
From: Olivier Parcollet <olivier.parcollet@cea.fr>
Date: Tue, 10 Sep 2013 11:18:04 +0200
Subject: [PATCH] arrays : revised version of map & fold

- also cleaned a bit reinterpretview example
---
 test/triqs/arrays/fold_expression.cpp     |   2 +-
 test/triqs/arrays/group_indices.cpp       |   9 +-
 test/triqs/arrays/mapped_functions.cpp    |  18 +-
 test/triqs/arrays/mapped_functions.output | 197 ++++++++++++----------
 test/triqs/arrays/reinterpret_view.cpp    |  40 ++++-
 triqs/arrays/algorithms.hpp               |   8 +-
 triqs/arrays/asserts.hpp                  |   2 +-
 triqs/arrays/blas_lapack/getri.hpp        |   4 +-
 triqs/arrays/functional/fold.hpp          |  35 ++--
 triqs/arrays/functional/map.hpp           | 149 +++++++---------
 triqs/arrays/mapped_functions.hpp         |  86 +++-------
 triqs/utility/concept_tools.hpp           |   9 +-
 12 files changed, 274 insertions(+), 285 deletions(-)

diff --git a/test/triqs/arrays/fold_expression.cpp b/test/triqs/arrays/fold_expression.cpp
index 7480707a..5d5b1be1 100644
--- a/test/triqs/arrays/fold_expression.cpp
+++ b/test/triqs/arrays/fold_expression.cpp
@@ -28,7 +28,7 @@
 
 using namespace triqs::arrays;
 
-auto Abs = map( boost::function< double (double)> ( static_cast< double (*)(double)> (std::abs)) );
+auto Abs = map( std::function< double (double)> ( static_cast< double (*)(double)> (std::abs)) );
 
 int main(int argc, char **argv) {
 
diff --git a/test/triqs/arrays/group_indices.cpp b/test/triqs/arrays/group_indices.cpp
index d48b5dbc..a88cdcc2 100644
--- a/test/triqs/arrays/group_indices.cpp
+++ b/test/triqs/arrays/group_indices.cpp
@@ -5,7 +5,6 @@
 
 namespace tqa=triqs::arrays;
 namespace tql=triqs::clef;
-namespace mpl=boost::mpl;
 
 using tqa::m_index;
 
@@ -57,7 +56,9 @@ int main () {
  test<TRAVERSAL_ORDER_FORTRAN>();
  test<TRAVERSAL_ORDER_C>();
 
- // an alternative syntax ? Why to use the npl here ??
- //auto V = tqa::group_indices( A(i_,j_,k_,l_), m_index(i_,j_), m_index(k_,l_));
-
+ // template dans le namespace clef pour ADL et enable sur le concept de A ...
+ // an alternative syntax ? 
+ //auto V = group_indices( A(i_,j_,k_,l_), var(i_,j_), var(k_,l_));
+ // match A(...), deduce the int of hte ph, and make a finder for a ph in it
+ // finder<Expr>::find<i>::value --> m_index<...>
 }
diff --git a/test/triqs/arrays/mapped_functions.cpp b/test/triqs/arrays/mapped_functions.cpp
index 124f29dd..faf8ab72 100644
--- a/test/triqs/arrays/mapped_functions.cpp
+++ b/test/triqs/arrays/mapped_functions.cpp
@@ -22,13 +22,14 @@
 
 #include "./src/mapped_functions.hpp"
 #include "./src/matrix.hpp"
+#include "./src/array.hpp"
 #include <iostream>
 
-namespace tqa = triqs::arrays;
+using namespace triqs::arrays;
 
 template<typename T> void test( T val=1 ) { 
- tqa::matrix<T> A(3,3, FORTRAN_LAYOUT),B(3,3, FORTRAN_LAYOUT);
-
+ matrix<T> A(3,3, FORTRAN_LAYOUT),B(3,3, FORTRAN_LAYOUT);
+ typedef array_view<T,2> V;
  for (int i =0; i<3; ++i)
   for (int j=0; j<3; ++j)
   { A(i,j) = (i+2*j+1); B(i,j) = (i-j);}
@@ -36,8 +37,9 @@ template<typename T> void test( T val=1 ) {
  A *=val; B*=val;
  T s = 10;
  TEST(A);
- TEST(make_matrix(pow(A,2)));
- TEST(make_matrix(cosh(A)));
+
+ TEST(make_matrix(pow(V(A),2)));
+ TEST(make_matrix(cosh(V(A))));
  TEST(B);
  TEST(abs(B));
  TEST(make_matrix(abs(B)));
@@ -48,6 +50,12 @@ template<typename T> void test( T val=1 ) {
  TEST(make_matrix(real(B)));
  TEST(make_matrix(imag(B)));
 
+ auto aa = array<T,2>{ { 1,2}, {3,4}};
+ TEST(make_matrix(exp(aa)));
+
+ // does not compile, since exp is only element wise at the moment
+ // to do : implement it for matrix...
+ //TEST(make_matrix(exp( matrix<double>{{ 1,2}, {3,4}} )));
 }
 
 int main(int argc, char **argv) {
diff --git a/test/triqs/arrays/mapped_functions.output b/test/triqs/arrays/mapped_functions.output
index 3fa85a11..b6584857 100644
--- a/test/triqs/arrays/mapped_functions.output
+++ b/test/triqs/arrays/mapped_functions.output
@@ -2,93 +2,11 @@
 [[1,3,5]
  [2,4,6]
  [3,5,7]]
-(make_matrix(pow(A,2))) ---> 
+(make_matrix(pow(V(A),2))) ---> 
 [[1,9,25]
  [4,16,36]
  [9,25,49]]
-(make_matrix(cosh(A))) ---> 
-[[1,10,74]
- [3,27,201]
- [10,74,548]]
-(B) ---> 
-[[0,-1,-2]
- [1,0,-1]
- [2,1,0]]
-(abs(B)) ---> lazy matrix resulting of a mapping
-(make_matrix(abs(B))) ---> 
-[[0,1,2]
- [1,0,1]
- [2,1,0]]
-(make_matrix(abs(B+B))) ---> 
-[[0,2,4]
- [2,0,2]
- [4,2,0]]
-(make_matrix(A+ s*B)) ---> 
-[[1,-7,-15]
- [12,4,-4]
- [23,15,7]]
-(make_matrix(abs(A+s*B))) ---> 
-[[1,7,15]
- [12,4,4]
- [23,15,7]]
-(make_matrix(real(B))) ---> 
-[[0,-1,-2]
- [1,0,-1]
- [2,1,0]]
-(make_matrix(imag(B))) ---> 
-[[0,0,0]
- [0,0,0]
- [0,0,0]]
-(A) ---> 
-[[1,3,5]
- [2,4,6]
- [3,5,7]]
-(make_matrix(pow(A,2))) ---> 
-[[1,9,25]
- [4,16,36]
- [9,25,49]]
-(make_matrix(cosh(A))) ---> 
-[[1,10,74]
- [3,27,201]
- [10,74,548]]
-(B) ---> 
-[[0,-1,-2]
- [1,0,-1]
- [2,1,0]]
-(abs(B)) ---> lazy matrix resulting of a mapping
-(make_matrix(abs(B))) ---> 
-[[0,1,2]
- [1,0,1]
- [2,1,0]]
-(make_matrix(abs(B+B))) ---> 
-[[0,2,4]
- [2,0,2]
- [4,2,0]]
-(make_matrix(A+ s*B)) ---> 
-[[1,-7,-15]
- [12,4,-4]
- [23,15,7]]
-(make_matrix(abs(A+s*B))) ---> 
-[[1,7,15]
- [12,4,4]
- [23,15,7]]
-(make_matrix(real(B))) ---> 
-[[0,-1,-2]
- [1,0,-1]
- [2,1,0]]
-(make_matrix(imag(B))) ---> 
-[[0,0,0]
- [0,0,0]
- [0,0,0]]
-(A) ---> 
-[[1,3,5]
- [2,4,6]
- [3,5,7]]
-(make_matrix(pow(A,2))) ---> 
-[[1,9,25]
- [4,16,36]
- [9,25,49]]
-(make_matrix(cosh(A))) ---> 
+(make_matrix(cosh(V(A)))) ---> 
 [[1.54308,10.0677,74.2099]
  [3.7622,27.3082,201.716]
  [10.0677,74.2099,548.317]]
@@ -96,7 +14,7 @@
 [[0,-1,-2]
  [1,0,-1]
  [2,1,0]]
-(abs(B)) ---> lazy matrix resulting of a mapping
+(abs(B)) ---> mapping result
 (make_matrix(abs(B))) ---> 
 [[0,1,2]
  [1,0,1]
@@ -121,15 +39,106 @@
 [[0,0,0]
  [0,0,0]
  [0,0,0]]
+(make_matrix(exp(aa))) ---> 
+[[2.71828,7.38906]
+ [20.0855,54.5982]]
+(A) ---> 
+[[1,3,5]
+ [2,4,6]
+ [3,5,7]]
+(make_matrix(pow(V(A),2))) ---> 
+[[1,9,25]
+ [4,16,36]
+ [9,25,49]]
+(make_matrix(cosh(V(A)))) ---> 
+[[1.54308,10.0677,74.2099]
+ [3.7622,27.3082,201.716]
+ [10.0677,74.2099,548.317]]
+(B) ---> 
+[[0,-1,-2]
+ [1,0,-1]
+ [2,1,0]]
+(abs(B)) ---> mapping result
+(make_matrix(abs(B))) ---> 
+[[0,1,2]
+ [1,0,1]
+ [2,1,0]]
+(make_matrix(abs(B+B))) ---> 
+[[0,2,4]
+ [2,0,2]
+ [4,2,0]]
+(make_matrix(A+ s*B)) ---> 
+[[1,-7,-15]
+ [12,4,-4]
+ [23,15,7]]
+(make_matrix(abs(A+s*B))) ---> 
+[[1,7,15]
+ [12,4,4]
+ [23,15,7]]
+(make_matrix(real(B))) ---> 
+[[0,-1,-2]
+ [1,0,-1]
+ [2,1,0]]
+(make_matrix(imag(B))) ---> 
+[[0,0,0]
+ [0,0,0]
+ [0,0,0]]
+(make_matrix(exp(aa))) ---> 
+[[2.71828,7.38906]
+ [20.0855,54.5982]]
+(A) ---> 
+[[1,3,5]
+ [2,4,6]
+ [3,5,7]]
+(make_matrix(pow(V(A),2))) ---> 
+[[1,9,25]
+ [4,16,36]
+ [9,25,49]]
+(make_matrix(cosh(V(A)))) ---> 
+[[1.54308,10.0677,74.2099]
+ [3.7622,27.3082,201.716]
+ [10.0677,74.2099,548.317]]
+(B) ---> 
+[[0,-1,-2]
+ [1,0,-1]
+ [2,1,0]]
+(abs(B)) ---> mapping result
+(make_matrix(abs(B))) ---> 
+[[0,1,2]
+ [1,0,1]
+ [2,1,0]]
+(make_matrix(abs(B+B))) ---> 
+[[0,2,4]
+ [2,0,2]
+ [4,2,0]]
+(make_matrix(A+ s*B)) ---> 
+[[1,-7,-15]
+ [12,4,-4]
+ [23,15,7]]
+(make_matrix(abs(A+s*B))) ---> 
+[[1,7,15]
+ [12,4,4]
+ [23,15,7]]
+(make_matrix(real(B))) ---> 
+[[0,-1,-2]
+ [1,0,-1]
+ [2,1,0]]
+(make_matrix(imag(B))) ---> 
+[[0,0,0]
+ [0,0,0]
+ [0,0,0]]
+(make_matrix(exp(aa))) ---> 
+[[2.71828,7.38906]
+ [20.0855,54.5982]]
 (A) ---> 
 [[(1,0),(3,0),(5,0)]
  [(2,0),(4,0),(6,0)]
  [(3,0),(5,0),(7,0)]]
-(make_matrix(pow(A,2))) ---> 
+(make_matrix(pow(V(A),2))) ---> 
 [[(1,0),(9,0),(25,0)]
  [(4,0),(16,0),(36,0)]
  [(9,0),(25,0),(49,0)]]
-(make_matrix(cosh(A))) ---> 
+(make_matrix(cosh(V(A)))) ---> 
 [[(1.54308,0),(10.0677,0),(74.2099,0)]
  [(3.7622,0),(27.3082,0),(201.716,0)]
  [(10.0677,0),(74.2099,0),(548.317,0)]]
@@ -137,7 +146,7 @@
 [[(0,0),(-1,0),(-2,0)]
  [(1,0),(0,0),(-1,0)]
  [(2,0),(1,0),(0,0)]]
-(abs(B)) ---> lazy matrix resulting of a mapping
+(abs(B)) ---> mapping result
 (make_matrix(abs(B))) ---> 
 [[0,1,2]
  [1,0,1]
@@ -162,15 +171,18 @@
 [[0,0,0]
  [0,0,0]
  [0,0,0]]
+(make_matrix(exp(aa))) ---> 
+[[(2.71828,0),(7.38906,0)]
+ [(20.0855,0),(54.5982,0)]]
 (A) ---> 
 [[(-2,3),(-6,9),(-10,15)]
  [(-4,6),(-8,12),(-12,18)]
  [(-6,9),(-10,15),(-14,21)]]
-(make_matrix(pow(A,2))) ---> 
+(make_matrix(pow(V(A),2))) ---> 
 [[(-5,-12),(-45,-108),(-125,-300)]
  [(-20,-48),(-80,-192),(-180,-432)]
  [(-45,-108),(-125,-300),(-245,-588)]]
-(make_matrix(cosh(A))) ---> 
+(make_matrix(cosh(V(A)))) ---> 
 [[(-3.72455,-0.511823),(-183.789,-83.1297),(-8366.62,-7161.77)]
  [(26.2206,7.62523),(1257.75,799.751),(53734.9,61113.4)]
  [(-183.789,-83.1297),(-8366.62,-7161.77),(-329351,-503083)]]
@@ -178,7 +190,7 @@
 [[(0,0),(2,-3),(4,-6)]
  [(-2,3),(0,0),(2,-3)]
  [(-4,6),(-2,3),(0,0)]]
-(abs(B)) ---> lazy matrix resulting of a mapping
+(abs(B)) ---> mapping result
 (make_matrix(abs(B))) ---> 
 [[0,3.60555,7.2111]
  [3.60555,0,3.60555]
@@ -203,3 +215,6 @@
 [[0,-3,-6]
  [3,0,-3]
  [6,3,0]]
+(make_matrix(exp(aa))) ---> 
+[[(2.71828,0),(7.38906,0)]
+ [(20.0855,0),(54.5982,0)]]
diff --git a/test/triqs/arrays/reinterpret_view.cpp b/test/triqs/arrays/reinterpret_view.cpp
index 7bfc6abd..545ea067 100644
--- a/test/triqs/arrays/reinterpret_view.cpp
+++ b/test/triqs/arrays/reinterpret_view.cpp
@@ -2,7 +2,7 @@
  *
  * TRIQS: a Toolbox for Research in Interacting Quantum Systems
  *
- * Copyright (C) 2011 by O. Parcollet
+ * 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
@@ -24,18 +24,46 @@
 
 namespace triqs { namespace arrays { 
 
- template <typename V, int R, ull_t Opt, ull_t To, typename ... I>  
+/* template <typename V, int R, ull_t Opt, ull_t To, typename ... I>  
   array_view<V, sizeof...(I), Opt> reinterpret_array_view (array<V,R,Opt,To> const & a, I ... index) { 
-   static int constexpr rank = sizeof...(I);
-   typedef array_view<V, rank, Opt, indexmaps::mem_layout::c_order(rank)> return_t;
-   return return_t (typename return_t::indexmap_type (mini_vector<size_t,rank>(index...)) , a.storage());
+   //static int constexpr rank = sizeof...(I);
+   typedef array_view<V, sizeof...(I), Opt, indexmaps::mem_layout::c_order(rank)> return_t;
+   //typedef array_view<V, rank, Opt, indexmaps::mem_layout::c_order(rank)> return_t;
+   return return_t { make_shape(index...), a.storage() };
+   //return return_t (typename return_t::indexmap_type (mini_vector<size_t,rank>(index...)) , a.storage());
   }
+*/
+
+ template <typename V, int R, ull_t Opt, ull_t To, typename ... I>  
+  array_view<V, sizeof...(I), Opt> reinterpret (array<V,R,Opt,To> const & a, I ... index) { 
+   return { {make_shape(index...)}, a.storage() };
+  }
+
+ // wrong for views
+ template <typename V, int R, ull_t Opt, ull_t To, bool B, typename ... I>  
+  array_view<V, sizeof...(I), Opt,indexmaps::mem_layout::c_order(sizeof...(I))> reinterpret_array_view (array_view<V,R,Opt,To,B> const & a, I ... index) { 
+   if (!has_contiguous_data(a)) TRIQS_RUNTIME_ERROR << "reinterpretation failure : data of the view are not contiguous";
+   return { {make_shape(index...)}, a.storage() };
+  }
+/*
+ template <typename V, int R, ull_t Opt, ull_t To, typename ... I>  
+  array_view<V, sizeof...(I), Opt,indexmaps::mem_layout::c_order(sizeof...(I))> reinterpret_array_view_c (matrix<V,Opt,To> const & a, I ... index) { 
+   return { {make_shape(index...)}, a.storage() };
+  }
+
+template <typename V, int R, ull_t Opt, ull_t To, typename ... I>  
+  array_view<V, sizeof...(I), Opt,indexmaps::mem_layout::c_order(sizeof...(I))> reinterpret_array_view_c (vector<V,Opt> const & a, I ... index) { 
+   return { {make_shape(index...)}, a.storage() };
+  }
+*/
 
 }}
 
+using namespace triqs::arrays;
 int main(int argc, char **argv) {
  triqs::arrays::array<long,1> A = {1,2,3,4,5,6};
- std::cout << reinterpret_array_view(A, 2,3);
+ std::cout << reinterpret(A, 2,3)<< std::endl;
+ //std::cout << reinterpret_array_view(A(range(1,5)), 2,3)<< std::endl;
  return 0;
 }
 
diff --git a/triqs/arrays/algorithms.hpp b/triqs/arrays/algorithms.hpp
index 8779c6ca..18727c24 100644
--- a/triqs/arrays/algorithms.hpp
+++ b/triqs/arrays/algorithms.hpp
@@ -39,14 +39,14 @@ namespace triqs { namespace arrays {
 
  template<class A>
   typename A::value_type max_element(A const &a) { 
-   //typedef typename A::value_type const & T; return fold ( static_cast<T(*)(T,T)>(std::max<T>))  ( a, get_first_element(a)); 
-   return fold (std::max<typename A::value_type const & >)  (a, get_first_element(a)); 
+   typedef typename A::value_type T;
+   return fold ([](T const & a, T const & b) { return std::max(a,b);} )  (a, get_first_element(a)); 
   }
 
  template<class A>
   typename A::value_type min_element(A const &a) { 
-   // typedef typename A::value_type const & T; return fold ( static_cast<T(*)(T,T)>(std::min<T>))  ( a, get_first_element(a)); 
-   return fold (std::min<typename A::value_type const & >)  (a, get_first_element(a)); 
+   typedef typename A::value_type T;
+   return fold ([](T const & a, T const & b) { return std::min(a,b);} )  (a, get_first_element(a)); 
   }
 
  template <class A>
diff --git a/triqs/arrays/asserts.hpp b/triqs/arrays/asserts.hpp
index 32766494..0f732cc4 100644
--- a/triqs/arrays/asserts.hpp
+++ b/triqs/arrays/asserts.hpp
@@ -30,7 +30,7 @@ namespace triqs { namespace arrays {
  template<class ArrayType1, class ArrayType2 >
   void assert_all_close( ArrayType1 const & A, ArrayType2 const & B, double precision, bool relative = false) {
    typedef typename ArrayType1::value_type F;
-   auto Abs =  map( boost::function<double(F)> (assert_abs<F>) );
+   auto Abs =  map( std::function<double(F)> (assert_abs<F>) );
    auto r =  max_element (Abs(A-B));
    auto r2 =  max_element (Abs(A) + Abs(B));
    if ( r > (relative ? precision * r2 : precision) )
diff --git a/triqs/arrays/blas_lapack/getri.hpp b/triqs/arrays/blas_lapack/getri.hpp
index c76690ce..45a6a0df 100644
--- a/triqs/arrays/blas_lapack/getri.hpp
+++ b/triqs/arrays/blas_lapack/getri.hpp
@@ -42,8 +42,8 @@ namespace triqs { namespace arrays { namespace lapack {
   }
  }
 
- inline size_t r_round(double x) { return round(x);}
- inline size_t r_round(std::complex<double> x) { return round(real(x));}
+ inline size_t r_round(double x) { return std::round(x);}
+ inline size_t r_round(std::complex<double> x) { return std::round(std::real(x));}
  /**
   * Calls getri on a matrix or view
   * Takes care of making temporary copies if necessary
diff --git a/triqs/arrays/functional/fold.hpp b/triqs/arrays/functional/fold.hpp
index be11eaa2..126a2ae7 100644
--- a/triqs/arrays/functional/fold.hpp
+++ b/triqs/arrays/functional/fold.hpp
@@ -28,37 +28,24 @@
 namespace triqs { namespace arrays {
 
  template<class F> 
-  class fold_worker  { 
+  struct fold_worker  { 
    F f;
-   typedef typename boost::remove_const<typename boost::remove_reference<typename F::result_type >::type>::type result_type;
    
-   template<class A>
-    struct fold_func_adaptor { 
-     A const & b;
-     F f; result_type r;
-     fold_func_adaptor(F f_, A const & a_,typename A::value_type init ):b(a_), f(f_) { r= init;}
-     template<typename ... Args> void operator()(Args const & ... args) { r = f(r,b(args...)); }
+   template<class A, class R> struct fold_func_adaptor { 
+     F const & f; A const & a; R & r;
+     template<typename ... Args> void operator()(Args const & ... args) { r = f(r,a(args...)); }
     };
 
-   public:
-
-   fold_worker ( F const & f_):f(f_) {} 
-
-   template<class A>   
-    result_type operator() (A const & a, typename A::value_type init = typename A::value_type() )  const { 
-     fold_func_adaptor<A> func(f,a,init);
-     //fold_func_adaptor<typename A::value_type> func(f,init);
-     foreach(a,std::ref(func));
-     //foreach(a,boost::ref(func));
-     return func.r;
+   template<class A, class R>   
+    R operator() (A const & a, R init)  const { 
+     foreach(a, fold_func_adaptor<A,R> {f,a,init});
+     return init;
     }
+
+   template<class A> typename A::value_type operator() (A const & a)  const { return (*this)(a, typename A::value_type{});}
   };
 
- template<class F> fold_worker<F> fold (F const & f) { return fold_worker<F>(f);}
-
- template<typename R, typename A1> 
-  fold_worker< boost::function<R(A1,A1)> >
-  fold (R (*f)(A1,A1)) { return fold( boost::function<R(A1,A1)>(f)); } 
+ template<class F> fold_worker<F> fold (F f) { return {std::move(f)};}
 
 }}//namespace triqs::arrays
 
diff --git a/triqs/arrays/functional/map.hpp b/triqs/arrays/functional/map.hpp
index 23e3becb..73e2d93f 100644
--- a/triqs/arrays/functional/map.hpp
+++ b/triqs/arrays/functional/map.hpp
@@ -22,6 +22,8 @@
 #define TRIQS_ARRAYS_EXPRESSION_MAP_H
 #include "../impl/common.hpp"
 #include <functional>
+//#include "../../utility/function_arg_ret_type.hpp"
+
 namespace triqs { namespace arrays { 
  
  template<class F, int arity=F::arity> class map_impl;
@@ -30,105 +32,80 @@ namespace triqs { namespace arrays {
   * Given a function f : arg_type -> result_type, map(f) is the function promoted to arrays
   * map(f) : array<arg_type, N, Opt> --> array<result_type, N, Opt> 
   */
- template<class F> map_impl<F,1> map (F const & f) { return map_impl<F,1>(f,true); }
- template<class F> map_impl<F,2> map2 (F const & f) { return map_impl<F,2>(f,true); }
+ //template<class F> map_impl<F,utility::function_arg_ret_type<F>::arity> map  (F f) { return {std::move(f),true}; }
 
+ template<class F> map_impl<F,1> map  (F f) { return {std::move(f),true}; }
+ template<class F> map_impl<F,2> map2 (F f) { return {std::move(f),true}; }
+ 
  // ----------- implementation  -------------------------------------
 
+ template<typename F, int arity, bool is_vec, typename ... A> struct map_impl_result;
+
+ template<typename F, bool is_vec, typename A> struct map_impl_result<F,1,is_vec,A> { 
+  typedef typename std::result_of<F(typename std::remove_reference<A>::type::value_type)>::type value_type;
+  typedef typename std::remove_reference<A>::type::domain_type domain_type;
+  F f; 
+  typename std::add_const<A>::type a; // A is a T or a T& : add const to them.
+  domain_type domain() const { return a.domain(); } 
+  template<typename ... Args> value_type operator() (Args && ... args) const { return f(a(std::forward<Args>(args)...)); }
+  friend std::ostream & operator<<(std::ostream & out, map_impl_result const & x){ return out<<"mapping result";}
+  // rest is only for vector 
+  template<bool vec = is_vec> 
+   TYPE_ENABLE_IFC(size_t,vec) size() const { return a.size();}
+  template<typename Args, bool vec=is_vec>  
+   TYPE_ENABLE_IFC(value_type,vec) operator[] (Args && args) const { return f(a[std::forward<Args>(args)]);}
+ };
+
+ // possible to generalize to N order using tuple techniques ...
+ template<typename F, bool is_vec, typename A, typename B> struct map_impl_result<F,2,is_vec,A,B> { 
+  typedef typename std::result_of<F(typename remove_cv_ref<A>::type::value_type, typename remove_cv_ref<B>::type::value_type)>::type value_type;
+  typedef typename remove_cv_ref<A>::type::domain_type domain_type;
+  F f; 
+  typename std::add_const<A>::type a;
+  typename std::add_const<B>::type b;
+  domain_type domain() const { return a.domain(); } 
+  template<typename ... Args> value_type operator() (Args && ... args) const { return f(a(std::forward<Args>(args)...),b(std::forward<Args>(args)...)); }
+  friend std::ostream & operator<<(std::ostream & out, map_impl_result const & x){ return out<<"mapping result";}
+  // rest is only for vector 
+  template<bool vec = is_vec> 
+   TYPE_ENABLE_IFC(size_t,vec) size() const { return a.size();}
+  template<typename Args, bool vec=is_vec>  
+   TYPE_ENABLE_IFC(value_type,vec) operator[] (Args && args) const { return f(a[std::forward<Args>(args)],b[std::forward<Args>(args)]);}
+ };
+
+ template<typename ... T> struct _and; 
+ template<typename T0, typename ... T> struct _and<T0, T...> : std::integral_constant<bool, T0::value && _and<T...>::value>{}; 
+ template<typename T> struct _and<T> : T{};
+
+ template<typename F, int arity, bool b, typename ... A> struct ImmutableCuboidArray<map_impl_result<F,arity,b,A...>> : std::true_type{};
+ 
+ template<typename F, int arity, bool b, typename ... A> struct ImmutableArray <map_impl_result<F,arity,b,A...>> : _and<typename ImmutableArray <typename std::remove_reference<A>::type>::type...>{};
+ template<typename F, int arity, bool b, typename ... A> struct ImmutableMatrix<map_impl_result<F,arity,b,A...>> : _and<typename ImmutableMatrix<typename std::remove_reference<A>::type>::type...>{};
+ template<typename F, int arity, bool b, typename ... A> struct ImmutableVector<map_impl_result<F,arity,b,A...>> : _and<typename ImmutableVector<typename std::remove_reference<A>::type>::type...>{};
+ 
+ //template<typename F, int arity, bool b, typename A> struct ImmutableArray <map_impl_result<F,arity,b,A>> : ImmutableArray <A>{};
+ //template<typename F, int arity, bool b, typename A> struct ImmutableMatrix<map_impl_result<F,arity,b,A>> : ImmutableMatrix<A>{};
+ //template<typename F, int arity, bool b, typename A> struct ImmutableVector<map_impl_result<F,arity,b,A>> : ImmutableVector<A>{};
+
  // NB The bool is to make constructor not ambiguous
  // clang on os X with lib++ has a pb otherwise (not clear what the pb is)
- template<class F> class map_impl<F,1>  { 
+ template<class F, int arity> class map_impl  { 
   F f;
   public :   
-//  map_impl(F const & f_):f(f_) {}
-  map_impl(F const & f_, bool):f(f_) {}
-  map_impl(F && f_, bool):f(f_) {}
+  map_impl(F f_, bool):f(std::move(f_)) {}
   map_impl(map_impl const &) = default;
-  //map_impl(map_impl &&) = default;
+  map_impl(map_impl &&)      = default;
+  map_impl & operator = (map_impl const &) = default;
+  map_impl & operator = (map_impl &&)      = default;
 
-  template<typename A, typename Enable = void> class m_result;
-
-  template<typename A> class m_result<A,typename boost::enable_if<ImmutableArray<A> >::type> : TRIQS_CONCEPT_TAG_NAME(ImmutableArray) { 
-    public:
-     typedef typename std::result_of<F(typename A::value_type)>::type value_type;
-     typedef typename A::domain_type domain_type;
-     A const & a; F f;
-     m_result(F const & f_, A const & a_):a(a_),f(f_) {}
-     domain_type domain() const { return a.domain(); } 
-     template<typename ... Args> value_type operator() (Args const & ... args) const { return f(a(args...)); }
-     friend std::ostream & operator<<(std::ostream & out, m_result const & x){ return out<<"lazy matrix resulting of a mapping";}
-   };
-  
-  template<typename A> class m_result<A,typename boost::enable_if<ImmutableMatrix<A> >::type> : TRIQS_CONCEPT_TAG_NAME(ImmutableMatrix) { 
-    public:
-     typedef typename std::result_of<F(typename A::value_type)>::type value_type;
-     typedef typename A::domain_type domain_type;
-     A const & a; F f;
-     m_result(F const & f_, A const & a_):a(a_),f(f_) {}
-     domain_type domain() const { return a.domain(); } 
-     template<typename ... Args> value_type operator() (Args const & ... args) const { return f(a(args...)); }
-     friend std::ostream & operator<<(std::ostream & out, m_result const & x){ return out<<"lazy matrix resulting of a mapping";}
-   };
-
-  template<typename A> class m_result<A,typename boost::enable_if<ImmutableVector<A> >::type> : TRIQS_CONCEPT_TAG_NAME(ImmutableVector) { 
-    public:
-     typedef typename std::result_of<F(typename A::value_type)>::type value_type;
-     typedef typename A::domain_type domain_type;
-     A const & a; F f;
-     m_result(F const & f_, A const & a_):a(a_),f(f_) {}
-     domain_type domain() const { return a.domain(); } 
-     size_t size() const { return a.size();}
-     template<typename ... Args> value_type operator() (Args const & ... args) const { return f(a(args...)); }
-     value_type operator() ( size_t i) const { return f(a(i)); }
-     friend std::ostream & operator<<(std::ostream & out, m_result const & x){ return out<<"lazy matrix resulting of a mapping";}
-   };
-  
-
-  template<typename Sig> struct result;
-  template<typename This, typename A> struct result<This(A)> { typedef m_result<A> type;};
-
-  template< class A > m_result<A> operator()(A const & a) const { 
-    //static_assert( (ImmutableCuboidArray<A>::value), "map : A does not model ImmutableCuboidArray");
-    return m_result<A>(f,a);
-   } 
+  template<typename ... A> map_impl_result<F,arity,_and<typename ImmutableVector<A>::type...>::value,A...> 
+   operator()(A&&... a) const { return {f,std::forward<A>(a)...}; } 
 
   friend std::ostream & operator<<(std::ostream & out, map_impl const & x){ return out<<"map("<<"F"<<")";}
  };
 
- // ----------- // TO DO : make this with preprocessor ....
- template<class F> class map_impl<F,2>  { 
-  F f;
-  public : 
-  map_impl(F const & f_, bool):f(f_) {}
-  map_impl(F && f_, bool):f(f_) {}
-  map_impl(map_impl const &) = default;
-  //map_impl(map_impl &&) = default;
-
-  template<class A, class B> class m_result : TRIQS_CONCEPT_TAG_NAME(ImmutableCuboidArray) { 
-    static_assert( (std::is_same<typename  A::domain_type, typename  B::domain_type>::value), "type mismatch");
-   public:
-    typedef typename std::result_of<F(typename A::value_type,typename B::value_type)>::type value_type;
-    typedef typename A::domain_type domain_type;
-    A const & a; B const & b; F f;
-    m_result(F const & f_, A const & a_, B const & b_):a(a_),b(b_),f(f_) {
-     if (a.domain() != b.domain()) TRIQS_RUNTIME_ERROR<<"map2 : domain mismatch";
-    }
-    domain_type domain() const { return a.domain(); } 
-    template<typename ... Args> value_type operator() (Args const & ... args) const { return f(a(args...),b(args...)); }
-    friend std::ostream & operator<<(std::ostream & out, m_result const & x){ return out<<"lazy matrix resulting of a mapping";}
-  };
-
-  template<typename Sig> struct result;
-  template<typename This, typename A,typename B> struct result<This(A,B)> { typedef m_result<A,B> type;};
- 
-  template< class A, class B > m_result<A,B> operator()(A const & a, B const & b) { 
-    static_assert( (ImmutableCuboidArray<A>::value), "map1 : A does not model ImmutableCuboidArray");
-    static_assert( (ImmutableCuboidArray<B>::value), "map1 : B does not model ImmutableCuboidArray");
-    return m_result<A,B>(f,a,b);
-   } 
-
-  friend std::ostream & operator<<(std::ostream & out, map_impl const & x){ return out<<"map("<<"F"<<")";}
- };
 }}//namespace triqs::arrays
 #endif
 
+
+
diff --git a/triqs/arrays/mapped_functions.hpp b/triqs/arrays/mapped_functions.hpp
index 514b82c2..2b68ba47 100644
--- a/triqs/arrays/mapped_functions.hpp
+++ b/triqs/arrays/mapped_functions.hpp
@@ -22,84 +22,52 @@
 #define TRIQS_ARRAYS_MAPPED_FNT_H
 #include "./functional/map.hpp"
 #include <boost/preprocessor/seq/for_each.hpp>
-#include <boost/utility/result_of.hpp>
-
 namespace triqs { namespace arrays {
 
- using std::abs;
- struct abs_wrap { 
-   template<typename Sig> struct result;
-   template<typename This, typename A> struct result<This(A)> { typedef A type;};
-   template<typename This, typename A> struct result<This(std::complex<A>)> { typedef A type;};
-   static const int arity =1;
-   template<typename A> typename result<abs_wrap(A)>::type operator()( A const & a) const { return abs(a);}
- };
-
- template <typename A>
-  typename std::result_of<map_impl<abs_wrap>(A)>::type abs(A const & a) { return map(abs_wrap())(a); }
+ //C++14 will simply be ... 
+ //template <typename A> decltype(auto) abs(A && a) { return map( [](auto const &x) { using std::abs; return abs(a);}, std::forward<A>(a));}
 
  using std::pow;
  template<typename T> TYPE_ENABLE_IF(T, std::is_integral<T>) pow(T x, int n) { return (n==0 ? 1 : pow(x,n-1)*x);}
  struct pow_wrap { 
-   template<typename Sig> struct result;
-   template<typename This, typename A> struct result<This(A)> { typedef A type;};
-   static const int arity =1;
-   int n; pow_wrap(int n_): n(n_){}
-   template<typename A> typename result<abs_wrap(A)>::type operator()( A const & a) const { return pow(a,n);}
+   int n; 
+   template<typename A> A operator()( A const & a) const { return pow(a,n);}
  };
 
  template <typename A>
-  typename boost::result_of<map_impl<pow_wrap>(A)>::type pow(A const & a, int n) { return map(pow_wrap(n))(a); }
+  typename boost::lazy_enable_if_c<!ImmutableMatrix<A>::value,std::result_of<map_impl<pow_wrap,1>(A)>>::type
+  pow(A && a, int n) { return map(pow_wrap{n})(std::forward<A>(a)); }
 
 #define MAP_IT(FNT) \
  using std::FNT;\
- struct FNT##_wrap {\
-   template<typename Sig> struct result;\
-   template<typename This, typename A> struct result<This(A)> { typedef A type;};\
-   static const int arity =1;\
-   template<typename A> typename result<FNT##_wrap(A)>::type operator()( A const & a) const { return FNT(a);}\
- };\
- template <typename A>\
-  typename boost::result_of<map_impl<FNT##_wrap>(A)>::type FNT(A const & a) { return map(FNT##_wrap())(a); }
-
-#define TRIQS_ARRAYS_MATH_FNT1 (cos)(sin)(tan)(cosh)(sinh)(tanh)(acos)(asin)(atan)(exp)(log)(sqrt)(floor)
+ struct __triqs_##FNT##_wrap { template<typename A> auto operator()(A const & a) const DECL_AND_RETURN(FNT(a)); };\
+ template <typename A> \
+ typename boost::lazy_enable_if_c<ImmutableCuboidArray<A>::value,std::result_of<map_impl<__triqs_##FNT##_wrap,1>(A)>>::type\
+ FNT(A && a) { return map(__triqs_##FNT##_wrap{})(std::forward<A>(a)); }
+ 
+#define TRIQS_ARRAYS_MATH_FNT (abs)(real)(imag)(floor)
 
 #define AUX(r, data, elem) MAP_IT(elem)
- BOOST_PP_SEQ_FOR_EACH(AUX , nil , TRIQS_ARRAYS_MATH_FNT1);
+ BOOST_PP_SEQ_FOR_EACH(AUX , nil , TRIQS_ARRAYS_MATH_FNT);
 #undef AUX  
 #undef MAP_IT
+#undef TRIQS_ARRAYS_MATH_FNT
 
-
- using std::real;
- struct real_wrap { 
-   template<typename Sig> struct result;
-   template<typename This, typename A> struct result<This(A)> { typedef A type;};
-   template<typename This, typename A> struct result<This(std::complex<A>)> { typedef A type;};
-   static const int arity =1;
-   //   template<typename A> auto operator()( A const & a) const -> decltype(real(a)) { return real(a);}
-   template<typename A> typename result<real_wrap(A)>::type operator()( A const & a) const { return real(a);}
- };
-
- template <typename A>
-  //auto real(A const & a) -> decltype( map(real_wrap())(a)) { return map(real_wrap())(a); }
-  //  typename std::result_of<map(real_wrap)(A)>::type real(A const & a) { return map(real_wrap())(a); }
-  typename boost::result_of<map_impl<real_wrap>(A)>::type real(A const & a) { return map(real_wrap())(a); }
+// Functions only defined for vector, array but NOT for matrix. They act element wise.
+#define MAP_IT(FNT) \
+ using std::FNT;\
+ struct __triqs_##FNT##_wrap { template<typename A> auto operator()(A const & a) const DECL_AND_RETURN(FNT(a)); };\
+ template <typename A> \
+ typename boost::lazy_enable_if_c<ImmutableArray<A>::value||ImmutableVector<A>::value,std::result_of<map_impl<__triqs_##FNT##_wrap,1>(A)>>::type\
+ FNT(A && a) { return map(__triqs_##FNT##_wrap{})(std::forward<A>(a)); }
  
- using std::imag;
- struct imag_wrap { 
-   template<typename Sig> struct result;
-   template<typename This, typename A> struct result<This(A)> { typedef A type;};
-   template<typename This, typename A> struct result<This(std::complex<A>)> { typedef A type;};
-   static const int arity =1;
-   //   template<typename A> auto operator()( A const & a) const -> decltype(imag(a)) { return imag(a);}
-   template<typename A> typename result<imag_wrap(A)>::type operator()( A const & a) const { return imag(a);}
- };
-
- template <typename A>
-  //auto imag(A const & a) -> decltype( map(imag_wrap())(a)) { return map(imag_wrap())(a); }
-  //  typename std::result_of<map(imag_wrap)(A)>::type imag(A const & a) { return map(imag_wrap())(a); }
-  typename boost::result_of<map_impl<imag_wrap>(A)>::type imag(A const & a) { return map(imag_wrap())(a); }
+#define TRIQS_ARRAYS_MATH_FNT (exp)(cos)(sin)(tan)(cosh)(sinh)(tanh)(acos)(asin)(atan)(log)(sqrt)
 
+#define AUX(r, data, elem) MAP_IT(elem)
+ BOOST_PP_SEQ_FOR_EACH(AUX , nil , TRIQS_ARRAYS_MATH_FNT);
+#undef AUX  
+#undef MAP_IT
+#undef TRIQS_ARRAYS_MATH_FNT
 
 }}//namespace triqs::arrays 
 #endif
diff --git a/triqs/utility/concept_tools.hpp b/triqs/utility/concept_tools.hpp
index 3d20100b..99b9ccb6 100644
--- a/triqs/utility/concept_tools.hpp
+++ b/triqs/utility/concept_tools.hpp
@@ -31,13 +31,18 @@
 
 #define TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT(MyBeautifulConcept) \
  struct TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept) {};\
- template<typename T> struct MyBeautifulConcept : boost::is_base_of<TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept) , T> {};
+ template<typename T> struct MyBeautifulConcept : boost::is_base_of<TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept) , T> {};\
+ template<typename T> struct MyBeautifulConcept<T&> : boost::is_base_of<TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept), T> {};\
+ template<typename T> struct MyBeautifulConcept<T&&> : boost::is_base_of<TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept), T> {};
+
 
 #define TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT_R_AUX(r, data, i, elem) BOOST_PP_COMMA_IF(i) TRIQS_CONCEPT_TAG_NAME(elem) 
 
 #define TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT_R(MyBeautifulConcept,Rs) \
  struct TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept) : BOOST_PP_SEQ_FOR_EACH_I (TRIQS_DEFINE_CONCEPT_AND_ASSOCIATED_TRAIT_R_AUX,nil,Rs) {};\
- template<typename T> struct MyBeautifulConcept : boost::is_base_of<TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept), T> {};
+ template<typename T> struct MyBeautifulConcept : boost::is_base_of<TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept), T> {};\
+ template<typename T> struct MyBeautifulConcept<T&> : boost::is_base_of<TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept), T> {};\
+ template<typename T> struct MyBeautifulConcept<T&&> : boost::is_base_of<TRIQS_CONCEPT_TAG_NAME(MyBeautifulConcept), T> {};
 
 #ifdef TRIQS_COMPILE_TIME_DEBUG
 #define TRIQS_ASSERT_MODEL_CONCEPT(MyBeautifulConcept,T)  BOOST_CONCEPT_ASSERT((BCC_##MyBeautifulConcept<T>));