MonteCarlo.cc

Go to the documentation of this file.

#include "casm/monte_carlo/MonteCarlo.hh"
 
#include "casm/monte_carlo/MonteCounter.hh"
#include "casm/monte_carlo/MonteSampler.hh"
#include "casm/monte_carlo/MonteSettings.hh"
 
namespace CASM {
namespace Monte {
 
bool SamplerNameCompare::operator()(const std::string &A,
                                    const std::string &B) const {
  std::string::size_type Apos1 = A.find_first_of("([");
  std::string::size_type Bpos1 = B.find_first_of("([");
  if (A.substr(0, Apos1) == B.substr(0, Bpos1)) {
    std::string::size_type Apos2 = A.find_first_of("])");
    std::string::size_type Bpos2 = B.find_first_of("])");
 
    std::string Aindex = A.substr(Apos1 + 1, Apos2 - Apos1 - 1);
    std::string Bindex = B.substr(Bpos1 + 1, Bpos2 - Bpos1 - 1);
 
    for (int i = 0; i < Aindex.size(); i++) {
      if (!std::isdigit(Aindex[i])) {
        return Aindex < Bindex;
      }
    }
 
    return std::stoi(Aindex) < std::stoi(Bindex);
  }
  return A.substr(0, Apos1) < B.substr(0, Bpos1);
}
 
void MonteCarlo::sample_data(const MonteCounter &counter) {
  // call Sampler::sample(*this) for all samplers
  for (auto it = m_sampler.begin(); it != m_sampler.end(); ++it) {
    it->second->sample(*this, counter);
  }
  m_sample_time.push_back(std::make_pair(counter.pass(), counter.step()));
 
  if (m_write_trajectory) {
    m_trajectory.push_back(configdof());
  }
 
  m_is_equil_uptodate = false;
  m_is_converged_uptodate = false;
}
 
void MonteCarlo::clear_samples() {
  for (auto it = m_sampler.begin(); it != m_sampler.end(); ++it) {
    it->second->clear();
  }
  m_trajectory.clear();
  m_sample_time.clear();
 
  m_is_equil_uptodate = false;
  m_is_converged_uptodate = false;
  m_next_convergence_check = m_convergence_check_period;
}
 
std::pair<bool, Monte::size_type> MonteCarlo::is_equilibrated() const {
  if (m_is_equil_uptodate) {
    return m_is_equil;
  }
 
  m_is_equil_uptodate = true;
 
  _log().check<Log::verbose>("Equilibration");
  _log() << std::boolalpha << std::setw(24) << "quantity" << std::setw(20)
         << "is_equilibrated" << std::setw(16) << "at_sample" << std::endl;
  // check if all samplers that must converge have equilibrated, and find the
  // maximum
  //   of the number of samples needed to equilibrate
  size_type max_equil_samples = 0;
  for (auto it = m_sampler.cbegin(); it != m_sampler.cend(); ++it) {
    if (it->second->must_converge()) {
      // is_equilibrated returns std::pair(is_equilibrated?, equil_samples)
      auto equil = it->second->is_equilibrated();
 
      _log() << std::setw(24) << it->second->name() << std::setw(20)
             << (equil.first ? "true"
                             : "false");  // why isn't boolapha working?
      if (equil.first) {
        _log() << std::setw(16) << equil.second << std::endl;
      } else {
        _log() << std::setw(16) << "unknown" << std::endl;
      }
 
      if (!equil.first) {
        return m_is_equil = std::make_pair(false, size_type(0));
      }
      if (equil.second > max_equil_samples) {
        max_equil_samples = equil.second;
      }
    } else {
      //        _log() << std::setw(24) << it->second->name()
      //               << std::setw(20) << "unknown"
      //               << std::setw(16) << "unknown" << std::endl;
    }
  }
 
  _log() << "Overall equilibration at sample: " << max_equil_samples << "\n"
         << std::endl;
 
  return m_is_equil = std::make_pair(true, max_equil_samples);
}
 
bool MonteCarlo::is_converged() const {
  // if we've already calculated convergence, return result
  if (m_is_converged_uptodate) {
    return m_is_converged;
  }
 
  m_is_converged_uptodate = true;
 
  // set the next time for a convergence check
  _set_check_convergence_time();
 
  if (!m_sampler.size()) {
    m_is_converged = false;
    return m_is_converged;
  }
 
  // check if required equilibration has occured, and get the number of samples
  // required to equilibrate
  auto equil = is_equilibrated();
  if (!equil.first) {
    m_is_converged = false;
    return m_is_converged;
  }
 
  _log().check<Log::verbose>("Convergence");
  _log() << std::boolalpha << std::setw(24) << "quantity" << std::setw(16)
         << "mean" << std::setw(16) << "req_prec" << std::setw(16)
         << "calc_prec" << std::setw(16) << "is_converged" << std::endl;
 
  // check if all samplers that must converge have converged using data in range
  // [max_equil_samples, end)
  m_is_converged = std::all_of(
      m_sampler.cbegin(), m_sampler.cend(),
      [=](const SamplerMap::value_type &val) {
        if (val.second->must_converge()) {
          bool result = val.second->is_converged(equil.second);
          _log() << std::setw(24) << val.second->name() << std::setw(16)
                 << val.second->mean(equil.second) << std::setw(16)
                 << val.second->requested_precision() << std::setw(16)
                 << val.second->calculated_precision(equil.second)
                 << std::setw(16)
                 << (result ? "true" : "false")  // why isn't boolapha working?
                 << std::endl;
          return result;
        } else {
          //        _log() << std::setw(24) << val.second->name()
          //               << std::setw(16) << val.second->mean(equil.second)
          //               << std::setw(16) << "none"
          //               << std::setw(16) << "unknown"
          //               << std::setw(16) << "unknown"
          //               << std::endl;
          return true;
        }
      });
 
  _log() << "Overall convergence?: " << std::boolalpha << m_is_converged << "\n"
         << std::endl;
 
  return m_is_converged;
}
 
bool MonteCarlo::check_convergence_time() const {
  if (m_sample_time.size() >= m_next_convergence_check) {
    return true;
  }
  return false;
}
 
void MonteCarlo::_set_check_convergence_time() const {
  m_next_convergence_check = m_sample_time.size() + m_convergence_check_period;
}
 
}  // namespace Monte
}  // namespace CASM