Skip to content
Snippets Groups Projects
Commit 12dd3035 authored by Gareth Tribello's avatar Gareth Tribello
Browse files

Using new weight actions in analysis

parent 20cd869f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 71 additions and 220 deletions
regtest/analysis/rt-mds/analysis.0.list_embed.reference
+ 1
− 1
View file @ 12dd3035
#! FIELDS @17.1 @17.2
#! FIELDS @18.1 @18.2
0.0873 0.0013
0.0535 0.0044
0.0150 0.0044
......
regtest/analysis/rt-mds/list_embed.reference
+ 1
− 1
View file @ 12dd3035
#! FIELDS @17.1 @17.2
#! FIELDS @18.1 @18.2
0.0594 0.0053
0.0683 0.0034
0.0705 0.0006
......
regtest/analysis/rt-mds/plumed.dat
+ 4
− 1
View file @ 12dd3035
......@@ -16,13 +16,16 @@ DISTANCE ATOMS=7,com LABEL=d7
UPPER_WALLS ARG=d7 AT=2.0 KAPPA=100.
COORDINATIONNUMBER SPECIES=1-7 MOMENTS=2-3 SWITCH={RATIONAL R_0=1.5 NN=8} LABEL=c1
ww: REWEIGHT_TEMP REWEIGHT_TEMP=0.1 TEMP=0.2
CLASSICAL_MDS ...
ARG=c1.moment-2,c1.moment-3
REWEIGHT_TEMP=0.1 TEMP=0.2
STRIDE=10
RUN=1000
NLOW_DIM=2
FMT=%8.4f
OUTPUT_FILE=list_embed
EMBEDDING_OFILE=embed
LOGWEIGHTS=ww
... CLASSICAL_MDS
src/analysis/Analysis.cpp
+ 40
− 142
View file @ 12dd3035
......@@ -57,47 +57,28 @@ P(s',t) = \frac{ \sum_{t'}^t \delta( s(x) - s' ) \exp\left( +\frac{V(x,t')}{k_B
//+ENDPLUMEDOC
void Analysis::registerKeywords( Keywords& keys ){
Action::registerKeywords( keys );
ActionPilot::registerKeywords( keys );
ActionAtomistic::registerKeywords( keys );
ActionWithArguments::registerKeywords( keys );
vesselbase::ActionWithAveraging::registerKeywords( keys );
keys.use("ARG"); keys.reset_style("ARG","optional");
keys.add("atoms","ATOMS","the atoms whose positions we are tracking for the purpose of analysing the data");
keys.add("compulsory","METRIC","EUCLIDEAN","how are we measuring the distances between configurations");
keys.add("compulsory","STRIDE","1","the frequency with which data should be stored for analysis");
keys.addFlag("USE_ALL_DATA",false,"use the data from the entire trajectory to perform the analysis");
keys.add("compulsory","RUN","the frequency with which to run the analysis algorithm. This is not required if you specify USE_ALL_DATA");
keys.add("compulsory","RUN","0","the frequency with which to run the analysis algorithm. The default value of zero assumes you want to analyse the whole trajectory");
keys.add("optional","FMT","the format that should be used in analysis output files");
keys.addFlag("REWEIGHT_BIAS",false,"reweight the data using all the biases acting on the dynamics. For more information see \\ref reweighting.");
keys.add("optional","TEMP","the system temperature. This is required if you are reweighting or doing free energies.");
keys.add("optional","REWEIGHT_TEMP","reweight data from a trajectory at one temperature and output the probability "
"distribution at a second temperature. For more information see \\ref reweighting. "
"This is not possible during postprocessing.");
keys.addFlag("WRITE_CHECKPOINT",false,"write out a checkpoint so that the analysis can be restarted in a later run");
keys.add("hidden","REUSE_DATA_FROM","eventually this will allow you to analyse the same set of data multiple times");
keys.add("hidden","IGNORE_REWEIGHTING","this allows you to ignore any reweighting factors");
keys.use("RESTART");
keys.use("UPDATE_FROM");
keys.use("UPDATE_UNTIL");
ActionWithVessel::registerKeywords( keys ); keys.remove("TOL");
keys.use("RESTART"); keys.use("UPDATE_FROM"); keys.use("UPDATE_UNTIL"); keys.remove("TOL");
}
Analysis::Analysis(const ActionOptions&ao):
Action(ao),
ActionPilot(ao),
ActionAtomistic(ao),
ActionWithArguments(ao),
ActionWithVessel(ao),
single_run(true),
ActionWithAveraging(ao),
nomemory(true),
write_chq(false),
reusing_data(false),
ignore_reweight(false),
freq(0),
needeng(false),
idata(0),
firstAnalysisDone(false),
old_norm(0.0),
//firstAnalysisDone(false),
//old_norm(0.0),
ofmt("%f"),
current_args(getNumberOfArguments()),
argument_names(getNumberOfArguments())
......@@ -132,59 +113,18 @@ argument_names(getNumberOfArguments())
if( ignore_reweight ) log.printf(" reusing data stored by %s but ignoring all reweighting\n",prev_analysis.c_str() );
else log.printf(" reusing data stored by %s\n",prev_analysis.c_str() );
} else {
if( keywords.exists("REWEIGHT_BIAS") ){
bool dobias; parseFlag("REWEIGHT_BIAS",dobias);
if( dobias ){
std::vector<ActionWithValue*> all=plumed.getActionSet().select<ActionWithValue*>();
if( all.empty() ) error("your input file is not telling plumed to calculate anything");
std::vector<Value*> arg( getArguments() );
log.printf(" reweigting using the following biases ");
for(unsigned j=0;j<all.size();j++){
std::string flab; flab=all[j]->getLabel() + ".rbias";
if( all[j]->exists(flab) ){
biases.push_back( all[j]->copyOutput(flab) );
arg.push_back( all[j]->copyOutput(flab) );
log.printf(" %s",flab.c_str());
} else {
std::string flab2; flab2=all[j]->getLabel() + ".bias";
if( all[j]->exists(flab2) ){
biases.push_back( all[j]->copyOutput(flab2) );
arg.push_back( all[j]->copyOutput(flab2) );
log.printf(" %s",flab2.c_str());
}
}
}
log.printf("\n");
if( biases.empty() ) error("you are asking to reweight bias but there does not appear to be a bias acting on your system");
requestArguments( arg );
}
}
rtemp=0;
if( keywords.exists("REWEIGHT_TEMP") ) parse("REWEIGHT_TEMP",rtemp);
if( rtemp!=0 ){
needeng=true;
log.printf(" reweighting simulation to probabilities at temperature %f\n",rtemp);
rtemp*=plumed.getAtoms().getKBoltzmann();
parse("RUN",freq);
if( freq==0 ){
log.printf(" analyzing all data in trajectory\n");
} else {
if( freq%getStride()!=0 ) error("frequncy of running is not a multiple of the stride");
log.printf(" running analysis every %u steps\n",freq);
ndata=freq/getStride(); data.resize( ndata ); logweights.resize( ndata );
for(unsigned i=0;i<ndata;++i){
data[i]=metricRegister().create<ReferenceConfiguration>( metricname );
data[i]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
}
}
simtemp=0.; parse("TEMP",simtemp);
if(simtemp>0) simtemp*=plumed.getAtoms().getKBoltzmann();
else simtemp=plumed.getAtoms().getKbT();
if( rtemp>0 || !biases.empty() ){
if(simtemp==0) error("The MD engine does not pass the temperature to plumed so you have to specify it using TEMP");
}
if( keywords.exists("USE_ALL_DATA") ){
parseFlag("USE_ALL_DATA",single_run);
if( !single_run ){
unsigned astride; parse("RUN",astride);
log.printf(" running analysis every %u steps\n",astride);
setAnalysisStride( false, astride );
} else {
log.printf(" analyzing all data in trajectory\n");
}
}
parseFlag("WRITE_CHECKPOINT",write_chq);
if( write_chq && single_run ){
write_chq=false;
......@@ -208,26 +148,12 @@ argument_names(getNumberOfArguments())
rfile.open( filename.c_str() ); // In overwrite mode automatically because there is no restart
}
if( write_chq ){
rfile.addConstantField("old_normalization");
//rfile.addConstantField("old_normalization");
for(unsigned i=0;i<getNumberOfArguments();++i) rfile.setupPrintValue( getPntrToArgument(i) );
}
}
}
void Analysis::setAnalysisStride( const bool& use_all, const unsigned& astride ){
if( freq>0 && astride!=freq ) error("each histogram can only be output with one stride");
else if( use_all ) return;
freq=astride; single_run=false;
if( astride%getStride()!= 0 ) error("Frequncy of running is not a multiple of the stride");
ndata=freq/getStride(); data.resize( ndata );
for(unsigned i=0;i<ndata;++i){
data[i]=metricRegister().create<ReferenceConfiguration>( metricname );
data[i]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
}
logweights.resize( ndata );
}
void Analysis::readDataFromFile( const std::string& filename ){
FILE* fp=fopen(filename.c_str(),"r");
if(fp!=NULL){
......@@ -243,7 +169,7 @@ void Analysis::readDataFromFile( const std::string& filename ){
error("frequency of data storage in " + filename + " is not equal to frequency of data storage plumed.dat file");
}
data[idata]->parse("LOG_WEIGHT",logweights[idata]);
data[idata]->parse("OLD_NORM",old_norm);
//data[idata]->parse("OLD_NORM",old_norm);
data[idata]->checkRead();
idata++; first=false; oldtstep=tstep;
} else{
......@@ -252,7 +178,7 @@ void Analysis::readDataFromFile( const std::string& filename ){
}
fclose(fp);
}
if(old_norm>0) firstAnalysisDone=true;
// if(old_norm>0) firstAnalysisDone=true;
}
void Analysis::parseOutputFile( const std::string& key, std::string& filename ){
......@@ -266,32 +192,11 @@ void Analysis::parseOutputFile( const std::string& key, std::string& filename ){
}
}
void Analysis::prepare(){
if(needeng) plumed.getAtoms().setCollectEnergy(true);
}
void Analysis::calculate(){
// do nothing
}
void Analysis::accumulate(){
// Don't store the first step (also don't store if we are getting data from elsewhere)
if( (!single_run && getStep()==0) || reusing_data ) return;
// This is used when we have a full quota of data from the first run
if( !single_run && idata==logweights.size() ) return;
// Retrieve the bias
double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get();
double ww=0;
if(needeng){
double energy=plumed.getAtoms().getEnergy()+bias;
// Reweighting because of temperature difference
ww=-( (1.0/rtemp) - (1.0/simtemp) )*(energy+bias);
}
// Reweighting because of biases
if( !biases.empty() ) ww += bias/ simtemp;
// Get the arguments ready to transfer to reference configuration
for(unsigned i=0;i<getNumberOfArguments();++i) current_args[i]=getArgument(i);
......@@ -300,17 +205,17 @@ void Analysis::accumulate(){
plumed_dbg_assert( data.size()==idata+1 );
data[idata]->setNamesAndAtomNumbers( getAbsoluteIndexes(), argument_names );
data[idata]->setReferenceConfig( getPositions(), current_args, getMetric() );
logweights.push_back(ww);
logweights.push_back(lweight);
} else {
// Get the arguments and store them in a vector of vectors
data[idata]->setReferenceConfig( getPositions(), current_args, getMetric() );
logweights[idata] = ww;
logweights[idata] = lweight;
}
// Write data to checkpoint file
if( write_chq ){
rfile.rewind();
data[idata]->print( rfile, getTime(), logweights[idata], atoms.getUnits().getLength()/0.1, old_norm );
data[idata]->print( rfile, getTime(), logweights[idata], atoms.getUnits().getLength()/0.1, 1.0 ); //old_norm );
rfile.flush();
}
// Increment data counter
......@@ -344,7 +249,7 @@ void Analysis::finalizeWeights( const bool& ignore_weights ){
// Check that we have the correct ammount of data
if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong. Am trying to run analysis but I don't have sufficient data");
norm=0; // Reset normalization constant
double norm=0; // Reset normalization constant
if( ignore_weights ){
for(unsigned i=0;i<logweights.size();++i){
data[i]->setWeight(1.0); norm+=1.0;
......@@ -365,18 +270,20 @@ void Analysis::finalizeWeights( const bool& ignore_weights ){
}
// Calculate normalized weights (with memory)
} else {
if( !firstAnalysisDone ) finalizeWeightsNoLogSums( 1.0 );
else finalizeWeightsNoLogSums( old_norm );
plumed_merror("analysis can now only support block averages");
// if( !firstAnalysisDone )
// finalizeWeightsNoLogSums( 1.0 );
// else finalizeWeightsNoLogSums( old_norm );
}
}
void Analysis::finalizeWeightsNoLogSums( const double& onorm ){
if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong. Am trying to run analysis but I don't have sufficient data");
// Calculate normalization constant
norm=0; for(unsigned i=0;i<logweights.size();++i) norm+=exp( logweights[i] );
// Calculate weights (with memory)
for(unsigned i=0;i<logweights.size();++i) data[i]->setWeight( exp( logweights[i] ) / onorm );
}
// void Analysis::finalizeWeightsNoLogSums( const double& onorm ){
// if( !reusing_data && idata!=logweights.size() ) error("something has gone wrong. Am trying to run analysis but I don't have sufficient data");
// // Calculate normalization constant
// double norm=0; for(unsigned i=0;i<logweights.size();++i) norm+=exp( logweights[i] );
// // Calculate weights (with memory)
// for(unsigned i=0;i<logweights.size();++i) data[i]->setWeight( exp( logweights[i] ) / norm );
// }
void Analysis::getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const {
plumed_dbg_assert( getNumberOfAtoms()==0 );
......@@ -399,27 +306,18 @@ void Analysis::runAnalysis(){
finalizeWeights( ignore_reweight );
} else {
mydatastash->finalizeWeights( ignore_reweight );
norm=mydatastash->retrieveNorm();
// norm=mydatastash->retrieveNorm();
}
// This ensures everything is set up to run the calculation
if( single_run ) setAnalysisStride( single_run, freq );
// if( single_run ) setAnalysisStride( single_run, freq );
// And run the analysis
performAnalysis(); idata=0;
// Update total normalization constant
old_norm+=norm; firstAnalysisDone=true;
}
double Analysis::getNormalization() const {
if( nomemory || !firstAnalysisDone ) return norm;
return ( 1. + norm/old_norm );
}
// old_norm+=norm; firstAnalysisDone=true;
double Analysis::getTemp() const {
return simtemp;
}
void Analysis::update(){
void Analysis::performOperations( const bool& from_update ){
accumulate();
if( !single_run ){
if( getStep()>0 && getStep()%freq==0 ) runAnalysis();
......@@ -434,7 +332,7 @@ bool Analysis::getPeriodicityInformation(const unsigned& i, std::string& dmin, s
}
void Analysis::runFinalJobs() {
if( !single_run ) return;
if( freq>0 ) return;
if( getNumberOfDataPoints()==0 ) error("no data is available for analysis");
runAnalysis();
}
......
src/analysis/Analysis.h
+ 6
− 68
View file @ 12dd3035
......@@ -22,10 +22,7 @@
#ifndef __PLUMED_analysis_Analysis_h
#define __PLUMED_analysis_Analysis_h
#include "core/ActionPilot.h"
#include "core/ActionAtomistic.h"
#include "core/ActionWithArguments.h"
#include "vesselbase/ActionWithVessel.h"
#include "vesselbase/ActionWithAveraging.h"
#define PLUMED_ANALYSIS_INIT(ao) Action(ao),Analysis(ao)
......@@ -42,12 +39,7 @@ is information as to how to go about implementing a new analysis method.
*/
class Analysis :
public ActionPilot,
public ActionAtomistic,
public ActionWithArguments,
public vesselbase::ActionWithVessel
{
class Analysis : public vesselbase::ActionWithAveraging {
private:
/// Are we running only once for the whole trajectory
bool single_run;
......@@ -71,16 +63,14 @@ private:
double rtemp;
/// Do we need the energy (are we reweighting at a different temperature)
bool needeng;
/// The biases we are using in reweighting and the args we store them separately
std::vector<Value*> biases;
/// The piece of data we are inserting
unsigned idata;
/// The weights of all the data points
std::vector<double> logweights;
/// Have we analyzed the data for the first time
bool firstAnalysisDone;
// bool firstAnalysisDone;
/// The value of the old normalization constant
double norm, old_norm;
// double norm, old_norm;
/// The format to use in output files
std::string ofmt;
/// Tempory vector to store values of arguments
......@@ -91,11 +81,6 @@ private:
OFile rfile;
/// Read in data from a file
void readDataFromFile( const std::string& filename );
/// This retrieves the value of norm from the analysis action.
/// It is only used to transfer data from one analysis action to
/// another. You should never need to use it. If you think you need it
/// you probably need getNormalization()
double retrieveNorm() const ;
/// Get the metric if we are using malonobius distance and flexible hill
std::vector<double> getMetric() const ;
protected:
......@@ -109,8 +94,6 @@ protected:
std::vector<ReferenceConfiguration*> data;
/// Get the name of the metric we are using to measure distances
std::string getMetricName() const ;
/// Return the number of arguments (this overwrites the one in ActionWithArguments)
unsigned getNumberOfArguments() const;
/// Return the number of data points
unsigned getNumberOfDataPoints() const;
/// Return the weight of the ith point
......@@ -119,37 +102,22 @@ protected:
void getDataPoint( const unsigned& idata, std::vector<double>& point, double& weight ) const ;
/// Returns true if argument i is periodic together with the domain
bool getPeriodicityInformation(const unsigned& i, std::string& dmin, std::string& dmax);
/// Return the normalization constant
double getNormalization() const;
/// Return the set temperature (N.B. k_B T is what is returned by this function)
double getTemp () const;
/// Are we analyzing each data block separately (if we are not this also returns the old normalization )
bool usingMemory() const;
/// Overwrite ActionWithArguments getArguments() so that we don't return
/// the bias
std::vector<Value*> getArguments();
/// Return the format to use for numbers in output files
std::string getOutputFormat() const ;
/// Finalize the weights without using the log sums
void finalizeWeightsNoLogSums( const double& onorm );
// void finalizeWeightsNoLogSums( const double& onorm );
public:
static void registerKeywords( Keywords& keys );
explicit Analysis(const ActionOptions&);
~Analysis();
void prepare();
void calculate();
void update();
void accumulate();
void performOperations( const bool& from_update );
virtual void performAnalysis()=0;
void apply(){}
void runFinalJobs();
void runAnalysis();
void lockRequests();
void unlockRequests();
void calculateNumericalDerivatives( ActionWithValue* a=NULL ){ plumed_error(); }
bool isPeriodic(){ plumed_error(); return false; }
unsigned getNumberOfDerivatives(){ return 0; }
virtual void setAnalysisStride( const bool& use_all, const unsigned& astride );
/// Convert the stored log weights to proper weights
virtual void finalizeWeights( const bool& ignore_weights );
};
......@@ -159,18 +127,6 @@ std::string Analysis::getMetricName() const {
return metricname;
}
inline
void Analysis::lockRequests(){
ActionAtomistic::lockRequests();
ActionWithArguments::lockRequests();
}
inline
void Analysis::unlockRequests(){
ActionAtomistic::unlockRequests();
ActionWithArguments::unlockRequests();
}
inline
unsigned Analysis::getNumberOfDataPoints() const {
if( !reusing_data ){
......@@ -190,24 +146,6 @@ bool Analysis::usingMemory() const {
}
}
inline
unsigned Analysis::getNumberOfArguments() const {
unsigned nargs=ActionWithArguments::getNumberOfArguments();
return nargs - biases.size();
}
inline
double Analysis::retrieveNorm() const {
return norm;
}
inline
std::vector<Value*> Analysis::getArguments(){
std::vector<Value*> arg_vals( ActionWithArguments::getArguments() );
for(unsigned i=0;i<biases.size();++i) arg_vals.erase(arg_vals.end()-1);
return arg_vals;
}
inline
std::string Analysis::getOutputFormat() const {
return ofmt;
......
src/bias/ReweightTemperature.cpp
+ 5
− 1
View file @ 12dd3035
......@@ -44,6 +44,7 @@ private:
public:
static void registerKeywords(Keywords&);
explicit ReweightTemperature(const ActionOptions&ao);
void prepare();
double getLogWeight() const ;
};
......@@ -59,7 +60,6 @@ ReweightTemperature::ReweightTemperature(const ActionOptions&ao):
Action(ao),
ReweightBase(ao)
{
plumed.getAtoms().setCollectEnergy(true);
parse("REWEIGHT_TEMP",rtemp);
log.printf(" reweighting simulation to probabilities at temperature %f\n",rtemp);
rtemp*=plumed.getAtoms().getKBoltzmann();
......@@ -67,6 +67,10 @@ ReweightBase(ao)
retrieveAllBiases( "bias", biases );
}
void ReweightTemperature::prepare(){
plumed.getAtoms().setCollectEnergy(true);
}
double ReweightTemperature::getLogWeight() const {
// Retrieve the bias
double bias=0.0; for(unsigned i=0;i<biases.size();++i) bias+=biases[i]->get();
......
src/vesselbase/ActionWithAveraging.cpp
+ 3
− 2
View file @ 12dd3035
......@@ -43,6 +43,7 @@ ActionAtomistic(ao),
ActionWithArguments(ao),
ActionWithValue(ao),
ActionWithVessel(ao),
lweight(0),cweight(0),
myaverage(NULL),
useRunAllTasks(false),
clearstride(0)
......@@ -100,9 +101,9 @@ void ActionWithAveraging::update(){
// Calculate the weight for all reweighting
if ( weights.size()>0 ){
double sum=0; for(unsigned i=0;i<weights.size();++i) sum+=weights[i]->get();
cweight = exp( sum );
lweight=sum; cweight = exp( sum );
} else {
cweight=1.0;
lweight=0; cweight=1.0;
}
// Prepare to do the averaging
prepareForAveraging();
......
src/vesselbase/ActionWithAveraging.h
+ 11
− 2
View file @ 12dd3035
......@@ -62,8 +62,8 @@ private:
/// Are we accumulated the unormalized quantity
bool unormalised;
protected:
/// The current weight
double cweight;
/// The current weight and its logarithm
double lweight, cweight;
/// Set the averaging action
void setAveragingAction( AveragingVessel* av_vessel, const bool& usetasks );
public:
......@@ -74,6 +74,8 @@ public:
void calculateNumericalDerivatives(PLMD::ActionWithValue*);
virtual unsigned getNumberOfDerivatives(){ return 0; }
unsigned getNumberOfArguments() const ;
/// Overwrite ActionWithArguments getArguments() so that we don't return the bias
std::vector<Value*> getArguments();
void calculate(){}
void apply(){}
void update();
......@@ -92,6 +94,13 @@ unsigned ActionWithAveraging::getNumberOfArguments() const {
return ActionWithArguments::getNumberOfArguments() - weights.size();
}
inline
std::vector<Value*> ActionWithAveraging::getArguments(){
std::vector<Value*> arg_vals( ActionWithArguments::getArguments() );
for(unsigned i=0;i<weights.size();++i) arg_vals.erase(arg_vals.end()-1);
return arg_vals;
}
}
}
#endif
src/vesselbase/ActionWithVessel.h
+ 0
− 2
View file @ 12dd3035
......@@ -188,8 +188,6 @@ public:
Vessel* getVesselWithName( const std::string& mynam );
/// Does the weight have derivatives
bool weightWithDerivatives() const ;
/// This is used to set the stride for the output of histograms
virtual void setAnalysisStride( const bool& use_all, const unsigned& astride ){}
/// Return the position in the current task list
unsigned getPositionInCurrentTaskList( const unsigned& myind ) const ;
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment