Work on gf

gf: security in the case beta<0 added in the mesh construction
gf: inline added in slice
test/triqs/gf: test of on_mesh() added
gfs: scalar for two-real_times
test/triqs/gf/ renamed in gfs, test gf_retw.cpp completed
gfs: evaluator homogeneised
two_times: evaluator corrected
test/triqs/gf/ renamed in gfs, test gf_retw.cpp completed

+ Correction after rebase

Fix a test : gf_re_im_freq_time
There is an issue with the last point.
To be fixed.

test/triqs/gfs/gf_re_im_freq_time.cpp

@@ -0,0 +1,82 @@
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+
+#include <triqs/gfs/re_im_freq.hpp> 
+#include <triqs/gfs/re_im_time.hpp> 
+#include <triqs/gfs/refreq_imtime.hpp> 
+
+#include <triqs/gfs/local/fourier_real.hpp> 
+#include <triqs/arrays.hpp>
+
+namespace tql= triqs::clef;
+using namespace triqs::gfs;
+
+int main() {
+
+ double precision=10e-9;
+ double beta =1.;
+  
+ double tmin=0.;
+ double tmax=1.0;
+ int n_re_time=100;
+ int n_im_time=100;
+ 
+ double wmin=0.;
+ double wmax=1.0;
+ 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_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);
+ 
+ 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_tau(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];
+
+ //identical functions
+ G_w(w_) << exp( -2*tau ) / (w_*w_ + 1 );
+ //the singularity must be removed as it is inexistent in re_im_time, to give the same TF. 
+ 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_mesh_imtime(G_w_tau, index);
+ for(auto& w:G_w.mesh())
+  if ( std::abs(G_w(w)-G_w2(w)) > precision) TRIQS_RUNTIME_ERROR<<" fourier_slice error : w="<< w <<" ,G_w="<< G_w(w)<<" ,G_w2="<< G_w2(w) <<"\n";
+ 
+ //test of the interpolation
+ std::cout << G_t_tau(0.789,0.123) << std::endl; 
+ std::cout << G_w_wn( 0.789,0.123) << std::endl; 
+ std::cout << G_w_tau(0.789,0.123) << std::endl; 
+ 
+ //test of on_mesh()
+ std::cout << G_w_tau.on_mesh(n_re_freq/2,n_im_time/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_tau", G_w_tau);
+
+ 
+ // try to slice it
+ auto gt = slice_mesh_imtime(G_t_tau, 1);
+ h5_write(file, "gt0", gt);
+ auto gw = slice_mesh_imtime(G_w_tau, 1);
+ h5_write(file, "gw0", gw);
+
+ //comparison of the TF of the one time and sliced two times GF's
+ auto G_t = inverse_fourier(G_w);
+ 
+ auto G_t2 = inverse_fourier(slice_mesh_imtime(G_w_tau, index) );
+ for(auto& t:G_t.mesh())
+  { 
+   // BUG HERE the last point is badly rounded
+   if ( (t< G_t.mesh().size()-1) && (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";  
+  }
+}

test/triqs/gfs/gf_retw.cpp

@@ -2,16 +2,19 @@
 
 #include <triqs/gfs/refreq.hpp> 
 #include <triqs/gfs/retime.hpp> 
+#include <triqs/gfs/imfreq.hpp> 
+#include <triqs/gfs/imtime.hpp> 
 
-namespace tql= triqs::clef;
-namespace tqa= triqs::arrays;
-using tqa::range;
-using triqs::arrays::make_shape;
-using triqs::gfs::scalar_valued;
-using triqs::gfs::Fermion;
 using triqs::gfs::refreq;
 using triqs::gfs::retime;
+using triqs::gfs::imfreq;
+using triqs::gfs::imtime;
 using triqs::gfs::make_gf;
+using triqs::gfs::scalar_valued;
+using triqs::gfs::Fermion;
+using triqs::arrays::make_shape;
+using triqs::arrays::range;
+double precision=10e-12;
 
 #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl;
 
@@ -19,29 +22,56 @@ int main() {
 
  double beta =1;
  double tmax=10;
+ double wmax=10;
  int N=1000;
- auto Gw = make_gf<refreq> (-10, 10, N, make_shape(2,2));
+ 
+ auto Gw = make_gf<refreq> (-wmax, wmax, N, make_shape(2,2));
  auto Gt = make_gf<retime> (0, tmax, N, make_shape(2,2));
- auto Gw2 = make_gf<refreq,scalar_valued> (-10, 10, N);
+ auto Gw2 = make_gf<refreq,scalar_valued> (-wmax, wmax, N);
  auto Gt2 = make_gf<retime,scalar_valued> (0, tmax, N);
+ auto Giw = make_gf<imfreq> (beta, Fermion, make_shape(2,2), N);
+ auto Git = make_gf<imtime> (beta, Fermion, make_shape(2,2), N);
+ auto Giw2 = make_gf<imfreq,scalar_valued> (beta, Fermion, N);
+ auto Git2 = make_gf<imtime,scalar_valued> (beta, Fermion, N);
 
  int i =0;
  for (auto & t : Gt.mesh()) Gt[t] = 1.0*t;
  for (auto & w : Gw.mesh()) Gw[w] = 1.0*w;
+ for (auto & it : Git.mesh()) Git[it] = 1.0*it;
+ for (auto & iw : Giw.mesh()) Giw[iw] = 1.0*iw;
  
  triqs::clef::placeholder<0> t_;
  triqs::clef::placeholder<1> w_;
  Gt2( t_) << 1.0*t_;
  Gw2( w_) << 1.0*w_;
+ Git2( t_) << 1.0*t_;
+ Giw2( w_) << 1.0*w_;
  
- std::cout << Gt(3.255)(0,0) << std::endl; 
+ std::cout << "Gt(0.3352*tmax)(0,0)=" << Gt(0.3352*tmax)(0,0) << std::endl;
- std::cout << Gt2(3.255) << std::endl; 
+ if ( std::abs(Gt(0.3352*tmax)(0,0)-0.3352*tmax) > precision) TRIQS_RUNTIME_ERROR<< "evaluator error in Gt\n";  
- std::cout << Gw(3.255)(0,0) << std::endl; 
+ std::cout << "Gt2(0.3352*tmax)=" << Gt2(0.3352*tmax) << std::endl;
- std::cout << Gw2(3.255) << std::endl; 
+ if ( std::abs(Gt2(0.3352*tmax)-0.3352*tmax) > precision) TRIQS_RUNTIME_ERROR<< "evaluator error in Gt2\n";  
+ std::cout << "Gw(0.3352*wmax)(0,0)=" << Gw(0.3352*wmax)(0,0) << std::endl;
+ if ( std::abs(Gw(0.3352*wmax)(0,0)-0.3352*wmax) > precision) TRIQS_RUNTIME_ERROR<< "evaluator error in Gw\n";  
+ std::cout << "Gw2(0.3352*wmax)=" << Gw2(0.3352*wmax) << std::endl;
+ if ( std::abs(Gw2(0.3352*wmax)-0.3352*wmax) > precision) TRIQS_RUNTIME_ERROR<< "evaluator error in Gw2\n";  
+ 
+ std::cout << "Git(0.3352*beta)(0,0)=" << Git(0.3352*beta)(0,0) << std::endl;
+ if ( std::abs(Git(0.3352*beta)(0,0)-0.3352*beta) > precision) TRIQS_RUNTIME_ERROR<< "evaluator error in Git\n"; 
+ std::cout << "Git2(0.3352*beta)=" << Git2(0.3352*beta) << std::endl; 
+ if ( std::abs(Git2(0.3352*beta)-0.3352*beta) > precision) TRIQS_RUNTIME_ERROR<< "evaluator error in Git2\n";   
+ 
+ if ( std::abs(Gw2.on_mesh(N/3) -Gw2[N/3] ) > precision) TRIQS_RUNTIME_ERROR<< "error in on_mesh()\n";  
+ if ( std::abs(Gt2.on_mesh(N/3) -Gt2[N/3] ) > precision) TRIQS_RUNTIME_ERROR<< "error in on_mesh()\n";  
+ if ( std::abs(Git2.on_mesh(N/3)-Git2[N/3]) > precision) TRIQS_RUNTIME_ERROR<< "error in on_mesh()\n";  
  
  // test hdf5 
  H5::H5File file("ess_gfre.h5", H5F_ACC_TRUNC );
  h5_write(file, "gt",  Gt);
  h5_write(file, "gw",  Gw);
+ h5_write(file, "git",  Git);
+ h5_write(file, "giw",  Giw);
+ 
+ 
 
 }

triqs/gfs/domains/matsubara.hpp

@@ -30,7 +30,9 @@ namespace triqs { namespace gfs {
   typedef typename mpl::if_c<IsComplex, std::complex<double>, double>::type point_t;
   double beta;
   statistic_enum statistic;
-  matsubara_domain (double Beta=1, statistic_enum s = Fermion): beta(Beta), statistic(s){}
+  matsubara_domain (double Beta=1, statistic_enum s = Fermion): beta(Beta), statistic(s){
+    if(beta<0)TRIQS_RUNTIME_ERROR<<"Matsubara domain construction :  beta <0 : beta ="<< beta <<"\n";
+  }
   bool operator == (matsubara_domain const & D) const { return ((std::abs(beta - D.beta)<1.e-15) && (statistic == D.statistic));}
 
   /// Write into HDF5

triqs/gfs/re_im_time.hpp

@@ -57,18 +57,16 @@ namespace triqs { namespace gfs {
     static constexpr int arity = 2;
     template<typename G>
      std::complex<double> operator() (G const * g, double t, double tau)  const {
-      //auto & data = g->data();
-      //auto & mesh = g->mesh();
       double beta = std::get<1>(g->mesh().components()).domain().beta;
       int p = std::floor(tau/beta);
       tau -= p*beta;
       size_t nr,ni; double wr,wi; bool in;
       std::tie(in, nr, wr) = windowing( std::get<0>(g->mesh().components()),t);
-      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds, tmax=" << std::get<0>(g->mesh().components()).x_max() << ", tmin=" << std::get<0>(g->mesh().components()).x_min() << "here, t=" <<t;
       std::tie(in, ni, wi) = windowing( std::get<1>(g->mesh().components()),tau);
-      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds, taumax=" << std::get<1>(g->mesh().components()).x_max()<< ", taumin=" << std::get<1>(g->mesh().components()).x_min() << "here, tau=" <<tau;
       auto gg = on_mesh(*g); //[g]( size_t nr, size_t ni) {return g->on_mesh(nr,ni);}; //data( g->mesh().index_to_linear(nr,ni));
-      auto res =  wr *( wi*gg(nr,ni) + (1-wi)*gg(nr,ni+1)) + (1-wr) * ( wi*gg(nr+1,ni) + (1-wi)*gg(nr+1,ni+1)); 
+      auto res =  (1-wr) * ( (1-wi) * gg(nr,ni) + wi * gg(nr,ni+1)) + wr * ( (1-wi) * gg(nr+1,ni) + wi * gg(nr+1,ni+1)); 
       return ((std::get<1>(g->mesh().components()).domain().statistic == Fermion) && (p%2==1) ? -res : res);
      } 
    };

triqs/gfs/two_real_times.hpp

@@ -1,5 +1,5 @@
 /*******************************************************************************
- *
+ * 
  * TRIQS: a Toolbox for Research in Interacting Quantum Systems
  *
  * Copyright (C) 2012 by M. Ferrero, O. Parcollet
@@ -26,7 +26,7 @@
 #include "./meshes/product.hpp"
 
 namespace triqs { namespace gfs { 
-
+ 
  struct two_real_times {};
 
  // the mesh
@@ -40,9 +40,6 @@ namespace triqs { namespace gfs {
 
  namespace gfs_implementation { 
 
-  // h5 name
-  template<typename Opt> struct h5_name<two_real_times,matrix_valued,Opt> { static std::string invoke(){ return  "GfTwoRealTime";}};
-
   /// ---------------------------  closest mesh point on the grid ---------------------------------
 
   template<typename Opt>
@@ -57,75 +54,62 @@ namespace triqs { namespace gfs {
 
    };
 
+   // h5 name
+  template<typename Opt> struct h5_name<two_real_times,matrix_valued,Opt> { static std::string invoke(){ return  "GfTwoRealTime";}};
+  template<typename Opt> struct h5_name<two_real_times,scalar_valued,Opt> { static std::string invoke(){ return  "GfTwoRealTime_s";}};
+  
   /// ---------------------------  evaluator ---------------------------------
-
+  
-  template<typename Opt>
+  template<typename Opt, typename Target>
-   struct evaluator<two_real_times,matrix_valued,Opt> {
+  struct evaluator<two_real_times,Target,Opt> {
-    static constexpr int arity = 2;
+   static constexpr int arity = 2;
-    template<typename G>
+   typedef typename std::conditional < std::is_same<Target, matrix_valued>::value, arrays::matrix<std::complex<double>>, std::complex<double>>::type rtype; 
-     arrays::matrix<std::complex<double> > operator() (G const * g, double t0, double t1)  const {
+   template<typename G>
-      auto & data = g->data();
+   rtype operator() (G const * g, double t0, double t1)  const {
-      auto & m = std::get<0>(g->mesh().components()); 
+    size_t n0,n1; double w0,w1; bool in;
-      size_t n0,n1; double w0,w1; bool in;
+    std::tie(in, n0, w0) = windowing(std::get<0>(g->mesh().components()),t0);
-      std::tie(in, n0, w0) = windowing(m,t0);
+    if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
-      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+    std::tie(in, n1, w1) = windowing(std::get<1>(g->mesh().components()),t1);
-      std::tie(in, n1, w1) = windowing(m,t1);
+    if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
-      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+    auto gg = on_mesh(*g);
-      auto gg = [g,data]( size_t n0, size_t n1) {return data(g->mesh().index_to_linear(std::tuple<size_t,size_t>{n0,n1}), arrays::ellipsis());};
+    return  (1-w0) * ( (1-w1) * gg(n0, n1) + w1 * gg(n0, n1+1) ) + w0 * ( (1-w1) * gg(n0+1, n1) + w1 * gg(n0+1, n1+1)); 
-      return  w0 * ( w1*gg(n0,n1) + (1-w1)*gg(n0,n1+1) ) + (1-w0) * ( w1*gg(n0+1,n1) + (1-w1)*gg(n0+1,n1+1)); 
+   } 
-     } 
+  };
-   };
+  
-
   /// ---------------------------  data access  ---------------------------------
-
+  
   template<typename Opt> struct data_proxy<two_real_times,matrix_valued,Opt> : data_proxy_array<std::complex<double>,3> {};
-
+  template<typename Opt> struct data_proxy<two_real_times,scalar_valued,Opt> : data_proxy_array<std::complex<double>,1> {};
+  
   // -------------------------------   Factories  --------------------------------------------------
-
+  
+  //matrix_valued
   template<typename Opt> struct factories<two_real_times, matrix_valued,Opt> {
    typedef gf<two_real_times, matrix_valued,Opt> gf_t;
    typedef gf_mesh<two_real_times, Opt> mesh_t;
-
    static gf_t make_gf(double tmax, double n_time_slices, tqa::mini_vector<size_t,2> shape) {
     auto m =  gf_mesh<two_real_times,Opt>(tmax, n_time_slices);
     typename gf_t::data_regular_t A(shape.front_append(m.size())); A() =0;
     return gf_t (m, std::move(A), nothing(), nothing() ) ;
    }
   };
-
+  
-  // -------------------------------   Path  --------------------------------------------------
+  //scalar_valued
-  /*
+  template<typename Opt> struct factories<two_real_times, scalar_valued,Opt> {
-     struct path { 
+   typedef gf<two_real_times, scalar_valued,Opt> gf_t;
-
+   typedef gf_mesh<two_real_times, Opt> mesh_t;
-     typedef typename mesh_t::index_t mesh_pt_t; 
+   
-     typedef triqs::arrays::mini_vector<long,2> delta_t;
+   static gf_t make_gf(double tmax, double n_time_slices) {
-
+    auto m =  gf_mesh<two_real_times,Opt>(tmax, n_time_slices);
-     delta_t pt, delta;
+    typename gf_t::data_non_view_t A(m.size()); A() =0;
-     size_t L;
+    return gf_t (m, std::move(A), nothing(), nothing() ) ;
-
+   }
-     path( mesh_t const & m, pt_t const & start_pt, delta_t const & d_) : pt(start_pt), delta(d_), L(std::get<1>(m.components()).size()){}
+  };
-
+  
-     void advance() { pt += delta;}
-
-     bool out_of_mesh () const { return (! ( (pt[1]>=0) && ( pt[0] >= pt[1]) && (pt[0]<= L)));} 
-
-     typedef mesh_pt_generator<path> iterator;
-     iterator begin() const { return {this, false};}
-     iterator end()   const { return {this, true};}
-
-     };
-
-     path make_path ( mesh_t const & m, typename mesh_t::index_t starting_point, delta) {
-     return path(m, starting_point,delta);
-     }
-
-  // for (auto & p : make_path(G.mesh(), make_tuple(i,j), make_tuple(di,dj) )) G(p) +=0;
-  */
-
  } // gfs_implementation
-
+ 
  // -------------------------------   Additionnal free function for this gf  --------------------------------------------------
-
+ 
  // from g(t,t') and t, return g(t-t') for any t'>t 
  // 
  gf<retime> slice (gf_view<two_real_times> const & g, double t) { 
@@ -133,7 +117,7 @@ namespace triqs { namespace gfs {
   long it = get_closest_mesh_pt_index(m, t);              //index of t on this mesh
   long nt = m.size() - it;
   if (it+1 < nt) nt = it+1 ;                              //nt=length of the resulting GF's mesh
-  double dt = m.delta();
+    double dt = m.delta();
   auto res = make_gf<retime>(0, 2*(nt-1)*dt, nt, g(t,t).shape());
   res() = 0;
   auto _ = arrays::range();// everyone
@@ -142,11 +126,11 @@ namespace triqs { namespace gfs {
   }
   return res;
  }
-
+ 
  // Get the 1 time mesh from the 2 times cartesian product (for cython interface mainly)
  template<typename M> 
-  auto get_1d_mesh_from_2times_mesh(M const & m) DECL_AND_RETURN(std::get<0>(m.components()));
+ auto get_1d_mesh_from_2times_mesh(M const & m) DECL_AND_RETURN(std::get<0>(m.components()));
-
+ 
 }}
 #endif