diff --git a/c++/triqs/gfs/meshes/tail_fitter.hpp b/c++/triqs/gfs/meshes/tail_fitter.hpp index 462208168e..ec75ed84d3 100644 --- a/c++/triqs/gfs/meshes/tail_fitter.hpp +++ b/c++/triqs/gfs/meshes/tail_fitter.hpp @@ -30,10 +30,10 @@ namespace triqs::gfs { //---------------------------------------------------------------------------------------------- // construct the Vandermonde matrix - inline arrays::matrix vander(std::vector const &pts, int expansion_order) { + inline arrays::matrix vander(std::vector const &pts, int expansion_order, int error_order) { arrays::matrix 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,6 +72,7 @@ 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>, max_order + 1> _lss; std::array, max_order + 1> _lss_hermitian; @@ -79,9 +80,10 @@ namespace triqs::gfs { std::vector _fit_idx_lst; public: - tail_fitter(double tail_fraction, int n_tail_max, std::optional expansion_order = {}) + tail_fitter(double tail_fraction, int n_tail_max, std::optional expansion_order = {}, std::optional 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 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 expansion_order = {}) const { - _tail_fitter = std::make_shared(tail_fitter{tail_fraction, n_tail_max, expansion_order}); + std::optional expansion_order = {}, std::optional error_order = {}) const { + _tail_fitter = std::make_shared(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 expansion_order = {}) const { - set_tail_fit_parameters(tail_fraction, n_tail_max, expansion_order); + std::optional expansion_order = {}, std::optional error_order = {}) const { + set_tail_fit_parameters(tail_fraction, n_tail_max, expansion_order, error_order); return *_tail_fitter; } diff --git a/doc/documentation/manual/triqs/gfs/c++/tail.rst b/doc/documentation/manual/triqs/gfs/c++/tail.rst index ad23b83826..16de2b4c23 100644 --- a/doc/documentation/manual/triqs/gfs/c++/tail.rst +++ b/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) ----------------------------------- diff --git a/python/triqs/gf/meshes_desc.py b/python/triqs/gf/meshes_desc.py index 38b3d942e3..e49c2acf74 100644 --- a/python/triqs/gf/meshes_desc.py +++ b/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 expansion_order = {})""") +m.add_method("""void set_tail_fit_parameters(double tail_fraction, int n_tail_max = 30, std::optional expansion_order = {}, std::optional error_order = {})""") m.add_property(name = "beta", getter = cfunction(calling_pattern="double result = self_c.domain().beta", signature = "double()",