void inverse_fourier_impl  (gf_view<imtime,matrix_valued> gt , gf_view<imfreq,matrix_valued> const gw, matrix_valued){
  impl_worker w;
  for (size_t n1=0; n1<gw.data().shape()[1];n1++)
   for (size_t n2=0; n2<gw.data().shape()[2];n2++){
    auto gt_sl=slice_target_to_scalar(gt, n1, n2);
    auto gw_sl=slice_target_to_scalar(gw, n1, n2);
    w.inverse( gt_sl, gw_sl);
   }
 }

 void fit_tail(gf_view<imfreq> gf, tail_view known_moments, int n_moments, int n_min, int n_max,
   bool replace_by_fit ) {
  if (get_target_shape(gf) != known_moments.shape()) TRIQS_RUNTIME_ERROR << "shape of tail does not match shape of gf";
  gf.singularity() = fit_tail_impl(gf, known_moments, n_moments, n_min, n_max);
  if (replace_by_fit) { // replace data in the fitting range by the values from the fitted tail
   int i = 0;
   for (auto iw : gf.mesh()) { // (arrays::range(n_min,n_max+1)) {
    if (i >= n_min) gf[iw] = evaluate(gf.singularity(),iw);
    i++;
   }
   }
  }

 // compute a tail from the Legendre GF
 // this is Eq. 8 of our paper
 local::tail_view get_tail(gf_view<legendre> const & gl, int size = 10, int omin = -1) {

   auto sh = gl.data().shape().front_pop();
   local::tail t(sh, size, omin);
   t.data() = 0.0;

   for (int p=1; p<=t.order_max(); p++)
     for (auto l : gl.mesh())
       t(p) += (triqs::utility::legendre_t(l.index(),p)/pow(gl.domain().beta,p)) * gl[l];

   return t;

 }

 void legendre_matsubara_direct(gf_view<imfreq> gw, gf_const_view<legendre> gl) {

  gw() = 0.0;
  triqs::arrays::range R;

  // Use the transformation matrix
  for (auto om : gw.mesh()) {
   for (auto l : gl.mesh()) {
    gw[om] += legendre_T(om.index(), l.index()) * gl[l];
   }
  }

  gw.singularity() = get_tail(gl, gw.singularity().size(), gw.singularity().order_min());
 }

 tqa::matrix<double> density( gf_view<legendre> const & gl) { 

   auto sh = gl.data().shape().front_pop();
   tqa::matrix<double> res(sh);
   res() = 0.0;

   for (auto l : gl.mesh()) {
     res -= sqrt(2*l.index()+1) * gl[l];
   }
   res /= gl.domain().beta;

   return res;

 }

 // Impose a discontinuity G(\tau=0)-G(\tau=\beta)
 void enforce_discontinuity(gf_view<legendre> & gl, arrays::array_view<double,2> disc) {

   double norm = 0.0;
   arrays::vector<double> t(gl.data().shape()[0]);
   for (int i=0; i<t.size(); ++i) {
     t(i) = triqs::utility::legendre_t(i,1) / gl.domain().beta;
     norm += t(i)*t(i);
   }

   arrays::array<dcomplex, 2> corr(disc.shape());
   corr() = 0;
   for (auto const &l : gl.mesh()) corr += t(l.index()) * gl[l];

   auto _ = arrays::range{};
   for (auto const& l : gl.mesh()) gl.data()(l.index(), _, _) += (disc - corr) * t(l.index()) / norm;
 }

 auto curry_impl(gf_view<cartesian_product<Ms...>, Target, Singularity, Evaluator, IsConst> g) {
  // pick up the meshed corresponding to the curryed variables
  auto meshes_tuple = triqs::tuple::filter<pos...>(g.mesh().components());
  using var_t = cart_prod<triqs::tuple::filter_t<std::tuple<Ms...>, pos...>>;
  auto m = triqs::tuple::apply_construct<gf_mesh<var_t>>(meshes_tuple);
  auto l = [g](auto&&... x) { return partial_eval_linear_index<pos...>(g, std::make_tuple(x...)); };
  return make_gf_view_lambda_valued<var_t>(m, l);
 };

 static auto invoke(gf_view<cartesian_product<Ms...>, Target, Singularity, Evaluator, IsConst> g, XTuple const& x_tuple) {
  using var_t = cart_prod<triqs::tuple::filter_out_t<std::tuple<Ms...>, pos...>>;
  // meshes of the returned gf_view : just drop the mesh of the evaluated variables
  auto meshes_tuple_partial = triqs::tuple::filter_out<pos...>(g.mesh().components());
  // The mesh of the resulting function
  auto m = triqs::tuple::apply_construct<gf_mesh<var_t>>(meshes_tuple_partial);
  // rebuild a tuple of the size sizeof...(Ms), containing the linear indices and range at the position of evaluated variables.
  auto arr_args = triqs::tuple::inverse_filter<sizeof...(Ms), pos...>(x_tuple, arrays::range());
  // from it, we make a slice of the array of g, corresponding to the data of the returned gf_view
  auto arr2 = triqs::tuple::apply(g.data(), std::tuple_cat(arr_args, std::make_tuple(arrays::ellipsis{})));
  // We now also partial_eval the singularity
  auto singv = partial_eval_linear_index<pos...>(g.singularity(), x_tuple);
  using r_sing_t = typename decltype(singv)::regular_type;
  // finally, we build the view on this data.
  using r_t = gf_view<var_t, Target, r_sing_t, void, IsConst>;
  return r_t{m, arr2, singv, {}, {}};
 }

 void legendre_matsubara_inverse(gf_view<legendre> gl, gf_const_view<imtime> gt) {

  gl() = 0.0;
  legendre_generator L;
  auto N = gt.mesh().size() - 1;
  double coef;

  // Do the integral over imaginary time
  for (auto t : gt.mesh()) {
   if (t.index()==0 || t.index()==N) coef = 0.5;
   else coef = 1.0;
   L.reset(2 * t / gt.domain().beta - 1);
   for (auto l : gl.mesh()) {
    gl[l] += coef * sqrt(2 * l.index() + 1) * L.next() * gt[t];
   }
  }
  gl.data() *= gt.mesh().delta();
 }

 // Impose a discontinuity G(\tau=0)-G(\tau=\beta)
 void enforce_discontinuity(gf_view<legendre> & gl, tqa::array_view<double,2> disc) {

   double norm = 0.0;
   tqa::vector<double> t(gl.data().shape()[0]);
   for (int i=0; i<t.size(); ++i) {
     t(i) = triqs::utility::legendre_t(i,1) / gl.domain().beta;
     norm += t(i)*t(i);
   }

   tqa::array<double,2> corr(disc.shape()); corr() = 0;
   for (auto l : gl.mesh()) {
     corr += t(l.index()) * gl[l];
   }

   tqa::range R;
   for (auto l : gl.mesh()) {
     gl.data()(l.index(),R,R) += (disc - corr) * t(l.index()) / norm;
   }

 }

 void pade (gf_view<refreq> &gr, gf_view<imfreq> const &gw, int n_points, double freq_offset) {

  // make sure the GFs have the same structure
  //assert(gw.shape() == gr.shape());

  // copy the tail. it doesn't need to conform to the pade approximant
  gr.singularity() = gw.singularity();

  auto sh = gw.data().shape().front_pop();
  int N1 = sh[0], N2 = sh[1];
  for (int n1=0; n1<N1; n1++) {
    for (int n2=0; n2<N2; n2++) {

      arrays::vector<dcomplex> z_in(n_points); // complex points
      arrays::vector<dcomplex> u_in(n_points); // values at these points
      arrays::vector<dcomplex> a(n_points);    // corresponding Pade coefficients

      for (int i=0; i < n_points; ++i) z_in(i) = gw.mesh()[i];
      for (int i=0; i < n_points; ++i) u_in(i) = gw.on_mesh(i)(n1,n2);

      triqs::utility::pade_approximant PA(z_in,u_in);

      gr() = 0.0;
      for (auto om : gr.mesh()) {
        dcomplex e = om + dcomplex(0.0,1.0)*freq_offset;
        gr[om](n1,n2) = PA(e);
      }

    }
  }

 }

 template <typename M, typename T, typename S, typename E> gf<M, real_target_t<T>, S> real(gf_view<M, T, S, E> g) {
  return {g.mesh(), real(g.data()), g.singularity(), g.symmetry(), {}, {}};
 }

