Skip to content

Commit

Permalink
Merge branch 'isdb_opt_sigma_max'
Browse files Browse the repository at this point in the history
  • Loading branch information
@carlocamilloni
carlocamilloni committed Jul 15, 2020
2 parents 5580837 + 2ab5a63 commit c19727d
Show file tree
Hide file tree
Showing 3 changed files with 108 additions and 9 deletions.
60 changes: 55 additions & 5 deletions src/isdb/Metainference.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -241,6 +241,11 @@ class Metainference : public bias::Bias
// optimize sigma mean
vector< vector < vector <double> > > sigma_mean2_last_;
unsigned optsigmamean_stride_;
// optimize sigma max
unsigned N_optimized_step_;
unsigned optimized_step_;
bool sigmamax_opt_done_;
std::vector<double> sigma_max_est_;

// average weights
unsigned average_weights_stride_;
Expand Down Expand Up @@ -312,6 +317,7 @@ void Metainference::registerKeywords(Keywords& keys) {
keys.add("optional","DSIGMA","maximum MC move of the uncertainty parameter");
keys.add("compulsory","OPTSIGMAMEAN","NONE","Set to NONE/SEM to manually set sigma mean, or to estimate it on the fly");
keys.add("optional","SIGMA_MEAN0","starting value for the uncertainty in the mean estimate");
keys.add("optional","SIGMA_MAX_STEPS", "Number of steps used to optimise SIGMA_MAX, before that the SIGMA_MAX value is used");
keys.add("optional","TEMP","the system temperature - this is only needed if code doesn't pass the temperature to plumed");
keys.add("optional","MC_STEPS","number of MC steps");
keys.add("optional","MC_CHUNKSIZE","MC chunksize");
Expand Down Expand Up @@ -361,6 +367,9 @@ Metainference::Metainference(const ActionOptions&ao):
do_reweight_(false),
do_optsigmamean_(0),
optsigmamean_stride_(0),
N_optimized_step_(0),
optimized_step_(0),
sigmamax_opt_done_(false),
average_weights_stride_(1)
{
bool noensemble = false;
Expand Down Expand Up @@ -451,14 +460,21 @@ Metainference::Metainference(const ActionOptions&ao):
parse("OPTSIGMAMEAN", stringa_optsigma);
if(stringa_optsigma=="NONE") do_optsigmamean_=0;
else if(stringa_optsigma=="SEM") do_optsigmamean_=1;
else if(stringa_optsigma=="SEM_MAX") do_optsigmamean_=2;

unsigned aver_max_steps=0;
parse("SIGMA_MAX_STEPS", aver_max_steps);
if(aver_max_steps==0&&do_optsigmamean_==2) aver_max_steps=averaging*1000;
if(aver_max_steps>0&&do_optsigmamean_<2) error("SIGMA_MAX_STEPS can only be used together with OPTSIGMAMEAN=SEM_MAX");
if(aver_max_steps>0&&do_optsigmamean_==2) N_optimized_step_=aver_max_steps;

// resize vector for sigma_mean history
sigma_mean2_last_.resize(nsel);
for(unsigned i=0; i<nsel; i++) sigma_mean2_last_[i].resize(narg);

vector<double> read_sigma_mean_;
parseVector("SIGMA_MEAN0",read_sigma_mean_);
if(!do_optsigmamean_ && read_sigma_mean_.size()==0 && !getRestart())
if(do_optsigmamean_==0 && read_sigma_mean_.size()==0 && !getRestart())
error("If you don't use OPTSIGMAMEAN and you are not RESTARTING then you MUST SET SIGMA_MEAN0");

if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
Expand Down Expand Up @@ -572,8 +588,8 @@ Metainference::Metainference(const ActionOptions&ao):
Dsigma_.resize(read_dsigma.size());
Dsigma_=read_dsigma;
} else {
Dsigma_.resize(sigma_max_.size());
for(unsigned i=0; i<sigma_max_.size(); i++) Dsigma_[i] = 0.05*(sigma_max_[i] - sigma_min_[i]);
Dsigma_.resize(sigma_max_.size(), -1.);
/* in this case Dsigma is initialised after reading the restart file if present */
}

// monte carlo stuff
Expand Down Expand Up @@ -616,6 +632,8 @@ Metainference::Metainference(const ActionOptions&ao):
}
}

sigma_max_est_.resize(sigma_max_.size(), 0.);

