Skip to content

Split generate and thermalize #167

@Chrisgant

Description

@Chrisgant

The generate function of the HMC class is doing the job of two functions. For better readability it should be split up.

template<Chain chain, NSL::Concept::isNumber Type>
std::conditional_t<chain == Chain::AllStates, std::vector<NSL::MCMC::MarkovState<Type>>, NSL::MCMC::MarkovState<Type>>
generate(const NSL::MCMC::MarkovState<Type> & state, NSL::size_t Nconf, NSL::size_t saveFrequency = 1, std::string baseNode = "markovChain"){
// ensure that saveFrequency is at least 1.
if (saveFrequency <= 0) {
saveFrequency = 1;
}
std::size_t logFrequency = 1;
if(Nconf >= 100){
logFrequency = static_cast<NSL::size_t>( 0.01*Nconf );
}
if constexpr(chain == Chain::AllStates) {
// prepare some memory to all states
std::vector<NSL::MCMC::MarkovState<Type>> MC(Nconf);
// Put the initial configuration in the 0th element
// check if the file already contains configs and if we are allowed to overwrite them
NSL::size_t nstart = 0;
if (h5_.exist(baseNode)){
auto [minConfigID, maxConfigID] = h5_.getMinMaxConfigs(baseNode);
nstart = maxConfigID;
}
if (nstart > 0){
// this reader looks for the most recent markov state and reads it into the state
MC[nstart] = state;
h5_.read(MC[nstart], baseNode);
} else{
MC[nstart] = state;
NSL::Logger::info("HMC: Starting new Markov Chain");
h5_.write(MC[nstart],fmt::format("{}/{}",baseNode,nstart));
}
double runningAcceptance = 0;
// generate Nconf-1 configurations
auto mc_time = NSL::Logger::start_profile("HMC");
for(NSL::size_t n = nstart+1; n < Nconf; ++n){
auto tmp = MC[n-1];
// between each configuration generate saveFrequency which
// are not used for measurements
for(NSL::size_t m = 0; m < saveFrequency-1; ++m){
tmp = this->generate_(tmp);
}
MC[n] = this->generate_(tmp);
h5_.write(MC[n],fmt::format("{}/{}",baseNode,n));
runningAcceptance += static_cast<double>(MC[n].accepted);
// ToDo: have a proper hook being called here
if (n % logFrequency == 0){
NSL::Logger::info("HMC: {}/{}; Running Acceptence Rate: {:.6}%", n, Nconf, runningAcceptance*100./(n-nstart));
NSL::Logger::elapsed_profile(mc_time);
}
}
NSL::Logger::stop_profile(mc_time);
// return the Markov Chain
return MC;
} else {
// for Chain::LastState we only need a new state that becomes overwritten over and over again.
NSL::MCMC::MarkovState<Type> newState = state;
// generate Nconf-1 configurations
// As none is returned we just multiply the number of configurations
for(NSL::size_t n = 1; n < Nconf*saveFrequency; ++n){
newState = this->generate_(newState);
if (n % logFrequency == 0){
NSL::Logger::info("HMC: {}/{}", n, Nconf);
}
}
// return the Markov Chain
return newState;
}
}

Metadata

Metadata

Assignees

No one assigned

    Labels

    enhancementNew feature or request

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions