Added basic doc + additional test for tail fitting

- added test for a 'real-life' GF + corresponding output
- added basic usage documentation for tail fitting from c++. Full implementation details yet to be written

doc/reference/c++/gf/set_tail_from_fit.rst

@@ -0,0 +1,67 @@
+Fitting the tail of a Matsubara Green's function
+#################################################
+
+API
+~~~~~
+
+The tail of a given ``gf<imfreq>/gf<block_index, gf<imfreq>> gw`` can be fitted using the two following functions:
+
+ ``void set_tail_from_fit(gf<imfreq> &gf, tail_view known_moments, int n_moments, size_t wn_min, size_t wn_max, bool replace_by_fit = false);``
+
+ ``void set_tail_from_fit(gf<block_index, gf<imfreq>> &block_gf, tail_view known_moments, int n_moments, size_t wn_min,  size_t wn_max, bool replace_by_fit = false);``
+
+
+where
+
++-------------+----------------+----------------------------------------------------------------------+----------+
+|  type       | name           | description                                                          | default  |
++=============+================+======================================================================+==========+
+|  gf<imfreq> | gf             | Green's function to be fit                                           | no       |
++-------------+----------------+----------------------------------------------------------------------+----------+
+|  tail_view  | known_moments  | known part of the tail                                               | no       |
++-------------+----------------+----------------------------------------------------------------------+----------+
+|  int        | n_moments      | number of moments in the final tail (including known ones)           | no       |
++-------------+----------------+----------------------------------------------------------------------+----------+
+|  size_t     | wn_min         | frequency to start the fit                                           | no       |
++-------------+----------------+----------------------------------------------------------------------+----------+
+|  size_t     | wn_max         | final fitting frequency (included)                                   | no       |
++-------------+----------------+----------------------------------------------------------------------+----------+
+|  bool       | replace_by_fit | if true, replace the gf data in the fitting range by the tail values | true     |
++-------------+----------------+----------------------------------------------------------------------+----------+
+
+
+Example
+~~~~~~~~
+
+.. compileblock:: 
+
+  #include <triqs/gfs.hpp>
+  #include <triqs/gfs/local/fit_tail.hpp>
+  using namespace triqs::gfs; 
+  int main(){
+   triqs::clef::placeholder<0> iom_;
+   double beta =10;
+   int N=100;
+
+   auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}};
+   gw(iom_) << 1/(iom_-1);
+
+   size_t wn_min=50, wn_max=90;
+   int n_moments=4;  
+   int  size=1; //means that we know one moment
+   int order_min=1; //means that the first moment in the final tail will be the first moment
+   auto known_moments = local::tail(make_shape(1,1), size, order_min); //length is 0, first moment to fit is order_min
+   known_moments(1)=1.;//set the first moment
+
+   set_tail_from_fit(gw, known_moments, n_moments, wn_min, wn_max, true);
+   
+   std::cout << gw.singularity() << std::endl;
+  }
+
+
+
+Implementation
+~~~~~~~~~~~~~~~
+
+The fitting problem is written as a linear system of equations, which is solved by SVD.
+

doc/reference/c++/gf/tail.rst

@@ -38,6 +38,41 @@ In TRIQS, the tail is implemented as an object ``tail``. Here is a simple exampl
    std::cout << t << std::endl;
   }
 
+
+
+Fitting the tail of a Green's function
+---------------------------------------
+
+Given an imaginary-frequency Green's function, one can compute the moments of its high-frequency tail with the function ``set_tail_from_fit``:
+
+
+.. compileblock:: 
+
+  #include <triqs/gfs.hpp>
+  #include <triqs/gfs/local/fit_tail.hpp>
+  using namespace triqs::gfs; 
+  int main(){
+   triqs::clef::placeholder<0> iom_;
+   double beta =10;
+   int N=100;
+
+   auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}};
+   gw(iom_) << 1/(iom_-1);
+
+   size_t wn_min=50; //frequency to start the fit
+   size_t wn_max=90; //final fitting frequency (included)
+   int n_moments=4;  //number of moments in the final tail (including known ones)
+   int  size=1; //means that we know one moment
+   int order_min=1; //means that the first moment in the final tail will be the first moment
+   auto known_moments = local::tail(make_shape(1,1), size, order_min); //length is 0, first moment to fit is order_min
+   known_moments(1)=1.;//set the first moment
+   set_tail_from_fit(gw, known_moments, n_moments, wn_min, wn_max, true);//true replace the gf data in the fitting range by the tail values
+   std::cout << gw.singularity() << std::endl;
+  }
+
+
+The full documentation of ``set_tail_from_fit`` is :doc:`here<set_tail_from_fit>`.
+
 API
 ****
 Here are the main methods of the ``tail`` class:

test/triqs/gfs/test_fit_tail.cpp

@@ -7,8 +7,7 @@ using triqs::arrays::make_shape;
 using namespace triqs::gfs;
 using triqs::gfs::local::tail;
 #define TEST(X) std::cout << BOOST_PP_STRINGIZE((X)) << " ---> "<< (X) <<std::endl<<std::endl;
-
+void test_0(){
-int main() {
 
  double precision=10e-9;
 
@@ -29,9 +28,7 @@ int main() {
 
  gw(iom_) << c(0)/iom_ + c(1)/iom_/iom_ + c(2)/iom_/iom_/iom_; 
 
- //show tail
+ TEST(gw.singularity());
-// std::cout<< "before fitting:" <<std::endl;
-// for(auto &i : gw.singularity().data()) std::cout << i << std::endl;
 
  //erase tail
  for(auto &i : gw.singularity().data()) i = 0.0;
@@ -43,13 +40,12 @@ int main() {
  //restore tail
  set_tail_from_fit(gw, known_moments, n_moments, wn_min, wn_max);
 
-// std::cout<< "after fitting:" <<std::endl; 
+ TEST(gw.singularity());
-// for(auto &i : gw.singularity().data()) std::cout << i << std::endl;
 
  for(size_t i=0; i<first_dim(c); i++){
-   double diff = std::abs( c(i) - gw.singularity().data()(i,0,0) );
+  double diff = std::abs( c(i) - gw.singularity().data()(i,0,0) );
-   //std::cout<< "diff: " << diff <<std::endl;
+  //std::cout<< "diff: " << diff <<std::endl;
-   if (diff > precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="<<diff<<"\n";
+  if (diff > precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="<<diff<<"\n";
  }
 
 
@@ -64,14 +60,37 @@ int main() {
  known_moments(1)=1.;//set the first moment
  set_tail_from_fit(gw, known_moments, n_moments, wn_min, wn_max, true);//true replace the gf data in the fitting range by the tail values
 
+ TEST(gw.singularity());
 
  for(size_t i=0; i<first_dim(c); i++){
-   double diff = std::abs( c(i) - gw.singularity().data()(i,0,0) );
+  double diff = std::abs( c(i) - gw.singularity().data()(i,0,0) );
-   //std::cout<< "diff: " << diff <<std::endl;
+  //std::cout<< "diff: " << diff <<std::endl;
-   if (diff > precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="<<diff<<"\n";
+  if (diff > precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="<<diff<<"\n";
  }
+}
+void test_1(){
+ //real life test: find tails of 1/(iom -1)
+ triqs::clef::placeholder<0> iom_;
+ double beta =10;
+ int N=100;
 
+ auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}};
+ gw(iom_) << 1/(iom_-1);
 
+ size_t wn_min=50; //frequency to start the fit
+ size_t wn_max=90; //final fitting frequency (included)
+ int n_moments=4;  //number of moments in the final tail (including known ones)
+ int  size=1; //means that we know one moment
+ int order_min=1; //means that the first moment in the final tail will be the first moment
+ auto known_moments = tail(make_shape(1,1), size, order_min); //length is 0, first moment to fit is order_min
+ known_moments(1)=1.;//set the first moment
+ set_tail_from_fit(gw, known_moments, n_moments, wn_min, wn_max, true);//true replace the gf data in the fitting range by the tail values
+ TEST(gw.singularity());
+}
+
+int main() {
+ test_0();
+ test_1();
 }
 
 

test/triqs/gfs/test_fit_tail.output

@@ -0,0 +1,48 @@
+(gw.singularity()) ---> tail/tail_view: min/smallest/max = -1 1 8
+ ...  Order -1 = 
+[[(0,0)]]
+ ...  Order 0 = 
+[[(0,0)]]
+ ...  Order 1 = 
+[[(1,0)]]
+ ...  Order 2 = 
+[[(3,0)]]
+ ...  Order 3 = 
+[[(5,0)]]
+ ...  Order 4 = 
+[[(0,0)]]
+ ...  Order 5 = 
+[[(0,0)]]
+ ...  Order 6 = 
+[[(0,0)]]
+ ...  Order 7 = 
+[[(0,0)]]
+ ...  Order 8 = 
+[[(0,0)]]
+
+(gw.singularity()) ---> tail/tail_view: min/smallest/max = 1 1 3
+ ...  Order 1 = 
+[[(1,0)]]
+ ...  Order 2 = 
+[[(3,0)]]
+ ...  Order 3 = 
+[[(5,0)]]
+
+(gw.singularity()) ---> tail/tail_view: min/smallest/max = 1 1 3
+ ...  Order 1 = 
+[[(1,0)]]
+ ...  Order 2 = 
+[[(3,0)]]
+ ...  Order 3 = 
+[[(5,0)]]
+
+(gw.singularity()) ---> tail/tail_view: min/smallest/max = 1 1 4
+ ...  Order 1 = 
+[[(1,0)]]
+ ...  Order 2 = 
+[[(1,0)]]
+ ...  Order 3 = 
+[[(0.999251,0)]]
+ ...  Order 4 = 
+[[(0.998655,0)]]
+