void fill_data(gf_view<legendre> g_l, int i, vector<double> const& data) {
 g_l.data()(range(),i,i) = data;
}

template <typename T> size_t n_blocks(gf_view<block_index, T> const &g) { return g.mesh().size(); }

 tail fit_tail_impl(gf_view<imfreq> gf, const tail_view known_moments, int n_moments, int n_min, int n_max) {

  // precondition : check that n_max is not too large
  n_max = std::min(n_max, int(gf.mesh().size()-1));

  tail res(get_target_shape(gf));
  if (known_moments.size())
   for (int i = known_moments.order_min(); i <= known_moments.order_max(); i++) res(i) = known_moments(i);

  // if known_moments.size()==0, the lowest order to be obtained from the fit is determined by order_min in known_moments
  // if known_moments.size()==0, the lowest order is the one following order_max in known_moments

  int n_unknown_moments = n_moments - known_moments.size();
  if (n_unknown_moments < 1) return known_moments;

  // get the number of even unknown moments: it is n_unknown_moments/2+1 if the first
  // moment is even and n_moments is odd; n_unknown_moments/2 otherwise
  int omin = known_moments.size() == 0 ? known_moments.order_min() : known_moments.order_max() + 1; // smallest unknown moment
  int omin_even = omin % 2 == 0 ? omin : omin + 1;
  int omin_odd = omin % 2 != 0 ? omin : omin + 1;
  int size_even = n_unknown_moments / 2;
  if (n_unknown_moments % 2 != 0 && omin % 2 == 0) size_even += 1;
  int size_odd = n_unknown_moments - size_even;

  int size1 = n_max - n_min + 1;
  // size2 is the number of moments

  arrays::matrix<double> A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT);
  arrays::matrix<double> B(size1, 1, FORTRAN_LAYOUT);
  arrays::vector<double> S(std::max(size_even, size_odd));
  const double rcond = 0.0;
  int rank;

  for (int i = 0; i < get_target_shape(gf)[0]; i++) {
   for (int j = 0; j < get_target_shape(gf)[1]; j++) {

    // fit the odd moments
    S.resize(size_odd);
    A.resize(size1,size_odd); //when resizing, gelss segfaults
    for (int k = 0; k < size1; k++) {
     auto n = n_min + k;
     auto iw = std::complex<double>(gf.mesh().index_to_point(n));

     B(k, 0) = imag(gf.data()(gf.mesh().index_to_linear(n), i, j));
     // subtract known tail if present
     if (known_moments.size() > 0)
      B(k, 0) -= imag(evaluate(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)), iw)(0, 0));

     for (int l = 0; l < size_odd; l++) {
      int order = omin_odd + 2 * l;
      A(k, l) = imag(pow(iw, -1.0 * order)); // set design matrix for odd moments
     }
    }

    arrays::lapack::gelss(A, B, S, rcond, rank);
    for (int m = 0; m < size_odd; m++) {
     res(omin_odd + 2 * m)(i, j) = B(m, 0);
    }

    // fit the even moments
    S.resize(size_even);
    A.resize(size1,size_even); //when resizing, gelss segfaults
    for (int k = 0; k < size1; k++) {
     auto n = n_min + k;
     auto iw = std::complex<double>(gf.mesh().index_to_point(n));

     B(k, 0) = real(gf.data()(gf.mesh().index_to_linear(n), i, j));
     // subtract known tail if present
     if (known_moments.size() > 0)
      B(k, 0) -= real(evaluate(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)), iw)(0, 0));

     for (int l = 0; l < size_even; l++) {
      int order = omin_even + 2 * l;
      A(k, l) = real(pow(iw, -1.0 * order)); // set design matrix for odd moments
     }
    }

    arrays::lapack::gelss(A, B, S, rcond, rank);
    for (int m = 0; m < size_even; m++) {
     res(omin_even + 2 * m)(i, j) = B(m, 0);
    }
   }
  }
  res.mask()()=n_moments;
  return res; // return tail
 }

 void inverse(gf_view<imtime,scalar_valued> gt, gf_view<imfreq,scalar_valued> const gw){
  using namespace impl_local_matsubara;
  static bool Green_Function_Are_Complex_in_time = false;
  // If the Green function are NOT complex, then one use the symmetry property
  // fold the sum and get a factor 2
  auto ta = gw(freq_infty());
  //TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO
  dcomplex d= ta(1)(0,0), A= ta.get_or_zero(2)(0,0), B = ta.get_or_zero(3)(0,0);
  double b1, b2, b3;
  dcomplex a1, a2, a3;
  
  double beta=gw.domain().beta;
  size_t L= gt.mesh().size() - ( gt.mesh().kind() == full_bins ? 1 : 0); //L can be different from gt.mesh().size() (depending on the mesh kind) and is given to the FFT algorithm 
  dcomplex iomega = dcomplex(0.0,1.0) * std::acos(-1) / beta;
  dcomplex iomega2 = -iomega * 2 * gt.mesh().delta() * (gt.mesh().kind() == half_bins ? 0.5 : 0.0) ; 
  double fact = (Green_Function_Are_Complex_in_time ? 1 : 2)/beta;
  g_in.resize( gw.mesh().size());
  g_out.resize(gt.mesh().size());
  
  if (gw.domain().statistic == Fermion){
   b1 = 0; b2 =1; b3 =-1; 
   a1 = d-B; a2 = (A+B)/2; a3 = (B-A)/2; 
  }
  else { 
   b1 = -0.5; b2 =-1; b3 =1; 
   a1=4*(d-B)/3; a2=B-(d+A)/2; a3=d/6+A/2+B/3;
  }
  g_in() = 0;
  for (auto & w: gw.mesh()) {
   g_in[ w.index() ] =  fact * exp(w.index()*iomega2) * ( gw[w] - (a1/(w-b1) + a2/(w-b2) + a3/(w-b3)) );
  }
  // for bosons GF(w=0) is divided by 2 to avoid counting it twice
  if (gw.domain().statistic == Boson && !Green_Function_Are_Complex_in_time ) g_in(0) *= 0.5;
  
  details::fourier_base(g_in, g_out, L, false);
  
  // CORRECT FOR COMPLEX G(tau) !!!
  typedef double gt_result_type;
  //typedef typename gf<imtime>::mesh_type::gf_result_type gt_result_type;
  if (gw.domain().statistic == Fermion){
   for (auto & t : gt.mesh()){
    gt[t] = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() ) * exp(-iomega*t)
    + oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) );
   }
  }
  else {
   for (auto & t : gt.mesh()) 
    gt[t] = convert_green<gt_result_type> ( g_out( t.index() == L ? 0 : t.index() )
    + oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) );
  }
  if (gt.mesh().kind() == full_bins)
   gt.on_mesh(L) = -gt.on_mesh(0)-convert_green<gt_result_type>(ta(1)(0,0));
  // set tail
   gt.singularity() = gw.singularity();
 }

 void direct (gf_view<imfreq,scalar_valued> gw, gf_view<imtime,scalar_valued> const gt) {
  using namespace impl_local_matsubara;
  auto ta = gt(freq_infty());
  //TO BE MODIFIED AFTER SCALAR IMPLEMENTATION TODO
  dcomplex d= ta(1)(0,0), A= ta.get_or_zero(2)(0,0), B = ta.get_or_zero(3)(0,0);
  double b1=0, b2=0, b3=0;
  dcomplex a1, a2, a3;
  double beta=gt.mesh().domain().beta;
  auto L = ( gt.mesh().kind() == full_bins ? gt.mesh().size()-1 : gt.mesh().size());
  double fact= beta/ gt.mesh().size();
  dcomplex iomega = dcomplex(0.0,1.0) * std::acos(-1) / beta;
  dcomplex iomega2 = iomega * 2 * gt.mesh().delta() * ( gt.mesh().kind() == half_bins ? 0.5 : 0.0);
  g_in.resize(gt.mesh().size());
  g_out.resize(gw.mesh().size());
  if (gw.domain().statistic == Fermion){
   b1 = 0; b2 =1; b3 =-1; 
   a1 = d-B; a2 = (A+B)/2; a3 = (B-A)/2; 
  }
  else { 
   b1 = -0.5; b2 =-1; b3 =1; 
   a1 = 4*(d-B)/3; a2 = B-(d+A)/2; a3 = d/6+A/2+B/3; 
  }
  if (gw.domain().statistic == Fermion){
   for (auto & t : gt.mesh())  
    g_in[t.index()] = fact * exp(iomega*t) * ( gt[t] - ( oneFermion(a1,b1,t,beta) + oneFermion(a2,b2,t,beta)+ oneFermion(a3,b3,t,beta) ) );
  }
  else {
   for (auto & t : gt.mesh()) 
    g_in[t.index()] =  fact * ( gt[t] -  ( oneBoson(a1,b1,t,beta) + oneBoson(a2,b2,t,beta) + oneBoson(a3,b3,t,beta) ) );  
  }
  details::fourier_base(g_in, g_out, L, true);
  for (auto & w : gw.mesh()) {
   gw[w] = g_out( w.index() ) * exp(iomega2*w.index() ) + a1/(w-b1) + a2/(w-b2) + a3/(w-b3); 
  }
  gw.singularity() = gt.singularity();// set tail
 }

 std14::enable_if_t<!arrays::is_scalar<RHS>::value> triqs_gf_view_assign_delegation(gf_view<M, T, S, E> g, RHS const &rhs) {
  if (!(g.mesh() == rhs.mesh()))
   TRIQS_RUNTIME_ERROR << "Gf Assignment in View : incompatible mesh" << g.mesh() << " vs " << rhs.mesh();
  for (auto const &w : g.mesh()) g[w] = rhs[w];
  g.singularity() = rhs.singularity();
 }

 template <typename F, typename G> gf<block_index, std14::result_of_t<F(G)>> map(F &&f, gf_view<block_index, G> g) {
  return make_block_gf(get_block_names(g), _map(f, g.data()));
 }

 void triqs_clef_auto_assign(gf_view<M, T, S, E> g, RHS const &rhs) {
  triqs_clef_auto_assign_impl(g, rhs, typename std::is_base_of<tag::composite, gf_mesh<M>>::type());
  assign_singularity_from_function(g.singularity(), rhs);
  // access to the data . Beware, we view it as a *matrix* NOT an array... (crucial for assignment to scalars !)
  // if f is an expression, replace the placeholder with a simple tail.
 }