draft of multi-time gf (t,tau), etc...

test/triqs/gf/gf_re_im_freq_time.cpp

@@ -0,0 +1,64 @@
+#define TRIQS_ARRAYS_ENFORCE_BOUNDCHECK
+
+#include <triqs/gf/re_im_freq.hpp> 
+#include <triqs/gf/re_im_time.hpp> 
+#include <triqs/gf/refreq_imtime.hpp> 
+
+#include <triqs/gf/local/fourier_real.hpp> 
+#include <triqs/arrays.hpp>
+
+namespace tql= triqs::clef;
+using triqs::gf::Fermion;
+using triqs::gf::make_gf;
+using triqs::gf::retime;
+using triqs::gf::refreq;
+using triqs::gf::refreq_imtime;
+using triqs::gf::re_im_time;
+using triqs::gf::re_im_freq;
+using triqs::arrays::make_shape;
+using triqs::gf::scalar_valued;
+
+int main() {
+  
+ 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);
+ 
+ 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 );
+ 
+ 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 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);
+ std::cout  << gt.data()<< std::endl ;
+ h5_write(file, "gt0", gt);
+ auto gw = slice_mesh_imtime(G_w_tau, 1);
+ std::cout  << gw.data()<< std::endl ;
+ h5_write(file, "gw0", gw);
+
+
+}

triqs/gf/re_im_freq.hpp

