LCPQ/dft_tools
3
0
Fork 0
mirror of https://github.com/triqs/dft_tools synced 2024-12-26 06:14:14 +01:00
Code Issues Releases Wiki Activity

mc_tools : improve treatment of infinite mc ratios

Browse Source 
- when the ratio returned by an attempt of a move is infinite,
previous code was just throwing TRIQS_RUNTIME_ERROR.

- Now when the ratio is infinite, it is replaced by a large number
  (>1 is enough for metropolis), and the sign is properly updated
  using std::signbit.

- NB :

- a double/float in C++ can be :  normal/ zero/ nan/ infinite / subnormal.
  Here, the code will recover only from infinite case.

- std::signbit works for infinite (according to standard).
This commit is contained in:
Olivier Parcollet 2014-03-27 20:25:13 +01:00
parent 8045598534
commit c996c3ff7d
1 changed files with 13 additions and 7 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
triqs/mc_tools/mc_move_set.hpp
View File

@ -183,14 +183,20 @@ namespace triqs { namespace mc_tools {
std::cerr <<" Proposition probability = "<<proba<<std::endl; std::cerr <<" Proposition probability = "<<proba<<std::endl;
#endif #endif
MCSignType rate_ratio = current->attempt(); MCSignType rate_ratio = current->attempt();
double abs_rate_ratio;
if (std::isinf(rate_ratio)) { // in case the ratio is infinite
abs_rate_ratio = 100; //1.e30; // >1 for metropolis
try_sign_ratio = (std::signbit(rate_ratio) ? -1 : 1); // signbit -> true iif the number is negative
} else {
if (!std::isfinite(std::abs(rate_ratio))) if (!std::isfinite(std::abs(rate_ratio)))
TRIQS_RUNTIME_ERROR<<"Monte Carlo Error : the rate is not finite in move "<<name_of_currently_selected(); TRIQS_RUNTIME_ERROR << "Monte Carlo Error : the rate is not finite in move " << name_of_currently_selected();
double abs_rate_ratio = std::abs(rate_ratio); abs_rate_ratio = std::abs(rate_ratio);
#ifdef TRIQS_TOOLS_MC_DEBUG #ifdef TRIQS_TOOLS_MC_DEBUG
std::cerr << " Metropolis ratio " << rate_ratio<<". Abs(Metropolis ratio) " <<abs_rate_ratio << std::endl; std::cerr << " Metropolis ratio " << rate_ratio << ". Abs(Metropolis ratio) " << abs_rate_ratio << std::endl;
#endif #endif
assert ((abs_rate_ratio>=0)); assert((abs_rate_ratio >= 0));
try_sign_ratio = ( abs_rate_ratio> 1.e-14 ? rate_ratio/abs_rate_ratio : 1); // keep the sign try_sign_ratio = (abs_rate_ratio > 1.e-14 ? rate_ratio / abs_rate_ratio : 1); // keep the sign
}
return abs_rate_ratio; return abs_rate_ratio;
} }