Keep Crowd fit.

This commit is contained in:
Ye Luo 2019-10-24 18:47:01 -05:00
parent 87ce097dbd
commit 6d6a3f1618
5 changed files with 34 additions and 79 deletions

View File

@ -11,14 +11,6 @@
namespace qmcplusplus
{
void Crowd::clearResults()
{
// These were cleared to 1.0 each loop by VMCUpdatePbyP advance walker
// refactored code may depend on this initial value.
std::fill(log_gf_.begin(), log_gf_.end(), 1.0);
std::fill(log_gb_.begin(), log_gb_.end(), 1.0);
}
void Crowd::clearWalkers()
{
mcp_walkers_.clear();
@ -40,18 +32,6 @@ void Crowd::reserve(int crowd_size)
reserveCS(walker_elecs_);
reserveCS(walker_twfs_);
reserveCS(walker_hamiltonians_);
resizeResults(crowd_size);
}
void Crowd::resizeResults(int crowd_size) {
auto resizeCS = [crowd_size](auto& avector) { avector.resize(crowd_size); };
resizeCS(grads_now_);
resizeCS(grads_new_);
resizeCS(ratios_);
resizeCS(log_gf_);
resizeCS(log_gb_);
resizeCS(prob_);
}
void Crowd::addWalker(MCPWalker& walker, ParticleSet& elecs, TrialWaveFunction& twf, QMCHamiltonian& hamiltonian)
@ -61,8 +41,6 @@ void Crowd::addWalker(MCPWalker& walker, ParticleSet& elecs, TrialWaveFunction&
walker_elecs_.push_back(elecs);
walker_twfs_.push_back(twf);
walker_hamiltonians_.push_back(hamiltonian);
if(mcp_walkers_.size() != grads_now_.size())
resizeResults(mcp_walkers_.size());
};
void Crowd::loadWalkers()

View File

@ -54,14 +54,6 @@ public:
void loadWalkers();
/** clears log_gf and log_gb
*
* This is a legacy "call"
* TODO: likely remove log_gf and log_gb from crowd
* seems unlikely these should be persistent state.
*/
void clearResults();
/** Clears all walker vectors
*
* Unless you are _redistributing_ walkers to crowds don't
@ -94,12 +86,6 @@ public:
std::vector<std::reference_wrapper<TrialWaveFunction>>& get_walker_twfs() { return walker_twfs_; }
std::vector<std::reference_wrapper<QMCHamiltonian>>& get_walker_hamiltonians() { return walker_hamiltonians_; }
std::vector<GradType>& get_grads_now() { return grads_now_; }
std::vector<GradType>& get_grads_new() { return grads_new_; }
std::vector<TrialWaveFunction::PsiValueType>& get_ratios() { return ratios_; }
std::vector<RealType>& get_log_gf() { return log_gf_; }
std::vector<RealType>& get_log_gb() { return log_gb_; }
std::vector<RealType>& get_prob() { return prob_; }
RefVector<WFBuffer>& get_mcp_wfbuffers() { return mcp_wfbuffers_; }
const EstimatorManagerCrowd& get_estimator_manager_crowd() const { return estimator_manager_crowd_; }
int size() const { return mcp_walkers_.size(); }
@ -115,7 +101,6 @@ public:
}
unsigned long get_nonlocal_accept() { return n_nonlocal_accept_; }
private:
void resizeResults(int crowd_size);
/** @name Walker Vectors
*
* A single index into these ordered lists constitue a complete
@ -130,19 +115,6 @@ private:
/** }@ */
EstimatorManagerCrowd estimator_manager_crowd_;
/** @name Work Buffers
* @{
* There are many "local" variables in execution of the driver that are convenient to
* place in a stl containers to use <alogorithm> style calls with.
* Eventually this will also us to allow some sort of appropriate parallel policy for these calls.
*/
std::vector<GradType> grads_now_;
std::vector<GradType> grads_new_;
std::vector<TrialWaveFunction::PsiValueType> ratios_;
std::vector<RealType> log_gf_;
std::vector<RealType> log_gb_;
std::vector<RealType> prob_;
/** }@ */
/** @name Step State
*

View File

@ -120,11 +120,19 @@ void DMCBatched::advanceWalkers(const StateForThread& sft,
timers.buffer_timer.stop();
timers.movepbyp_timer.start();
int num_walkers = crowd.size();
const int num_walkers = crowd.size();
//This generates an entire steps worth of deltas.
step_context.nextDeltaRs();
auto it_delta_r = step_context.deltaRsBegin();
std::vector<TrialWaveFunction::GradType> grads_now(num_walkers);
std::vector<TrialWaveFunction::GradType> grads_new(num_walkers);
std::vector<TrialWaveFunction::PsiValueType> ratios(num_walkers);
std::vector<PosType> drifts(num_walkers);
std::vector<RealType> log_gf(num_walkers);
std::vector<RealType> log_gb(num_walkers);
std::vector<RealType> prob(num_walkers);
// local list to handle accept/reject
std::vector<std::reference_wrapper<ParticleSet>> elec_accept_list, elec_reject_list;
@ -160,15 +168,13 @@ void DMCBatched::advanceWalkers(const StateForThread& sft,
int end_index = step_context.getPtclGroupEnd(ig);
for (int iat = start_index; iat < end_index; ++iat)
{
// Correct for DMC as well as VMC?
crowd.clearResults();
ParticleSet::flex_setActive(crowd.get_walker_elecs(), iat);
auto delta_r_start = it_delta_r + iat * num_walkers;
auto delta_r_end = delta_r_start + num_walkers;
//get the displacement
TrialWaveFunction::flex_evalGrad(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, crowd.get_grads_now());
sft.drift_modifier.getDrifts(tauovermass, crowd.get_grads_now(), drifts);
TrialWaveFunction::flex_evalGrad(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, grads_now);
sft.drift_modifier.getDrifts(tauovermass, grads_now, drifts);
std::transform(drifts.begin(), drifts.end(), delta_r_start, drifts.begin(),
[sqrttau](PosType& drift, PosType& delta_r) { return drift + (sqrttau * delta_r); });
@ -188,8 +194,8 @@ void DMCBatched::advanceWalkers(const StateForThread& sft,
auto elecs = crowd.get_walker_elecs();
ParticleSet::flex_makeMove(crowd.get_walker_elecs(), iat, drifts);
TrialWaveFunction::flex_ratioGrad(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, crowd.get_ratios(),
crowd.get_grads_new());
TrialWaveFunction::flex_ratioGrad(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, ratios,
grads_new);
// This lambda is not nested thread safe due to the nreject, nnode_crossing updates
auto checkPhaseChanged = [&sft, &iat, &crowd, &nnode_crossing](TrialWaveFunction& twf, ParticleSet& elec,
@ -215,18 +221,18 @@ void DMCBatched::advanceWalkers(const StateForThread& sft,
//This is just convenient to do here
rr_proposed[iw] += rr[iw];
}
std::transform(delta_r_start, delta_r_end, crowd.get_log_gf().begin(),
std::transform(delta_r_start, delta_r_end, log_gf.begin(),
[mhalf](auto& delta_r) { return mhalf * dot(delta_r, delta_r); });
sft.drift_modifier.getDrifts(tauovermass, crowd.get_grads_new(), drifts);
sft.drift_modifier.getDrifts(tauovermass, grads_new, drifts);
std::transform(crowd.beginElectrons(), crowd.endElectrons(), drifts.begin(), drifts.begin(),
[iat](auto& elecs, auto& drift) { return elecs.get().R[iat] - elecs.get().activePos - drift; });
std::transform(drifts.begin(), drifts.end(), crowd.get_log_gb().begin(),
std::transform(drifts.begin(), drifts.end(), log_gb.begin(),
[oneover2tau](auto& drift) { return -oneover2tau * dot(drift, drift); });
std::transform(crowd.get_ratios().begin(), crowd.get_ratios().end(), crowd.get_prob().begin(),
std::transform(ratios.begin(), ratios.end(), prob.begin(),
[](auto ratio) { return std::norm(ratio); });
twf_accept_list.clear();
@ -237,19 +243,15 @@ void DMCBatched::advanceWalkers(const StateForThread& sft,
for (int iw = 0; iw < num_walkers; ++iw)
{
auto prob = crowd.get_prob()[iw];
auto log_gf = crowd.get_log_gf()[iw];
auto log_gb = crowd.get_log_gb()[iw];
if ((!rejects[iw]) && prob >= std::numeric_limits<RealType>::epsilon() &&
step_context.get_random_gen()() < prob * std::exp(log_gb - log_gf))
if ((!rejects[iw]) && prob[iw] >= std::numeric_limits<RealType>::epsilon() &&
step_context.get_random_gen()() < prob[iw] * std::exp(log_gb[iw] - log_gf[iw]))
{
did_walker_move[iw] += 1;
crowd.incAccept();
twf_accept_list.push_back(crowd.get_walker_twfs()[iw]);
elec_accept_list.push_back(crowd.get_walker_elecs()[iw]);
rr_accepted[iw] += rr[iw];
gf_acc[iw] *= prob;
gf_acc[iw] *= prob[iw];
}
else
{

View File

@ -86,14 +86,19 @@ void VMCBatched::advanceWalkers(const StateForThread& sft,
timers.buffer_timer.stop();
timers.movepbyp_timer.start();
int num_walkers = crowd.size();
const int num_walkers = crowd.size();
// Note std::vector<bool> is not like the rest of stl.
std::vector<bool> moved(num_walkers, false);
constexpr RealType mhalf(-0.5);
bool use_drift = sft.vmcdrv_input.get_use_drift();
const bool use_drift = sft.vmcdrv_input.get_use_drift();
//This generates an entire steps worth of deltas.
step_context.nextDeltaRs();
auto it_delta_r = step_context.deltaRsBegin();
std::vector<TrialWaveFunction::GradType> grads_now(num_walkers);
std::vector<TrialWaveFunction::GradType> grads_new(num_walkers);
std::vector<TrialWaveFunction::PsiValueType> ratios(num_walkers);
std::vector<PosType> drifts(num_walkers);
std::vector<RealType> log_gf(num_walkers);
std::vector<RealType> log_gb(num_walkers);
@ -117,7 +122,6 @@ void VMCBatched::advanceWalkers(const StateForThread& sft,
int end_index = step_context.getPtclGroupEnd(ig);
for (int iat = start_index; iat < end_index; ++iat)
{
crowd.clearResults();
ParticleSet::flex_setActive(crowd.get_walker_elecs(), iat);
// step_context.deltaRsBegin returns an iterator to a flat series of PosTypes
// fastest in walkers then particles
@ -126,8 +130,8 @@ void VMCBatched::advanceWalkers(const StateForThread& sft,
if (use_drift)
{
TrialWaveFunction::flex_evalGrad(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, crowd.get_grads_now());
sft.drift_modifier.getDrifts(tauovermass, crowd.get_grads_now(), drifts);
TrialWaveFunction::flex_evalGrad(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, grads_now);
sft.drift_modifier.getDrifts(tauovermass, grads_now, drifts);
std::transform(drifts.begin(), drifts.end(), delta_r_start, drifts.begin(),
[sqrttau](const PosType& drift, const PosType& delta_r) { return drift + (sqrttau * delta_r); });
@ -143,12 +147,12 @@ void VMCBatched::advanceWalkers(const StateForThread& sft,
// This is inelegant
if (use_drift)
{
TrialWaveFunction::flex_ratioGrad(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, crowd.get_ratios(),
crowd.get_grads_new());
TrialWaveFunction::flex_ratioGrad(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, ratios,
grads_new);
std::transform(delta_r_start, delta_r_end, log_gf.begin(),
[mhalf](const PosType& delta_r) { return mhalf * dot(delta_r, delta_r); });
sft.drift_modifier.getDrifts(tauovermass, crowd.get_grads_new(), drifts);
sft.drift_modifier.getDrifts(tauovermass, grads_new, drifts);
std::transform(crowd.beginElectrons(), crowd.endElectrons(), drifts.begin(), drifts.begin(),
[iat](const ParticleSet& elecs, const PosType& drift) { return elecs.R[iat] - elecs.activePos - drift; });
@ -158,10 +162,10 @@ void VMCBatched::advanceWalkers(const StateForThread& sft,
}
else
{
TrialWaveFunction::flex_calcRatio(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, crowd.get_ratios());
TrialWaveFunction::flex_calcRatio(crowd.get_walker_twfs(), crowd.get_walker_elecs(), iat, ratios);
}
std::transform(crowd.get_ratios().begin(), crowd.get_ratios().end(), prob.begin(),
std::transform(ratios.begin(), ratios.end(), prob.begin(),
[](auto ratio) { return std::norm(ratio); });
twf_accept_list.clear();

View File

@ -147,7 +147,6 @@ TEST_CASE("Crowd redistribute walkers")
crowd.addWalker(*crowd_with_walkers.walkers[iw], *crowd_with_walkers.psets[iw], *crowd_with_walkers.twfs[iw],
*crowd_with_walkers.hams[iw]);
REQUIRE(crowd.size() == 3);
REQUIRE(crowd.get_grads_new().size() == 3);
}
} // namespace qmcplusplus