@@ -0,0 +1,90 @@
+/*******************************************************************************
+ * 
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2012 by M. Ferrero, 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_RE_IM_FREQ_H
+#define TRIQS_GF_RE_IM_FREQ_H
+#include "./tools.hpp"
+#include "./gf.hpp"
+#include "./refreq.hpp"
+#include "./imfreq.hpp"
+#include "./meshes/product.hpp"
+
+namespace triqs { namespace gf { 
+  
+  struct re_im_freq {};
+  
+  namespace gf_implementation {
+    
+    // the mesh
+    template<typename Opt> struct mesh<re_im_freq,Opt>  { 
+      typedef typename mesh<refreq,Opt>::type m1_t; 
+      typedef typename mesh<imfreq,Opt>::type m2_t; 
+      typedef mesh_product<m1_t,m2_t> type;
+      static type make (double wmin, double wmax, size_t n_freq_re, double beta, statistic_enum S, size_t n_freq_im) { 
+        return {make_gf_mesh<refreq,Opt>(wmin,wmax,n_freq_re,full_bins), make_gf_mesh<imfreq,Opt>(beta, S, n_freq_im)};
+      }
+    };
+    
+    // singularity
+    template<typename Opt> struct singularity<re_im_freq,scalar_valued,Opt>  { typedef gf<refreq,scalar_valued> type;};
+    
+    // h5 name
+    template<typename Opt> struct h5_name<re_im_freq,scalar_valued,Opt> { static std::string invoke(){ return  "GfReImFreq";}};
+    
+    /// ---------------------------  data access  ---------------------------------
+    
+    template<typename Opt> struct data_proxy<re_im_freq,scalar_valued,Opt> : data_proxy_array<std::complex<double>,1> {};
+    
+    /// ---------------------------  evaluator ---------------------------------  
+    
+    template<typename Opt>
+    struct evaluator<re_im_freq,scalar_valued,Opt> {
+      static constexpr int arity = 2;
+      template<typename G>
+      std::complex<double> operator() (G const * g, double w, long n)  const {
+        auto & data = g->data();
+        auto & mesh = g->mesh();
+        size_t nr; double wr; bool in;
+        std::tie(in, nr, wr) = windowing( std::get<0>(g->mesh().components()), w);
+        if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+        auto gg = [g,data,mesh]( size_t nr, size_t n) {return data(mesh.index_to_linear(std::tuple<size_t,size_t>{nr,n}));};
+        return wr * gg(nr,n) + (1-wr) * gg(nr+1,n) ;
+      } 
+    };
+    
+    // -------------------------------   Factories  --------------------------------------------------
+    
+    template<typename Opt> struct factories<re_im_freq, scalar_valued,Opt> {
+      typedef gf<re_im_freq, scalar_valued,Opt> gf_t;
+//       typedef typename mesh<re_im_freq, Opt>::type mesh_t;
+      
+      static gf_t make_gf(double wmin, double wmax, size_t nw, double beta, statistic_enum S, size_t nwn) { 
+        auto m =  make_gf_mesh<re_im_freq,Opt>(wmin, wmax, nw, beta, S, nwn);
+        typename gf_t::data_non_view_t A(m.size()); 
+        A() =0;
+        return gf_t (m, std::move(A), triqs::gf::make_gf<refreq,scalar_valued>(wmin, wmax, nw), nothing() ) ;
+      }
+    };
+    
+} // gf_implementation
+
+}}
+#endif
+

triqs/gf/re_im_time.hpp

@@ -0,0 +1,116 @@
+/*******************************************************************************
+ * 
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2012 by M. Ferrero, 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_RE_IM_TIMES_H
+#define TRIQS_GF_RE_IM_TIMES_H
+#include "./tools.hpp"
+#include "./gf.hpp"
+#include "./retime.hpp"
+#include "./imtime.hpp"
+#include "./meshes/product.hpp"
+
+namespace triqs { namespace gf { 
+  
+ struct re_im_time {};
+
+ namespace gf_implementation { 
+
+  // the mesh
+  template<typename Opt> struct mesh<re_im_time,Opt>  { 
+   typedef typename mesh<retime,Opt>::type m1_t; 
+   typedef typename mesh<imtime,Opt>::type m2_t; 
+   typedef mesh_product<m1_t,m2_t> type;
+   static type make (double tmin, double tmax, size_t nt, double beta, statistic_enum S, size_t ntau, triqs::gf::mesh_kind mk=triqs::gf::full_bins) { 
+    return {make_gf_mesh<retime,Opt>(tmin,tmax,nt), make_gf_mesh<imtime,Opt>(beta,S, ntau, mk)};
+   }
+  };
+
+  // singularity
+  //template<typename Opt> struct singularity<re_im_time, scalar_valued, Opt>  { typedef gf<retime,scalar_valued> type;};
+
+  // h5 name
+  template<typename Opt> struct h5_name<re_im_time,scalar_valued,Opt> { static std::string invoke(){ return  "GfReImTime";}};
+
+  /// ---------------------------  data access  ---------------------------------
+
+  template<typename Opt> struct data_proxy<re_im_time,scalar_valued,Opt> : data_proxy_array<std::complex<double>,1> {};
+
+  /// ---------------------------  evaluator ---------------------------------
+  template<typename Opt>
+   struct evaluator<re_im_time,scalar_valued,Opt> {
+    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";
+      std::tie(in, ni, wi) = windowing( std::get<1>(g->mesh().components()),tau);
+      if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+      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)); 
+      return ((std::get<1>(g->mesh().components()).domain().statistic == Fermion) && (p%2==1) ? -res : res);
+     } 
+   };
+
+  // -------------------------------   Factories  --------------------------------------------------
+
+  template<typename Opt> struct factories<re_im_time, scalar_valued,Opt> {
+   typedef gf<re_im_time, scalar_valued,Opt> gf_t;
+   //       typedef typename mesh<re_im_time, Opt>::type mesh_t;
+
+   static gf_t make_gf(double tmin, double tmax, size_t nt, double beta, statistic_enum S, size_t ntau, mesh_kind mk=full_bins) { 
+    auto m =  make_gf_mesh<re_im_time,Opt>(tmin,tmax, nt, beta, S, ntau, mk);
+    typename gf_t::data_non_view_t A(m.size()); 
+    A() =0;
+    return gf_t (m, std::move(A), nothing(), nothing());
+    //return gf_t (m, std::move(A), triqs::gf::make_gf<retime,scalar_valued>(tmin, tmax, nt), nothing());
+   }
+  };
+
+ } // gf_implementation
+
+ // CHANGE THIS NAME !!! 
+ template<typename RHS, bool V, typename Variable, typename Target, typename Opt > 
+  void assign_from_expression (gf_impl<Variable,Target,Opt,V> const &, RHS) {}
+
+ //slices
+ gf_view<retime,scalar_valued> slice_mesh_imtime (gf_view<re_im_time,scalar_valued> g, size_t index) { 
+  auto arr = reinterpret_linear_array(g.mesh(),g.data()); // view it as a 2d array
+  return { std::get<0>(g.mesh().components()), arr(arrays::range(), index), local::tail(1,1), nothing() };
+ }
+
+ /* gf_view<imtime,scalar_valued> slice_mesh_retime ( gf_view<re_im_time,scalar_valued> g, size_t index) { 
+    auto arr = reinterpret_linear_array(g.mesh(),g.data()); // view it as a 2d array
+    return { std::get<1>(g.mesh().components()), arr(index, arrays::range()), g.singularity().singularity(), nothing() };
+    }
+    */
+ //   
+ //  gf_view<retime,scalar_valued>  slice_meshes ( gf_view<re_im_time,scalar_valued> g, size_t index) { 
+ //    return { std::get<0>(g.mesh().components()), g.data()(arrays::range(), index), tail ( g.singularity(.......) ), g.symmetry()}
+ //  }
+
+}}
+#endif
+

