Weighted Tailfit

c++/triqs/gfs/meshes/tail_fitter.hpp

@@ -30,10 +30,10 @@ namespace triqs::gfs {
 
   //----------------------------------------------------------------------------------------------
   // construct the Vandermonde matrix
-  inline arrays::matrix<dcomplex> vander(std::vector<dcomplex> const &pts, int expansion_order) {
+  inline arrays::matrix<dcomplex> vander(std::vector<dcomplex> const &pts, int expansion_order, int error_order) {
     arrays::matrix<dcomplex> V(pts.size(), expansion_order + 1);
     for (auto [i, p] : itertools::enumerate(pts)) {
-      dcomplex z = 1;
+      dcomplex z = std::abs(std::pow( p , error_order ));
       for (int n = 0; n <= expansion_order; ++n) {
         V(i, n) = z;
         z *= p;
@@ -72,16 +72,18 @@ namespace triqs::gfs {
     const int _n_tail_max;
     const bool _adjust_order;
     const int _expansion_order;
+    const int _error_order;
     const double _rcond = 1e-8;
     std::array<std::unique_ptr<const arrays::lapack::gelss_cache<dcomplex>>, max_order + 1> _lss;
     std::array<std::unique_ptr<const arrays::lapack::gelss_cache_hermitian>, max_order + 1> _lss_hermitian;
     arrays::matrix<dcomplex> _vander;
     std::vector<long> _fit_idx_lst;
 
     public:
-    tail_fitter(double tail_fraction, int n_tail_max, std::optional<int> expansion_order = {})
+    tail_fitter(double tail_fraction, int n_tail_max, std::optional<int> expansion_order = {}, std::optional<int> error_order = {})
        : _tail_fraction(tail_fraction),
          _n_tail_max(n_tail_max),
+         _error_order(error_order.has_value() ? *error_order : 0),
          _adjust_order(not expansion_order.has_value()),
          _expansion_order(_adjust_order ? max_order : *expansion_order) {}
     //----------------------------------------------------------------------------------------------
@@ -152,7 +154,7 @@ namespace triqs::gfs {
         std::vector<dcomplex> C;
         C.reserve(_fit_idx_lst.size());
         for (long n : _fit_idx_lst) C.push_back(om_max / m.index_to_point(n));
-        _vander = vander(C, _expansion_order);
+        _vander = vander(C, _expansion_order, _error_order);
       }
 
       if (n_fixed_moments + 1 > first_dim(_vander) / 2) TRIQS_RUNTIME_ERROR << "Insufficient data points for least square procedure";
@@ -229,6 +231,8 @@ namespace triqs::gfs {
           g_mat(i, 0) = g_data_swap_idx(m.index_to_linear(n));
         else
           for (auto [j, x] : enumerate(g_data_swap_idx(m.index_to_linear(n), ellipsis()))) g_mat(i, j) = x;
+
+        for(int j=0; j < ncols; j++) g_mat(i, j) *= std::abs( std::pow( std::abs(m.omega_max()) / m.index_to_point(n), _error_order) );
       }
 
       // If an array with known_moments was passed, flatten the array into a matrix
@@ -304,8 +308,8 @@ namespace triqs::gfs {
 
     // Adjust the parameters for the tail-fitting
     void set_tail_fit_parameters(double tail_fraction, int n_tail_max = tail_fitter::default_n_tail_max,
-                                 std::optional<int> expansion_order = {}) const {
-      _tail_fitter = std::make_shared<tail_fitter>(tail_fitter{tail_fraction, n_tail_max, expansion_order});
+                                 std::optional<int> expansion_order = {}, std::optional<int> error_order = {}) const {
+      _tail_fitter = std::make_shared<tail_fitter>(tail_fitter{tail_fraction, n_tail_max, expansion_order, error_order});
     }
 
     // The tail fitter is mutable, even if the mesh is immutable to cache some data
@@ -316,8 +320,8 @@ namespace triqs::gfs {
 
     // Adjust the parameters for the tail-fitting and return the fitter
     tail_fitter &get_tail_fitter(double tail_fraction, int n_tail_max = tail_fitter::default_n_tail_max,
-                                 std::optional<int> expansion_order = {}) const {
-      set_tail_fit_parameters(tail_fraction, n_tail_max, expansion_order);
+                                 std::optional<int> expansion_order = {}, std::optional<int> error_order = {}) const {
+      set_tail_fit_parameters(tail_fraction, n_tail_max, expansion_order, error_order);
       return *_tail_fitter;
     }
 

doc/documentation/manual/triqs/gfs/c++/tail.rst

@@ -63,7 +63,7 @@ The experienced user can adjust the parameters of the fitting procedure directly
 
 .. code-block:: c
 
-   G.mesh().set_tail_fit_parameters(tail_fraction, n_tail_max, expansion_order)
+   G.mesh().set_tail_fit_parameters(tail_fraction, n_tail_max, expansion_order, error_order)
 
 The fitting parameters are
 
@@ -77,6 +77,9 @@ The fitting parameters are
 
 * :code:`expansion_order` (integer, optional)
   This parameter will fix the largest order to be considered for the fit.
+
+* :code:`error_order` (integer, optional, default = 0)
+  Assume an omega-dependent error of the Green's function for the fit. Error_order :math:`\eta` is defined such that :math:`\Delta G(i\omega_n) \sim |\omega_n|^\eta`
 
 Adjusting the fit window (Advanced)
 -----------------------------------

python/triqs/gf/meshes_desc.py

@@ -78,7 +78,7 @@ def make_mesh(py_type, c_tag, index_type='long'):
 m.add_method("""int last_index()""")
 m.add_method("""int first_index()""")
 m.add_method("""bool positive_only()""")
-m.add_method("""void set_tail_fit_parameters(double tail_fraction, int n_tail_max = 30, std::optional<int> expansion_order = {})""")
+m.add_method("""void set_tail_fit_parameters(double tail_fraction, int n_tail_max = 30, std::optional<int> expansion_order = {}, std::optional<int> error_order = {})""")
 m.add_property(name = "beta",
                getter = cfunction(calling_pattern="double result = self_c.domain().beta",
                signature = "double()",