LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-25 13:53:40 +01:00
Code Issues Releases Wiki Activity

Extending fit_tail for scalar_valued (+ test)

Browse Source 
-> extension by using reinterpret_scalar_valued_as_matrix_valued
This commit is contained in:
tayral 2014-02-04 16:58:38 +00:00 committed by Olivier Parcollet
parent da7e7ec971
commit 708c47305c
3 changed files with 162 additions and 141 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -15,7 +15,8 @@ void test_0(){
double beta =10; double beta =10;
int N=100; int N=100;
auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}}; auto gw = gf<imfreq>{{beta, Fermion, N},{1,1}};
auto gw_s = gf<imfreq, scalar_valued>{{beta, Fermion, N}};
triqs::arrays::array<double,1> c(3); triqs::arrays::array<double,1> c(3);
triqs::clef::placeholder<1> i_; 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 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(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()); TEST(gw.singularity());
//erase tail //erase tail
for(auto &i : gw.singularity().data()) i = 0.0; 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_min=50; //frequency to start the fit
size_t wn_max=90; //final fitting frequency (included) size_t wn_max=90; //final fitting frequency (included)
@ -39,16 +42,17 @@ void test_0(){
//restore tail //restore tail
set_tail_from_fit(gw, known_moments, n_moments, wn_min, wn_max); 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.singularity());
TEST(gw_s.singularity());
/*
for(size_t i=0; i<first_dim(c); i++){ for(size_t i=0; i<first_dim(c); i++){
double diff = std::abs( c(i) - gw.singularity().data()(i,0,0) ); double diff = std::abs( c(i) - gw.singularity().data()(i,0,0) );
//std::cout<< "diff: " << diff <<std::endl; //std::cout<< "diff: " << diff <<std::endl;
if (diff > precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="<<diff<<"\n"; if (diff > precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="<<diff<<"\n";
} }
*/
//erase tail //erase tail
for(auto &i : gw.singularity().data()) i = 0.0; for(auto &i : gw.singularity().data()) i = 0.0;
@ -61,13 +65,15 @@ void test_0(){
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());
/*
for(size_t i=0; i<first_dim(c); i++){ for(size_t i=0; i<first_dim(c); i++){
double diff = std::abs( c(i) - gw.singularity().data()(i,0,0) ); double diff = std::abs( c(i) - gw.singularity().data()(i,0,0) );
//std::cout<< "diff: " << diff <<std::endl; //std::cout<< "diff: " << diff <<std::endl;
if (diff > precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="<<diff<<"\n"; if (diff > precision) TRIQS_RUNTIME_ERROR<<" fit_tail error : diff="<<diff<<"\n";
} }
*/
} }
void test_1(){ void test_1(){
//real life test: find tails of 1/(iom -1) //real life test: find tails of 1/(iom -1)
triqs::clef::placeholder<0> iom_; triqs::clef::placeholder<0> iom_;
@ -75,7 +81,9 @@ void test_1(){
int N=100; int N=100;
auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}}; auto gw = gf<imfreq>{{beta, Fermion, N}, {1, 1}};
auto gw_b = gf<imfreq>{{beta, Boson, N}, {1, 1}};
gw(iom_) << 1/(iom_-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_min=50; //frequency to start the fit
size_t wn_max=90; //final fitting frequency (included) 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 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 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, 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.singularity());
TEST(gw_b.singularity());
} }
void test_2(){ void test_2(){
//real life test: find tails of 1/(iom -1) -- with positive and negative matsubara //real life test: find tails of 1/(iom -1) -- with positive and negative matsubara
triqs::clef::placeholder<0> iom_; 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 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_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_2();
test_3();
} }

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

@ -20,7 +20,11 @@
... Order 8 = ... Order 8 =
[[(0,0)]] [[(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 = ... Order 1 =
[[(1,0)]] [[(1,0)]]
... Order 2 = ... Order 2 =
@ -28,7 +32,11 @@
... Order 3 = ... Order 3 =
[[(5,0)]] [[(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 = ... Order 1 =
[[(1,0)]] [[(1,0)]]
... Order 2 = ... Order 2 =
@ -36,7 +44,23 @@
... Order 3 = ... Order 3 =
[[(5,0)]] [[(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 = ... Order 1 =
[[(1,0)]] [[(1,0)]]
... Order 2 = ... Order 2 =
@ -46,23 +70,31 @@
... Order 4 = ... Order 4 =
[[(0.998655,0)]] [[(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 = ... Order 1 =
[[(1,0)]] [[(1,0)]]
... Order 2 = ... Order 2 =
[[(1,0)]] [[(1,0)]]
... Order 3 = ... Order 3 =
[[(0.999251,0)]] [[(0.999209,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 = ... Order 4 =
[[(0.998631,0)]] [[(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)]]

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

@ -26,139 +26,136 @@
#include <triqs/arrays/blas_lapack/gelss.hpp> #include <triqs/arrays/blas_lapack/gelss.hpp>
#include <triqs/python_tools/cython_proxy.hpp> #include <triqs/python_tools/cython_proxy.hpp>
namespace triqs { namespace triqs { namespace gfs { namespace local {
namespace gfs {
namespace local {
using triqs::gfs::imfreq; using triqs::gfs::imfreq;
using triqs::gfs::block_index; using triqs::gfs::block_index;
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 // routine for fitting the tail (singularity) of a Matsubara Green's function
// this is a *linear* least squares problem (with non-linear basis functions) // this is a *linear* least squares problem (with non-linear basis functions)
// which is solved by singular value decomposition of the design matrix // which is solved by singular value decomposition of the design matrix
// the routine fits the real part (even moments) and the imaginary part // the routine fits the real part (even moments) and the imaginary part
//(odd moments) separately, since this is more stable //(odd moments) separately, since this is more stable
// input: // input:
// 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): n_min, n_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 n_min, int n_max) { tail fit_tail_impl(gf_view<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));
if (known_moments.size()) if (known_moments.size())
for (int i = known_moments.order_min(); i <= known_moments.order_max(); i++) res(i) = known_moments(i); 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 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 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(); int n_unknown_moments = n_moments - known_moments.size();
if (n_unknown_moments < 1) return known_moments; 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 // 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 // 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 = 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_even = omin % 2 == 0 ? omin : omin + 1;
int omin_odd = omin % 2 != 0 ? omin : omin + 1; int omin_odd = omin % 2 != 0 ? omin : omin + 1;
int size_even = n_unknown_moments / 2; int size_even = n_unknown_moments / 2;
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 = n_max - n_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);
arrays::matrix<double, 2> B(size1, 1, FORTRAN_LAYOUT); arrays::matrix<double, 2> B(size1, 1, FORTRAN_LAYOUT);
arrays::vector<double> S(std::max(size_even, size_odd)); arrays::vector<double> S(std::max(size_even, size_odd));
const double rcond = 0.0; const double rcond = 0.0;
int rank; int rank;
for (size_t i = 0; i < get_target_shape(gf)[0]; i++) { 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 j = 0; j < get_target_shape(gf)[1]; j++) {
// fit the odd moments // fit the odd moments
// 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 = n_min + k; auto n = n_min + k;
auto iw = std::complex<double>(gf.mesh().index_to_point(n)); auto iw = std::complex<double>(gf.mesh().index_to_point(n));
B(k, 0) = imag(gf.data()(gf.mesh().index_to_linear(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));
for (int l = 0; l < size_odd; l++) { for (int l = 0; l < size_odd; l++) {
int order = omin_odd + 2 * l; int order = omin_odd + 2 * l;
A(k, l) = imag(pow(iw, -1.0 * order)); // set design matrix for odd moments 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<double>(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);
} }
} }
}
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<double>(gf.mesh().index_to_point(n));
void set_tail_from_fit(gf<imfreq> &gf, tail_view known_moments, int n_moments, size_t n_min, size_t n_max, B(k, 0) = real(gf.data()(gf.mesh().index_to_linear(n), i, j));
bool replace_by_fit = false) { // subtract known tail if present
if (get_target_shape(gf) != known_moments.shape()) TRIQS_RUNTIME_ERROR << "shape of tail does not match shape of gf"; if (known_moments.size() > 0)
gf.singularity() = fit_tail_impl(gf, known_moments, n_moments, n_min, n_max); 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 for (int l = 0; l < size_even; l++) {
size_t i = 0; int order = omin_even + 2 * l;
for (auto iw : gf.mesh()) { // (arrays::range(n_min,n_max+1)) { A(k, l) = real(pow(iw, -1.0 * order)); // set design matrix for odd moments
if ((i >= n_min) && (i <= n_max)) gf[iw] = gf.singularity().evaluate(iw); }
i++; }
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<imfreq> 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_index, gf<imfreq>> &block_gf, tail_view known_moments, int n_moments, size_t n_min, void set_tail_from_fit(gf_view<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) { 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(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, n_min, n_max, replace_by_fit); 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<imfreq, scalar_valued> 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 #endif