diff --git a/TODO.md b/TODO.md index 7163ce6fe..9a55f5570 100644 --- a/TODO.md +++ b/TODO.md @@ -1,3 +1,5 @@ * Test all kernels and close #3 accordingly -* Add automatic Python tests chi_* and chi_auto_* -* Make adding new kernels more straightforward +* Add automatic Python tests chi_*, chi_auto_* and zerotemp_* +* Fix tests to compare tails of g_w +* Make sure all tests run successfully with versions 0.9 and 0.95 +* Update year in copyright headers (2016-2017) diff --git a/c++/kernels/base.hpp b/c++/kernels/base.hpp index da0458718..8534ced0a 100644 --- a/c++/kernels/base.hpp +++ b/c++/kernels/base.hpp @@ -20,9 +20,12 @@ ******************************************************************************/ #pragma once +#include #include #include #include +#include +#include #include "../rectangle.hpp" #include "../configuration.hpp" @@ -31,7 +34,15 @@ namespace som { // Kinds of observables -enum observable_kind {FermionGf, BosonCorr, BosonAutoCorr, ZeroTemp}; +#define ALL_OBSERVABLES (FermionGf)(BosonCorr)(BosonAutoCorr)(ZeroTemp) + +enum observable_kind : unsigned int {BOOST_PP_SEQ_ENUM(ALL_OBSERVABLES)}; +constexpr const unsigned int n_observable_kinds = BOOST_PP_SEQ_SIZE(ALL_OBSERVABLES); + +// Is statistics defined for this observable? +bool is_stat_relevant(observable_kind kind) { + return kind != ZeroTemp; +} // Statistics of observables triqs::gfs::statistic_enum observable_statistics(observable_kind kind) { @@ -54,6 +65,47 @@ std::string observable_name(observable_kind kind) { } } +// Widest energy windows for observables +std::pair max_energy_window(observable_kind kind) { + switch(kind) { + case FermionGf: return std::make_pair(-HUGE_VAL,HUGE_VAL); + case BosonCorr: return std::make_pair(-HUGE_VAL,HUGE_VAL); + case BosonAutoCorr: return std::make_pair(0,HUGE_VAL); + case ZeroTemp: return std::make_pair(0,HUGE_VAL); + } +} + +// Construct a real-frequency GF from a configuration +void back_transform(observable_kind kind, + configuration const& conf, + cache_index & ci, + triqs::gfs::gf_view g_w) { + bool bosoncorr = kind == BosonCorr || kind == BosonAutoCorr; + + g_w() = 0; + auto & tail = g_w.singularity(); + + for(auto const& rect : conf) { + for(auto e : g_w.mesh()) g_w.data()(e.index()) += rect.hilbert_transform(double(e), bosoncorr); + tail.data()(range(),0,0) += rect.tail_coefficients(tail.order_min(),tail.order_max(), bosoncorr); + + if(kind == BosonAutoCorr) { + // Add a reflected rectangle + rectangle reflected_rect(-rect.center, rect.width, rect.height, ci); + for(auto e : g_w.mesh()) g_w.data()(e.index()) += reflected_rect.hilbert_transform(double(e), true); + tail.data()(range(),0,0) += reflected_rect.tail_coefficients(tail.order_min(),tail.order_max(), true); + } + + if(bosoncorr) { + g_w.data() *= -1.0/M_PI; + tail.data() *= -1.0/M_PI; + } + } +} + +// All meshes used with input data containers +#define ALL_INPUT_MESHES (imtime)(imfreq)(legendre) + // Mesh traits template struct mesh_traits; diff --git a/c++/som_core.cpp b/c++/som_core.cpp index 28f67f562..d6342769b 100644 --- a/c++/som_core.cpp +++ b/c++/som_core.cpp @@ -19,6 +19,8 @@ * ******************************************************************************/ #include +#include +#include #include "som_core.hpp" #include "kernels/fermiongf_imtime.hpp" @@ -119,10 +121,7 @@ som_core::som_core(gf_const_view g_tau, gf_const_view S_tau, mesh(g_tau.mesh()), kind(kind), norms(make_default_norms(norms,get_target_shape(g_tau)[0])), rhs(input_data_r_t()), error_bars(input_data_r_t()) { - if(kind != FermionGf && kind != BosonCorr && kind != BosonAutoCorr && kind != ZeroTemp) - fatal_error("unknown observable kind " + to_string(kind)); - - if(kind != ZeroTemp) check_gf_stat(g_tau, observable_statistics(kind)); + if(is_stat_relevant(kind)) check_gf_stat(g_tau, observable_statistics(kind)); check_input_gf(g_tau,S_tau); if(!is_gf_real(g_tau) || !is_gf_real(S_tau)) @@ -137,10 +136,7 @@ som_core::som_core(gf_const_view g_iw, gf_const_view S_iw, mesh(g_iw.mesh()), kind(kind), norms(make_default_norms(norms,get_target_shape(g_iw)[0])), rhs(input_data_c_t()), error_bars(input_data_c_t()) { - if(kind != FermionGf && kind != BosonCorr && kind != BosonAutoCorr && kind != ZeroTemp) - fatal_error("unknown observable kind " + to_string(kind)); - - if(kind != ZeroTemp) check_gf_stat(g_iw, observable_statistics(kind)); + if(is_stat_relevant(kind)) check_gf_stat(g_iw, observable_statistics(kind)); check_input_gf(g_iw,S_iw); if(!is_gf_real_in_tau(g_iw) || !is_gf_real_in_tau(S_iw)) @@ -157,10 +153,7 @@ som_core::som_core(gf_const_view g_l, gf_const_view S_l, mesh(g_l.mesh()), kind(kind), norms(make_default_norms(norms,get_target_shape(g_l)[0])), rhs(input_data_r_t()), error_bars(input_data_r_t()) { - if(kind != FermionGf && kind != BosonCorr && kind != BosonAutoCorr && kind != ZeroTemp) - fatal_error("unknown observable kind " + to_string(kind)); - - if(kind != ZeroTemp) check_gf_stat(g_l, observable_statistics(kind)); + if(is_stat_relevant(kind)) check_gf_stat(g_l, observable_statistics(kind)); check_input_gf(g_l,S_l); if(!is_gf_real(g_l) || !is_gf_real(S_l)) @@ -177,9 +170,17 @@ void som_core::run(run_parameters_t const& p) { params = p; - if((kind == BosonAutoCorr || kind == ZeroTemp) && params.energy_window.first < 0) { - params.energy_window.first = 0; - if(params.verbosity > 0) warning("left boundary of the energy window is reset to 0"); + double e_min, e_max; + std::tie(e_min,e_max) = max_energy_window(kind); + if(params.energy_window.first < e_min) { + params.energy_window.first = e_min; + if(params.verbosity > 0) + warning("left boundary of the energy window is reset to " + std::to_string(e_min)); + } + if(params.energy_window.second > e_max) { + params.energy_window.second = e_max; + if(params.verbosity > 0) + warning("right boundary of the energy window is reset to " + std::to_string(e_max)); } if(params.energy_window.first >= params.energy_window.second) @@ -202,20 +203,12 @@ void som_core::run(run_parameters_t const& p) { triqs::signal_handler::start(); run_status = 0; try { - #define RUN_IMPL_CASE(ok, mk) case (int(ok) + 4 * mesh_traits::index): run_impl>(); break; - switch(int(kind) + 4 * mesh.index()) { - RUN_IMPL_CASE(FermionGf,imtime); - RUN_IMPL_CASE(FermionGf,imfreq); - RUN_IMPL_CASE(FermionGf,legendre); - RUN_IMPL_CASE(BosonCorr,imtime); - RUN_IMPL_CASE(BosonCorr,imfreq); - RUN_IMPL_CASE(BosonCorr,legendre); - RUN_IMPL_CASE(BosonAutoCorr,imtime); - RUN_IMPL_CASE(BosonAutoCorr,imfreq); - RUN_IMPL_CASE(BosonAutoCorr,legendre); - RUN_IMPL_CASE(ZeroTemp,imtime); - RUN_IMPL_CASE(ZeroTemp,imfreq); - RUN_IMPL_CASE(ZeroTemp,legendre); + #define RUN_IMPL_CASE(r, okmk) \ + case (int(BOOST_PP_SEQ_ELEM(0,okmk)) + \ + n_observable_kinds * mesh_traits::index): \ + run_impl>(); break; + switch(int(kind) + n_observable_kinds * mesh.index()) { + BOOST_PP_SEQ_FOR_EACH_PRODUCT(RUN_IMPL_CASE, (ALL_OBSERVABLES)(ALL_INPUT_MESHES)) } #undef RUN_IMPL_CASE } catch(stopped & e) { @@ -461,65 +454,30 @@ template void triqs_gf_view_assign_delegation(gf_view g, som_core const& cont) { auto gf_dim = cont.results.size(); check_gf_dim(g, gf_dim); - if(cont.kind != ZeroTemp) check_gf_stat(g, observable_statistics(cont.kind)); + if(is_stat_relevant(cont.kind)) check_gf_stat(g, observable_statistics(cont.kind)); g() = 0; +#define FILL_DATA_CASE(r, data, ok) \ + case ok: { \ + kernel kern(g.mesh()); \ + for(int i : range(gf_dim)) fill_data(g, i, kern(cont.results[i])); \ + return; \ + } switch(cont.kind) { - case FermionGf: { - kernel kern(g.mesh()); - for(int i : range(gf_dim)) fill_data(g, i, kern(cont.results[i])); - return; - } - case BosonCorr: { - kernel kern(g.mesh()); - for(int i : range(gf_dim)) fill_data(g, i, kern(cont.results[i])); - return; - } - case BosonAutoCorr: { - kernel kern(g.mesh()); - for(int i : range(gf_dim)) fill_data(g, i, kern(cont.results[i])); - return; - } - case ZeroTemp: { - kernel kern(g.mesh()); - for(int i : range(gf_dim)) fill_data(g, i, kern(cont.results[i])); - return; - } - default: - fatal_error("unknown observable kind " + to_string(cont.kind)); + BOOST_PP_SEQ_FOR_EACH(FILL_DATA_CASE, _, ALL_OBSERVABLES) + default: fatal_error("unknown observable kind " + to_string(cont.kind)); } +#undef FILL_DATA_CASE } template<> void triqs_gf_view_assign_delegation(gf_view g_w, som_core const& cont) { - - bool bosoncorr = cont.kind == BosonCorr || cont.kind == BosonAutoCorr; - if(cont.kind != FermionGf && !bosoncorr && cont.kind != ZeroTemp) - fatal_error("unknown observable kind " + to_string(cont.kind)); auto gf_dim = cont.results.size(); check_gf_dim(g_w, gf_dim); - - g_w() = 0; - auto & tail = g_w.singularity(); - - for(int i : range(gf_dim)) { - auto const& conf = cont.results[i]; - for(auto const& rect : conf) { - for(auto e : g_w.mesh()) g_w.data()(e.index(),i,i) += rect.hilbert_transform(double(e), bosoncorr); - tail.data()(range(),i,i) += rect.tail_coefficients(tail.order_min(),tail.order_max(), bosoncorr); - - if(cont.kind == BosonAutoCorr) { - // Add a reflected rectangle - rectangle reflected_rect(-rect.center, rect.width, rect.height, const_cast(cont).ci); - for(auto e : g_w.mesh()) g_w.data()(e.index(),i,i) += reflected_rect.hilbert_transform(double(e), true); - tail.data()(range(),i,i) += reflected_rect.tail_coefficients(tail.order_min(),tail.order_max(), true); - } - } - - if(bosoncorr) { - g_w.data()(range(),i,i) *= -1.0/M_PI; - tail.data()(range(),i,i) *= -1.0/M_PI; - } - } + for(int i : range(gf_dim)) + back_transform(cont.kind, + cont.results[i], + const_cast(cont).ci, + slice_target_to_scalar(g_w, i, i)); } template void triqs_gf_view_assign_delegation(gf_view, som_core const&); diff --git a/test/python/gf_imfreq.py b/test/python/gf_imfreq.py index 651fed8e9..c52a6d65d 100644 --- a/test/python/gf_imfreq.py +++ b/test/python/gf_imfreq.py @@ -57,5 +57,6 @@ # arch['histograms'] = cont.histograms assert_gfs_are_close(g_rec_iw, arch['g_rec_iw']) assert_gfs_are_close(g_w, arch['g_w'], 6e-4) + assert_arrays_are_close(g_w.tail.data, arch['g_w'].tail.data, 6e-4) assert_arrays_are_close(cont.histograms[0].data, arch['histograms'][0].data) assert_arrays_are_close(cont.histograms[1].data, arch['histograms'][1].data) diff --git a/test/python/gf_imtime.py b/test/python/gf_imtime.py index f203b3f0f..a1a781dac 100644 --- a/test/python/gf_imtime.py +++ b/test/python/gf_imtime.py @@ -57,5 +57,6 @@ # arch['histograms'] = cont.histograms assert_gfs_are_close(g_rec_tau, arch['g_rec_tau']) assert_gfs_are_close(g_w, arch['g_w'], 1e-5) + assert_arrays_are_close(g_w.tail.data, arch['g_w'].tail.data, 1e-5) assert_arrays_are_close(cont.histograms[0].data, arch['histograms'][0].data) assert_arrays_are_close(cont.histograms[1].data, arch['histograms'][1].data) diff --git a/test/python/gf_legendre.py b/test/python/gf_legendre.py index b4d1ac076..679a494bd 100644 --- a/test/python/gf_legendre.py +++ b/test/python/gf_legendre.py @@ -57,5 +57,6 @@ # arch['histograms'] = cont.histograms assert_gfs_are_close(g_rec_l, arch['g_rec_l']) assert_gfs_are_close(g_w, arch['g_w'], 1e-5) + assert_arrays_are_close(g_w.tail.data, arch['g_w'].tail.data, 1e-5) assert_arrays_are_close(cont.histograms[0].data, arch['histograms'][0].data) assert_arrays_are_close(cont.histograms[1].data, arch['histograms'][1].data)