From da7e7ec971dc18acd43fba162c97deb44b50f507 Mon Sep 17 00:00:00 2001 From: tayral Date: Tue, 4 Feb 2014 16:13:14 +0000 Subject: [PATCH] Fixed fit_tail for pos. and neg. matsub + bosonic -> code was previously assuming mesh with only positive, fermionic matsubara freqs -> changed wn_min to n_min (was misleading, since it was an index, not a frequency) / same for <-> max -> changed doc accordingly --- doc/reference/c++/gf/set_tail_from_fit.rst | 38 ++++++++++---------- doc/reference/c++/gf/tail.rst | 6 ++-- test/triqs/gfs/test_fit_tail.cpp | 40 ++++++++++++++++++++++ test/triqs/gfs/test_fit_tail.output | 20 +++++++++++ triqs/gfs/local/fit_tail.hpp | 36 +++++++++---------- 5 files changed, 100 insertions(+), 40 deletions(-) diff --git a/doc/reference/c++/gf/set_tail_from_fit.rst b/doc/reference/c++/gf/set_tail_from_fit.rst index dbd27796..e10826dc 100644 --- a/doc/reference/c++/gf/set_tail_from_fit.rst +++ b/doc/reference/c++/gf/set_tail_from_fit.rst @@ -6,28 +6,28 @@ API The tail of a given ``gf/gf> gw`` can be fitted using the two following functions: - ``void set_tail_from_fit(gf &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 &gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max, bool replace_by_fit = false);`` - ``void set_tail_from_fit(gf> &block_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_gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max, bool replace_by_fit = false);`` where -+-------------+----------------+----------------------------------------------------------------------+----------+ -| type | name | description | default | -+=============+================+======================================================================+==========+ -| gf | 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 | -+-------------+----------------+----------------------------------------------------------------------+----------+ ++------------+----------------+----------------------------------------------------------------------+---------+ +| type | name | description | default | ++============+================+======================================================================+=========+ +| gf | gf | Green's function to be fit (bosonic/fermionic) | 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 | n_min | linear index on mesh to start the fit | no | ++------------+----------------+----------------------------------------------------------------------+---------+ +| size_t | n_max | final linear index for fit (included) | no | ++------------+----------------+----------------------------------------------------------------------+---------+ +| bool | replace_by_fit | if true, replace the gf data in the fitting range by the tail values | true | ++------------+----------------+----------------------------------------------------------------------+---------+ Example @@ -46,14 +46,14 @@ Example auto gw = gf{{beta, Fermion, N}, {1, 1}}; gw(iom_) << 1/(iom_-1); - size_t wn_min=50, wn_max=90; + size_t n_min=50, n_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); + set_tail_from_fit(gw, known_moments, n_moments, n_min, n_max, true); std::cout << gw.singularity() << std::endl; } diff --git a/doc/reference/c++/gf/tail.rst b/doc/reference/c++/gf/tail.rst index 71d4aef7..405737a6 100644 --- a/doc/reference/c++/gf/tail.rst +++ b/doc/reference/c++/gf/tail.rst @@ -59,14 +59,14 @@ Given an imaginary-frequency Green's function, one can compute the moments of it auto gw = gf{{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) + size_t n_min=50; //linear index on mesh to start the fit + size_t n_max=90; //final linear index for fit (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 + set_tail_from_fit(gw, known_moments, n_moments, n_min, n_max, true);//true replace the gf data in the fitting range by the tail values std::cout << gw.singularity() << std::endl; } diff --git a/test/triqs/gfs/test_fit_tail.cpp b/test/triqs/gfs/test_fit_tail.cpp index 65dc8001..3cf4b6b8 100644 --- a/test/triqs/gfs/test_fit_tail.cpp +++ b/test/triqs/gfs/test_fit_tail.cpp @@ -87,10 +87,50 @@ void test_1(){ 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()); } +void test_2(){ + //real life test: find tails of 1/(iom -1) -- with positive and negative matsubara + triqs::clef::placeholder<0> iom_; + double beta =10; + int N=200; + + auto gw = gf{{beta, Fermion, N, false}, {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()); +} +void test_3(){ + //real life test: find tails of 1/(iom -1) --> bosonic case + triqs::clef::placeholder<0> iom_; + double beta =10; + int N=100; + + auto gw = gf{{beta, Boson, 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_2(); + test_3(); } diff --git a/test/triqs/gfs/test_fit_tail.output b/test/triqs/gfs/test_fit_tail.output index 53ac69ae..371d7a66 100644 --- a/test/triqs/gfs/test_fit_tail.output +++ b/test/triqs/gfs/test_fit_tail.output @@ -46,3 +46,23 @@ ... Order 4 = [[(0.998655,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)]] + +(gw.singularity()) ---> tail/tail_view: min/smallest/max = 1 1 4 + ... Order 1 = +[[(1,0)]] + ... Order 2 = +[[(1,0)]] + ... Order 3 = +[[(0.999236,0)]] + ... Order 4 = +[[(0.998631,0)]] + diff --git a/triqs/gfs/local/fit_tail.hpp b/triqs/gfs/local/fit_tail.hpp index fc98151d..bcdff120 100644 --- a/triqs/gfs/local/fit_tail.hpp +++ b/triqs/gfs/local/fit_tail.hpp @@ -46,11 +46,11 @@ namespace gfs { // the input gf Green's function: gf // the known moments in the form of a tail(_view): known_moments // the TOTAL number of desired moments (including the known ones): n_moments - // the index of the first and last frequency to fit (the last one is included): wn_min, wn_max + // the index of the first and last frequency to fit (the last one is included): n_min, n_max // output: returns the tail obtained by fitting - tail fit_tail_impl(gf &gf, const tail_view known_moments, int n_moments, int wn_min, int wn_max) { + tail fit_tail_impl(gf &gf, const tail_view known_moments, int n_moments, int n_min, int n_max) { tail res(get_target_shape(gf), n_moments, known_moments.order_min()); if (known_moments.size()) @@ -73,7 +73,7 @@ namespace gfs { if (n_unknown_moments % 2 != 0 && omin % 2 == 0) size_even += 1; int size_odd = n_unknown_moments - size_even; - int size1 = wn_max - wn_min + 1; + int size1 = n_max - n_min + 1; // size2 is the number of moments arrays::matrix A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT); @@ -89,10 +89,10 @@ namespace gfs { // S.resize(size_odd); // A.resize(size1,size_odd); //when resizing, gelss segfaults for (int k = 0; k < size1; k++) { - auto n = wn_min + k; - auto iw = std::complex(0.0, (2 * n + 1) * M_PI / beta); + auto n = n_min + k; + auto iw = std::complex(gf.mesh().index_to_point(n)); - B(k, 0) = imag(gf.data()(n, i, j)); + B(k, 0) = imag(gf.data()(gf.mesh().index_to_linear(n), i, j)); // subtract known tail if present if (known_moments.size() > 0) B(k, 0) -= imag(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)).evaluate(iw)(0, 0)); @@ -112,10 +112,10 @@ namespace gfs { // S.resize(size_even); // A.resize(size1,size_even); //when resizing, gelss segfaults for (int k = 0; k < size1; k++) { - auto n = wn_min + k; - auto iw = std::complex(0.0, (2 * n + 1) * M_PI / beta); - - B(k, 0) = real(gf.data()(n, i, j)); + auto n = n_min + k; + auto iw = std::complex(gf.mesh().index_to_point(n)); + + B(k, 0) = real(gf.data()(gf.mesh().index_to_linear(n), i, j)); // subtract known tail if present if (known_moments.size() > 0) B(k, 0) -= real(slice_target(known_moments, arrays::range(i, i + 1), arrays::range(j, j + 1)).evaluate(iw)(0, 0)); @@ -137,26 +137,26 @@ namespace gfs { } - void set_tail_from_fit(gf &gf, tail_view known_moments, int n_moments, size_t wn_min, size_t wn_max, + void set_tail_from_fit(gf &gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max, bool replace_by_fit = false) { if (get_target_shape(gf) != known_moments.shape()) TRIQS_RUNTIME_ERROR << "shape of tail does not match shape of gf"; - gf.singularity() = fit_tail_impl(gf, known_moments, n_moments, wn_min, wn_max); + gf.singularity() = fit_tail_impl(gf, known_moments, n_moments, n_min, n_max); if (replace_by_fit) { // replace data in the fitting range by the values from the fitted tail size_t i = 0; - for (auto iw : gf.mesh()) { // (arrays::range(wn_min,wn_max+1)) { - if ((i >= wn_min) && (i <= wn_max)) gf[iw] = gf.singularity().evaluate(iw); + for (auto iw : gf.mesh()) { // (arrays::range(n_min,n_max+1)) { + if ((i >= n_min) && (i <= n_max)) gf[iw] = gf.singularity().evaluate(iw); i++; } } } - void set_tail_from_fit(gf> &block_gf, tail_view known_moments, int n_moments, size_t wn_min, - size_t wn_max, bool replace_by_fit = false) { - // for(auto &gf : block_gf) set_tail_from_fit(gf, known_moments, n_moments, wn_min, wn_max, replace_by_fit); + void set_tail_from_fit(gf> &block_gf, tail_view known_moments, int n_moments, size_t n_min, + size_t n_max, bool replace_by_fit = false) { + // for(auto &gf : block_gf) set_tail_from_fit(gf, known_moments, n_moments, n_min, n_max, replace_by_fit); for (size_t i = 0; i < block_gf.mesh().size(); i++) - set_tail_from_fit(block_gf[i], known_moments, n_moments, wn_min, wn_max, replace_by_fit); + set_tail_from_fit(block_gf[i], known_moments, n_moments, n_min, n_max, replace_by_fit); } } }