triqs/gf/refreq_imtime.hpp

@@ -0,0 +1,99 @@
+/*******************************************************************************
+ * 
+ * TRIQS: a Toolbox for Research in Interacting Quantum Systems
+ *
+ * Copyright (C) 2012 by M. Ferrero, 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_IMTIME_REFREQ_TIMES_H
+#define TRIQS_GF_IMTIME_REFREQ_TIMES_H
+#include "./tools.hpp"
+#include "./gf.hpp"
+#include "./imtime.hpp"
+#include "./refreq.hpp"
+#include "./meshes/product.hpp"
+
+namespace triqs { namespace gf {
+  
+  struct refreq_imtime {};
+  
+  namespace gf_implementation { 
+    
+    // the mesh
+    template<typename Opt> struct mesh<refreq_imtime,Opt>  {
+      typedef typename mesh<refreq,Opt>::type m1_t; 
+      typedef typename mesh<imtime,Opt>::type m2_t; 
+      typedef mesh_product<m1_t,m2_t> type;
+      static type make (double wmin, double wmax, size_t n_freq, double beta, statistic_enum S, size_t nt, mesh_kind mk=full_bins ) { 
+        return {make_gf_mesh<refreq,Opt>(wmin,wmax,n_freq), make_gf_mesh<imtime,Opt>(beta,S,nt,mk)};
+      }
+    };
+    
+    // singularity
+    //template<typename Opt> struct singularity<refreq_imtime,scalar_valued,Opt>  { typedef gf<retime,scalar_valued> type;};
+    
+    // h5 name
+    template<typename Opt> struct h5_name<refreq_imtime,scalar_valued,Opt> { static std::string invoke(){ return  "GfReFreqImTime";}};
+    
+    /// ---------------------------  data access  ---------------------------------
+    
+    template<typename Opt> struct data_proxy<refreq_imtime,scalar_valued,Opt> : data_proxy_array<std::complex<double>,1> {};
+    
+    /// ---------------------------  evaluator ---------------------------------
+    
+    template<typename Opt>
+    struct evaluator<refreq_imtime,scalar_valued,Opt> {
+      static constexpr int arity = 2;
+      template<typename G>
+      std::complex<double> operator() (G const * g, double omega, double tau)  const {
+       double beta = std::get<1>(g->mesh().components()).domain().beta;
+       int p = std::floor(tau/beta);
+       tau -= p*beta;
+       size_t n1,n2; double w1,w2; bool in;
+       std::tie(in, n1, w1) = windowing( std::get<0>(g->mesh().components()),omega);
+       if (!in) TRIQS_RUNTIME_ERROR <<" Evaluation out of bounds";
+       std::tie(in, n2, w2) = windowing( std::get<1>(g->mesh().components()),tau);
+       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));
+       return ((std::get<1>(g->mesh().components()).domain().statistic == Fermion) && (p%2==1) ? -res : res);
+      } 
+    };
+
+    // -------------------------------   Factories  --------------------------------------------------
+
+    template<typename Opt> struct factories<refreq_imtime, scalar_valued,Opt> {
+     typedef gf<refreq_imtime, scalar_valued,Opt> gf_t;
+     
+     static gf_t make_gf(double wmin, double wmax, size_t n_freq, double beta, statistic_enum S, size_t nt, mesh_kind mk=full_bins) { 
+      auto m =  make_gf_mesh<refreq_imtime,Opt>(wmin, wmax, n_freq, beta,S, nt, mk);
+      typename gf_t::data_non_view_t A(m.size()); 
+      A() =0;
+      return gf_t (m, std::move(A), nothing(), nothing() ) ;
+     }
+    };
+
+  } // gf_implementation
+
+  //slices
+ gf_view<refreq,scalar_valued> slice_mesh_imtime (gf_view<refreq_imtime,scalar_valued> g, size_t index) { 
+  auto arr = reinterpret_linear_array(g.mesh(),g.data()); // view it as a 2d array
+  return { std::get<0>(g.mesh().components()), arr(index,arrays::range()), local::tail(1,1), nothing() };
+ }
+
+}}
+#endif
+