IFile restart_sfile;
restart_sfile.link(*this);
if(getRestart()&&restart_sfile.FileExist(status_file_name_)) {
Expand Down Expand Up @@ -666,6 +684,12 @@ Metainference::Metainference(const ActionOptions&ao):
Tools::convert(i,msg);
restart_sfile.scanField("sigma_"+msg,sigma_[i]);
}
for(unsigned i=0; i<sigma_max_.size(); ++i) {
std::string msg;
Tools::convert(i,msg);
restart_sfile.scanField("sigma_max_"+msg,sigma_max_[i]);
sigmamax_opt_done_=true;
}
if(noise_type_==GENERIC) {
for(unsigned i=0; i<ftilde_.size(); ++i) {
std::string msg;
Expand Down Expand Up @@ -693,6 +717,9 @@ Metainference::Metainference(const ActionOptions&ao):
restart_sfile.close();
}

/* If DSIGMA is not yet initialised do it now */
for(unsigned i=0; i<sigma_max_.size(); i++) if(Dsigma_[i]==-1) Dsigma_[i] = 0.05*(sigma_max_[i] - sigma_min_[i]);

switch(noise_type_) {
case GENERIC:
log.printf(" with general metainference ");
Expand Down Expand Up @@ -1490,16 +1517,17 @@ void Metainference::get_weights(const unsigned iselect, double &fact, double &va
vector<double> bias(nrep_,0);
if(master) {
bias[replica_] = getArgument(narg);
//bias[replica_] = ((1.0/plumed.getAtoms().getKbT())- (1.0/kbt_) )*getArgument(narg);
if(nrep_>1) multi_sim_comm.Sum(&bias[0], nrep_);
}
comm.Sum(&bias[0], nrep_);

const double maxbias = *(std::max_element(bias.begin(), bias.end()));
for(unsigned i=0; i<nrep_; ++i) {
bias[i] = exp((bias[i]-maxbias)/kbt_);
//bias[i] = exp((bias[i]-maxbias));
norm += bias[i];
}

// accumulate weights
const double decay = 1./static_cast<double> (average_weights_stride_);
if(!firstTimeW[iselect]) {
Expand Down Expand Up @@ -1594,7 +1622,24 @@ void Metainference::get_sigma_mean(const unsigned iselect, const double fact, co
sigma_mean2_tmp[0] = max_now;
valueSigmaMean[0]->set(sqrt(sigma_mean2_tmp[0]));
}
// endif sigma optimization
// endif sigma mean optimization
// start sigma max optimization
if(do_optsigmamean_>1&&!sigmamax_opt_done_) {
for(unsigned i=0; i<sigma_max_.size(); i++) {
sigma_max_est_[i] += sqrt(sigma_mean2_tmp[i]);
// ready to set once and for all the value of sigma_max
if(optimized_step_==N_optimized_step_) {
double isteps = 1./static_cast<double>(optimized_step_);
sigmamax_opt_done_=true;
for(unsigned i=0; i<sigma_max_.size(); i++) {
sigma_max_[i]=sigma_max_est_[i]*isteps*sqrt(dnrep);
Dsigma_[i] = 0.05*(sigma_max_[i] - sigma_min_[i]);
}
}
}
optimized_step_++;
}
// end sigma max optimization
} else {
if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
for(unsigned i=0; i<narg; ++i) {
Expand Down Expand Up @@ -1726,6 +1771,11 @@ void Metainference::writeStatus()
Tools::convert(i,msg);
sfile_.printField("sigma_"+msg,sigma_[i]);
}
for(unsigned i=0; i<sigma_max_.size(); ++i) {
std::string msg;
Tools::convert(i,msg);
sfile_.printField("sigma_max_"+msg,sigma_max_[i]);
}
if(noise_type_==GENERIC) {
for(unsigned i=0; i<ftilde_.size(); ++i) {
std::string msg;
Expand Down
52 changes: 48 additions & 4 deletions src/isdb/MetainferenceBase.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@ void MetainferenceBase::registerKeywords( Keywords& keys ) {
keys.add("optional","DSIGMA","maximum MC move of the uncertainty parameter");
keys.add("compulsory","OPTSIGMAMEAN","NONE","Set to NONE/SEM to manually set sigma mean, or to estimate it on the fly");
keys.add("optional","SIGMA_MEAN0","starting value for the uncertainty in the mean estimate");
keys.add("optional","SIGMA_MAX_STEPS", "Number of steps used to optimise SIGMA_MAX, before that the SIGMA_MAX value is used");
keys.add("optional","TEMP","the system temperature - this is only needed if code doesn't pass the temperature to plumed");
keys.add("optional","MC_STEPS","number of MC steps");
keys.add("optional","MC_CHUNKSIZE","MC chunksize");
Expand Down Expand Up @@ -124,6 +125,9 @@ MetainferenceBase::MetainferenceBase(const ActionOptions&ao):
nsel_(1),
iselect(0),
optsigmamean_stride_(0),
N_optimized_step_(0),
optimized_step_(0),
sigmamax_opt_done_(false),
decay_w_(1.)
{
parseFlag("DOSCORE", doscore_);
Expand Down Expand Up @@ -195,10 +199,17 @@ MetainferenceBase::MetainferenceBase(const ActionOptions&ao):
parse("OPTSIGMAMEAN", stringa_optsigma);
if(stringa_optsigma=="NONE") do_optsigmamean_=0;
else if(stringa_optsigma=="SEM") do_optsigmamean_=1;
else if(stringa_optsigma=="SEM_MAX") do_optsigmamean_=2;

unsigned aver_max_steps=0;
parse("SIGMA_MAX_STEPS", aver_max_steps);
if(aver_max_steps==0&&do_optsigmamean_==2) aver_max_steps=averaging*1000;
if(aver_max_steps>0&&do_optsigmamean_<2) error("SIGMA_MAX_STEPS can only be used together with OPTSIGMAMEAN=SEM_MAX");
if(aver_max_steps>0&&do_optsigmamean_==2) N_optimized_step_=aver_max_steps;

vector<double> read_sigma_mean_;
parseVector("SIGMA_MEAN0",read_sigma_mean_);
if(!do_optsigmamean_ && read_sigma_mean_.size()==0 && !getRestart() && doscore_)
if(do_optsigmamean_==0 && read_sigma_mean_.size()==0 && !getRestart() && doscore_)
error("If you don't use OPTSIGMAMEAN and you are not RESTARTING then you MUST SET SIGMA_MEAN0");

if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
Expand Down Expand Up @@ -304,8 +315,8 @@ MetainferenceBase::MetainferenceBase(const ActionOptions&ao):
Dsigma_.resize(read_dsigma.size());
Dsigma_=read_dsigma;
} else {
Dsigma_.resize(sigma_max_.size());
for(unsigned i=0; i<sigma_max_.size(); i++) Dsigma_[i] = 0.05*(sigma_max_[i] - sigma_min_[i]);
Dsigma_.resize(sigma_max_.size(), -1.);
/* in this case Dsigma is initialised after reading the restart file if present */
}

// monte carlo stuff
Expand Down Expand Up @@ -417,6 +428,8 @@ void MetainferenceBase::Initialise(const unsigned input)
}
}

sigma_max_est_.resize(sigma_max_.size(), 0.);

IFile restart_sfile;
restart_sfile.link(*this);
if(getRestart()&&restart_sfile.FileExist(status_file_name_)) {
Expand Down Expand Up @@ -467,6 +480,12 @@ void MetainferenceBase::Initialise(const unsigned input)
Tools::convert(i,msg);
restart_sfile.scanField("sigma_"+msg,sigma_[i]);
}
for(unsigned i=0; i<sigma_max_.size(); ++i) {
std::string msg;
Tools::convert(i,msg);
restart_sfile.scanField("sigma_max_"+msg,sigma_max_[i]);
sigmamax_opt_done_=true;
}
if(noise_type_==GENERIC) {
for(unsigned i=0; i<ftilde_.size(); ++i) {
std::string msg;
Expand Down Expand Up @@ -494,6 +513,9 @@ void MetainferenceBase::Initialise(const unsigned input)
restart_sfile.close();
}

/* If DSIGMA is not yet initialised do it now */
for(unsigned i=0; i<sigma_max_.size(); i++) if(Dsigma_[i]==-1) Dsigma_[i] = 0.05*(sigma_max_[i] - sigma_min_[i]);

addComponentWithDerivatives("score");
componentIsNotPeriodic("score");
valueScore=getPntrToComponent("score");
Expand Down Expand Up @@ -1367,7 +1389,24 @@ void MetainferenceBase::get_sigma_mean(const double fact, const double var_fact,
sigma_mean2_tmp[0] = max_now;
valueSigmaMean[0]->set(sqrt(sigma_mean2_tmp[0]));
}
// endif sigma optimization
// endif sigma mean optimization
// start sigma max optimization
if(do_optsigmamean_>1&&!sigmamax_opt_done_) {
for(unsigned i=0; i<sigma_max_.size(); i++) {
sigma_max_est_[i] += sqrt(sigma_mean2_tmp[i]);
// ready to set once and for all the value of sigma_max
if(optimized_step_==N_optimized_step_) {
double isteps = 1./static_cast<double>(optimized_step_);
sigmamax_opt_done_=true;
for(unsigned i=0; i<sigma_max_.size(); i++) {
sigma_max_[i]=sigma_max_est_[i]*isteps*sqrt(dnrep);
Dsigma_[i] = 0.05*(sigma_max_[i] - sigma_min_[i]);
}
}
}
optimized_step_++;
}
// end sigma max optimization
} else {
if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
for(unsigned i=0; i<narg; ++i) {
Expand Down Expand Up @@ -1495,6 +1534,11 @@ void MetainferenceBase::writeStatus()
Tools::convert(i,msg);
sfile_.printField("sigma_"+msg,sigma_[i]);
}
for(unsigned i=0; i<sigma_max_.size(); ++i) {
std::string msg;
Tools::convert(i,msg);
sfile_.printField("sigma_max_"+msg,sigma_max_[i]);
}
if(noise_type_==GENERIC) {
for(unsigned i=0; i<ftilde_.size(); ++i) {
std::string msg;
Expand Down
5 changes: 5 additions & 0 deletions src/isdb/MetainferenceBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,6 +141,11 @@ class MetainferenceBase :
// optimize sigma mean
std::vector< std::vector < std::vector <double> > > sigma_mean2_last_;
unsigned optsigmamean_stride_;
// optimize sigma max
unsigned N_optimized_step_;
unsigned optimized_step_;
bool sigmamax_opt_done_;
std::vector<double> sigma_max_est_;

// average weights
double decay_w_;
Expand Down

0 comments on commit c19727d

Please sign in to comment.