diff --git a/test/triqs/gfs/test_fit_tail.cpp b/test/triqs/gfs/test_fit_tail.cpp index 3cf4b6b8..57a55acb 100644 --- a/test/triqs/gfs/test_fit_tail.cpp +++ b/test/triqs/gfs/test_fit_tail.cpp @@ -15,7 +15,8 @@ void test_0(){ double beta =10; int N=100; - auto gw = gf{{beta, Fermion, N}, {1, 1}}; + auto gw = gf{{beta, Fermion, N},{1,1}}; + auto gw_s = gf{{beta, Fermion, N}}; triqs::arrays::array c(3); triqs::clef::placeholder<1> i_; @@ -27,11 +28,13 @@ void test_0(){ auto known_moments = tail(make_shape(1,1), size, order_min); //length is 0, first moment to fit is order_min gw(iom_) << c(0)/iom_ + c(1)/iom_/iom_ + c(2)/iom_/iom_/iom_; + gw_s(iom_) << c(0)/iom_ + c(1)/iom_/iom_ + c(2)/iom_/iom_/iom_; TEST(gw.singularity()); //erase tail for(auto &i : gw.singularity().data()) i = 0.0; + for(auto &i : gw_s.singularity().data()) i = 0.0; size_t wn_min=50; //frequency to start the fit size_t wn_max=90; //final fitting frequency (included) @@ -39,16 +42,17 @@ void test_0(){ //restore tail set_tail_from_fit(gw, known_moments, n_moments, wn_min, wn_max); + set_tail_from_fit(gw_s, known_moments, n_moments, wn_min, wn_max); TEST(gw.singularity()); - + TEST(gw_s.singularity()); +/* for(size_t i=0; i precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="< precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="< iom_; @@ -75,7 +81,9 @@ void test_1(){ int N=100; auto gw = gf{{beta, Fermion, N}, {1, 1}}; + auto gw_b = gf{{beta, Boson, N}, {1, 1}}; gw(iom_) << 1/(iom_-1); + gw_b(iom_) << 1/(iom_-1); size_t wn_min=50; //frequency to start the fit size_t wn_max=90; //final fitting frequency (included) @@ -85,8 +93,11 @@ void test_1(){ 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 + set_tail_from_fit(gw_b, 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()); + TEST(gw_b.singularity()); } + void test_2(){ //real life test: find tails of 1/(iom -1) -- with positive and negative matsubara triqs::clef::placeholder<0> iom_; @@ -106,31 +117,12 @@ void test_2(){ 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 371d7a66..3d5225f7 100644 --- a/test/triqs/gfs/test_fit_tail.output +++ b/test/triqs/gfs/test_fit_tail.output @@ -20,7 +20,11 @@ ... Order 8 = [[(0,0)]] -(gw.singularity()) ---> tail/tail_view: min/smallest/max = 1 1 3 +(gw.singularity()) ---> tail/tail_view: min/smallest/max = -1 1 3 + ... Order -1 = +[[(0,0)]] + ... Order 0 = +[[(0,0)]] ... Order 1 = [[(1,0)]] ... Order 2 = @@ -28,7 +32,11 @@ ... Order 3 = [[(5,0)]] -(gw.singularity()) ---> tail/tail_view: min/smallest/max = 1 1 3 +(gw_s.singularity()) ---> tail/tail_view: min/smallest/max = -1 1 3 + ... Order -1 = +[[(0,0)]] + ... Order 0 = +[[(0,0)]] ... Order 1 = [[(1,0)]] ... Order 2 = @@ -36,7 +44,23 @@ ... Order 3 = [[(5,0)]] -(gw.singularity()) ---> tail/tail_view: min/smallest/max = 1 1 4 +(gw.singularity()) ---> tail/tail_view: min/smallest/max = -1 1 3 + ... Order -1 = +[[(0,0)]] + ... Order 0 = +[[(0,0)]] + ... 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 = +[[(0,0)]] + ... Order 0 = +[[(0,0)]] ... Order 1 = [[(1,0)]] ... Order 2 = @@ -46,23 +70,31 @@ ... Order 4 = [[(0.998655,0)]] -(gw.singularity()) ---> tail/tail_view: min/smallest/max = 1 1 4 +(gw_b.singularity()) ---> tail/tail_view: min/smallest/max = -1 1 4 + ... Order -1 = +[[(0,0)]] + ... Order 0 = +[[(0,0)]] ... 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)]] +[[(0.999209,0)]] ... Order 4 = [[(0.998631,0)]] +(gw.singularity()) ---> tail/tail_view: min/smallest/max = -1 1 4 + ... Order -1 = +[[(0,0)]] + ... Order 0 = +[[(0,0)]] + ... Order 1 = +[[(1,0)]] + ... Order 2 = +[[(1,0)]] + ... Order 3 = +[[(0.999251,0)]] + ... Order 4 = +[[(0.998655,0)]] + diff --git a/triqs/gfs/local/fit_tail.hpp b/triqs/gfs/local/fit_tail.hpp index bcdff120..cb9d704a 100644 --- a/triqs/gfs/local/fit_tail.hpp +++ b/triqs/gfs/local/fit_tail.hpp @@ -26,139 +26,136 @@ #include #include -namespace triqs { -namespace gfs { - namespace local { +namespace triqs { namespace gfs { namespace local { - using triqs::gfs::imfreq; - using triqs::gfs::block_index; - using triqs::gfs::block_index; + using triqs::gfs::imfreq; + using triqs::gfs::block_index; + using triqs::gfs::block_index; - namespace tgl = triqs::gfs::local; + namespace tgl = triqs::gfs::local; - // routine for fitting the tail (singularity) of a Matsubara Green's function - // this is a *linear* least squares problem (with non-linear basis functions) - // which is solved by singular value decomposition of the design matrix - // the routine fits the real part (even moments) and the imaginary part - //(odd moments) separately, since this is more stable + // routine for fitting the tail (singularity) of a Matsubara Green's function + // this is a *linear* least squares problem (with non-linear basis functions) + // which is solved by singular value decomposition of the design matrix + // the routine fits the real part (even moments) and the imaginary part + //(odd moments) separately, since this is more stable - // input: - // 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): n_min, n_max + // input: + // 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): n_min, n_max - // output: returns the tail obtained by fitting + // output: returns the tail obtained by fitting - tail fit_tail_impl(gf &gf, const tail_view known_moments, int n_moments, int n_min, int n_max) { + tail fit_tail_impl(gf_view 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()) - for (int i = known_moments.order_min(); i <= known_moments.order_max(); i++) res(i) = known_moments(i); + tail res(get_target_shape(gf)); + if (known_moments.size()) + for (int i = known_moments.order_min(); i <= known_moments.order_max(); i++) res(i) = known_moments(i); - // if known_moments.size()==0, the lowest order to be obtained from the fit is determined by order_min in known_moments - // if known_moments.size()==0, the lowest order is the one following order_max in known_moments + // if known_moments.size()==0, the lowest order to be obtained from the fit is determined by order_min in known_moments + // if known_moments.size()==0, the lowest order is the one following order_max in known_moments - const double beta = gf.mesh().domain().beta; + const double beta = gf.mesh().domain().beta; - int n_unknown_moments = n_moments - known_moments.size(); - if (n_unknown_moments < 1) return known_moments; + int n_unknown_moments = n_moments - known_moments.size(); + if (n_unknown_moments < 1) return known_moments; - // get the number of even unknown moments: it is n_unknown_moments/2+1 if the first - // moment is even and n_moments is odd; n_unknown_moments/2 otherwise - int omin = known_moments.size() == 0 ? known_moments.order_min() : known_moments.order_max() + 1; // smallest unknown moment - int omin_even = omin % 2 == 0 ? omin : omin + 1; - int omin_odd = omin % 2 != 0 ? omin : omin + 1; - int size_even = n_unknown_moments / 2; - if (n_unknown_moments % 2 != 0 && omin % 2 == 0) size_even += 1; - int size_odd = n_unknown_moments - size_even; + // get the number of even unknown moments: it is n_unknown_moments/2+1 if the first + // moment is even and n_moments is odd; n_unknown_moments/2 otherwise + int omin = known_moments.size() == 0 ? known_moments.order_min() : known_moments.order_max() + 1; // smallest unknown moment + int omin_even = omin % 2 == 0 ? omin : omin + 1; + int omin_odd = omin % 2 != 0 ? omin : omin + 1; + int size_even = n_unknown_moments / 2; + if (n_unknown_moments % 2 != 0 && omin % 2 == 0) size_even += 1; + int size_odd = n_unknown_moments - size_even; - int size1 = n_max - n_min + 1; - // size2 is the number of moments + 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); - arrays::matrix B(size1, 1, FORTRAN_LAYOUT); - arrays::vector S(std::max(size_even, size_odd)); - const double rcond = 0.0; - int rank; + arrays::matrix A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT); + arrays::matrix B(size1, 1, FORTRAN_LAYOUT); + arrays::vector S(std::max(size_even, size_odd)); + const double rcond = 0.0; + int rank; - for (size_t i = 0; i < get_target_shape(gf)[0]; i++) { - for (size_t j = 0; j < get_target_shape(gf)[1]; j++) { + for (size_t i = 0; i < get_target_shape(gf)[0]; i++) { + for (size_t j = 0; j < get_target_shape(gf)[1]; j++) { - // fit the odd moments - // S.resize(size_odd); - // A.resize(size1,size_odd); //when resizing, gelss segfaults - for (int k = 0; k < size1; k++) { - auto n = n_min + k; - auto iw = std::complex(gf.mesh().index_to_point(n)); + // fit the odd moments + // S.resize(size_odd); + // A.resize(size1,size_odd); //when resizing, gelss segfaults + for (int k = 0; k < size1; k++) { + auto n = n_min + k; + auto iw = std::complex(gf.mesh().index_to_point(n)); - 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)); + 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)); - for (int l = 0; l < size_odd; l++) { - int order = omin_odd + 2 * l; - A(k, l) = imag(pow(iw, -1.0 * order)); // set design matrix for odd moments - } - } - - arrays::lapack::gelss(A, B, S, rcond, rank); - for (int m = 0; m < size_odd; m++) { - res(omin_odd + 2 * m)(i, j) = B(m, 0); - } - - // fit the even moments - // S.resize(size_even); - // A.resize(size1,size_even); //when resizing, gelss segfaults - for (int k = 0; k < size1; k++) { - 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)); - - for (int l = 0; l < size_even; l++) { - int order = omin_even + 2 * l; - A(k, l) = real(pow(iw, -1.0 * order)); // set design matrix for odd moments - } - } - - arrays::lapack::gelss(A, B, S, rcond, rank); - for (int m = 0; m < size_even; m++) { - res(omin_even + 2 * m)(i, j) = B(m, 0); + for (int l = 0; l < size_odd; l++) { + int order = omin_odd + 2 * l; + A(k, l) = imag(pow(iw, -1.0 * order)); // set design matrix for odd moments } } - } - return res; // return tail - } + arrays::lapack::gelss(A, B, S, rcond, rank); + for (int m = 0; m < size_odd; m++) { + res(omin_odd + 2 * m)(i, j) = B(m, 0); + } + // fit the even moments + // S.resize(size_even); + // A.resize(size1,size_even); //when resizing, gelss segfaults + for (int k = 0; k < size1; k++) { + auto n = n_min + k; + auto iw = std::complex(gf.mesh().index_to_point(n)); - 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, n_min, n_max); + 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)); - 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(n_min,n_max+1)) { - if ((i >= n_min) && (i <= n_max)) gf[iw] = gf.singularity().evaluate(iw); - i++; + for (int l = 0; l < size_even; l++) { + int order = omin_even + 2 * l; + A(k, l) = real(pow(iw, -1.0 * order)); // set design matrix for odd moments + } + } + + arrays::lapack::gelss(A, B, S, rcond, rank); + for (int m = 0; m < size_even; m++) { + res(omin_even + 2 * m)(i, j) = B(m, 0); } } } + res.mask_view()=n_moments; + return res; // return tail + } + void set_tail_from_fit(gf_view 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, 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(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 n_min, - size_t n_max, bool replace_by_fit = false) { + void set_tail_from_fit(gf_view> 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, n_min, n_max, replace_by_fit); } + void set_tail_from_fit(gf_view gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max, bool replace_by_fit = false) { + set_tail_from_fit(reinterpret_scalar_valued_gf_as_matrix_valued(gf), known_moments, n_moments, n_min, n_max, replace_by_fit ); } -} -} // namespace + +}}} // namespace #endif