LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-11-01 03:33:50 +01:00
Code Issues Releases Wiki Activity

Fixed fit_tail for pos. and neg. matsub + bosonic

Browse Source 
-> 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
This commit is contained in:
tayral 2014-02-04 16:13:14 +00:00 committed by Olivier Parcollet
parent b129b3d17b
commit da7e7ec971
5 changed files with 100 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
doc/reference/c++/gf/set_tail_from_fit.rst
View File

@ -6,28 +6,28 @@ API
The tail of a given ``gf<imfreq>/gf<block_index, gf<imfreq>> gw`` can be fitted using the two following functions: 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<imfreq> &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_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);`` ``void set_tail_from_fit(gf<block_index, gf<imfreq>> &block_gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max, bool replace_by_fit = false);``
where where
+-------------+----------------+----------------------------------------------------------------------+----------+ +------------+----------------+----------------------------------------------------------------------+---------+
| type | name | description | default | | type | name | description | default |
+=============+================+======================================================================+==========+ +============+================+======================================================================+=========+
| gf<imfreq> | gf | Green's function to be fit | no | | gf<imfreq> | gf | Green's function to be fit (bosonic/fermionic) | no |
+-------------+----------------+----------------------------------------------------------------------+----------+ +------------+----------------+----------------------------------------------------------------------+---------+
| tail_view | known_moments | known part of the tail | 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 | | 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 | n_min | linear index on mesh to start the fit | no |
+-------------+----------------+----------------------------------------------------------------------+----------+ +------------+----------------+----------------------------------------------------------------------+---------+
| size_t | wn_max | final fitting frequency (included) | 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 | | bool | replace_by_fit | if true, replace the gf data in the fitting range by the tail values | true |
+-------------+----------------+----------------------------------------------------------------------+----------+ +------------+----------------+----------------------------------------------------------------------+---------+
Example Example
@ -46,14 +46,14 @@ Example
auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}}; auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}};
gw(iom_) << 1/(iom_-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 n_moments=4;
int size=1; //means that we know one moment 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 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 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 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; std::cout << gw.singularity() << std::endl;
} }

6
doc/reference/c++/gf/tail.rst
View File

@ -59,14 +59,14 @@ Given an imaginary-frequency Green's function, one can compute the moments of it
auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}}; auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}};
gw(iom_) << 1/(iom_-1); gw(iom_) << 1/(iom_-1);
size_t wn_min=50; //frequency to start the fit size_t n_min=50; //linear index on mesh to start the fit
size_t wn_max=90; //final fitting frequency (included) 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 n_moments=4; //number of moments in the final tail (including known ones)
int size=1; //means that we know one moment 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 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 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 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; std::cout << gw.singularity() << std::endl;
} }

40
test/triqs/gfs/test_fit_tail.cpp
View File

@ -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 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()); 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<imfreq>{{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<imfreq>{{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() { int main() {
test_0(); test_0();
test_1(); test_1();
test_2();
test_3();
} }

20
test/triqs/gfs/test_fit_tail.output
View File

@ -46,3 +46,23 @@
... Order 4 = ... Order 4 =
[[(0.998655,0)]] [[(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)]]

34
triqs/gfs/local/fit_tail.hpp
View File

@ -46,11 +46,11 @@ namespace gfs {
// the input gf<imfreq> Green's function: gf // the input gf<imfreq> Green's function: gf
// the known moments in the form of a tail(_view): known_moments // 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 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 // output: returns the tail obtained by fitting
tail fit_tail_impl(gf<imfreq> &gf, const tail_view known_moments, int n_moments, int wn_min, int wn_max) { tail fit_tail_impl(gf<imfreq> &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()); tail res(get_target_shape(gf), n_moments, known_moments.order_min());
if (known_moments.size()) if (known_moments.size())
@ -73,7 +73,7 @@ namespace gfs {
if (n_unknown_moments % 2 != 0 && omin % 2 == 0) size_even += 1; if (n_unknown_moments % 2 != 0 && omin % 2 == 0) size_even += 1;
int size_odd = n_unknown_moments - size_even; 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 // size2 is the number of moments
arrays::matrix<double, 2> A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT); arrays::matrix<double, 2> A(size1, std::max(size_even, size_odd), FORTRAN_LAYOUT);
@ -89,10 +89,10 @@ namespace gfs {
// S.resize(size_odd); // S.resize(size_odd);
// A.resize(size1,size_odd); //when resizing, gelss segfaults // A.resize(size1,size_odd); //when resizing, gelss segfaults
for (int k = 0; k < size1; k++) { for (int k = 0; k < size1; k++) {
auto n = wn_min + k; auto n = n_min + k;
auto iw = std::complex<double>(0.0, (2 * n + 1) * M_PI / beta); auto iw = std::complex<double>(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 // subtract known tail if present
if (known_moments.size() > 0) 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(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); // S.resize(size_even);
// A.resize(size1,size_even); //when resizing, gelss segfaults // A.resize(size1,size_even); //when resizing, gelss segfaults
for (int k = 0; k < size1; k++) { for (int k = 0; k < size1; k++) {
auto n = wn_min + k; auto n = n_min + k;
auto iw = std::complex<double>(0.0, (2 * n + 1) * M_PI / beta); auto iw = std::complex<double>(gf.mesh().index_to_point(n));
B(k, 0) = real(gf.data()(n, i, j)); B(k, 0) = real(gf.data()(gf.mesh().index_to_linear(n), i, j));
// subtract known tail if present // subtract known tail if present
if (known_moments.size() > 0) 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)); 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<imfreq> &gf, tail_view known_moments, int n_moments, size_t wn_min, size_t wn_max, void set_tail_from_fit(gf<imfreq> &gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max,
bool replace_by_fit = false) { 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"; 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 if (replace_by_fit) { // replace data in the fitting range by the values from the fitted tail
size_t i = 0; size_t i = 0;
for (auto iw : gf.mesh()) { // (arrays::range(wn_min,wn_max+1)) { for (auto iw : gf.mesh()) { // (arrays::range(n_min,n_max+1)) {
if ((i >= wn_min) && (i <= wn_max)) gf[iw] = gf.singularity().evaluate(iw); if ((i >= n_min) && (i <= n_max)) gf[iw] = gf.singularity().evaluate(iw);
i++; i++;
} }
} }
} }
void set_tail_from_fit(gf<block_index, gf<imfreq>> &block_gf, tail_view known_moments, int n_moments, size_t wn_min, void set_tail_from_fit(gf<block_index, gf<imfreq>> &block_gf, tail_view known_moments, int n_moments, size_t n_min,
size_t wn_max, bool replace_by_fit = false) { size_t n_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); // 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++) 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);
} }
} }
} }