Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "bias/Bias.h"
24 : #include "core/ActionRegister.h"
25 : #include "MetainferenceCore.h"
26 :
27 : namespace PLMD {
28 : namespace isdb {
29 :
30 : //+PLUMEDOC ISDB_BIAS METAINFERENCE
31 : /*
32 : Calculates the Metainference energy for a set of experimental data.
33 :
34 : The metainference method that is introduced in the first paper cited below is a Bayesian framework
35 : to model heterogeneous systems by integrating prior information with noisy, ensemble-averaged data.
36 : Metainference models a system and quantifies the level of noise in the data by considering a set of replicas of the system.
37 :
38 : Calculated experimental data are given in input as ARG while reference experimental values
39 : can be given either from fixed components of other actions using PARARG or as numbers using
40 : PARAMETERS. The default behavior is that of averaging the data over the available replicas,
41 : if this is not wanted the keyword NOENSEMBLE prevent this averaging.
42 :
43 : Metadynamics Metainference (see second paper cited below) more in general biased Metainference requires the knowledge of
44 : biasing potential in order to calculate the weighted average. In this case the value of the bias
45 : can be provided as the last argument in ARG and adding the keyword REWEIGHT. To avoid the noise
46 : resulting from the instantaneous value of the bias the weight of each replica can be averaged
47 : over a give time using the keyword AVERAGING.
48 :
49 : The data can be averaged by using multiple replicas and weighted for a bias if present.
50 : The functional form of Metainference can be chosen among four variants selected
51 : with NOISE=GAUSS,MGAUSS,OUTLIERS,MOUTLIERS,GENERIC which correspond to modelling the noise for
52 : the arguments as a single gaussian common to all the data points, a gaussian per data
53 : point, a single long-tailed gaussian common to all the data points, a log-tailed
54 : gaussian per data point or using two distinct noises as for the most general formulation of Metainference.
55 : In this latter case the noise of the replica-averaging is gaussian (one per data point) and the noise for
56 : the comparison with the experimental data can chosen using the keyword LIKELIHOOD
57 : between gaussian or log-normal (one per data point), furthermore the evolution of the estimated average
58 : over an infinite number of replicas is driven by DFTILDE.
59 :
60 : As for Metainference theory there are two sigma values: SIGMA_MEAN0 represent the
61 : error of calculating an average quantity using a finite set of replica and should
62 : be set as small as possible following the guidelines for replica-averaged simulations
63 : in the framework of the Maximum Entropy Principle. Alternatively, this can be obtained
64 : automatically using the internal sigma mean optimization as introduced in the third paper cited below
65 : (OPTSIGMAMEAN=SEM), in this second case sigma_mean is estimated from the maximum standard error
66 : of the mean either over the simulation or over a defined time using the keyword AVERAGING.
67 : SIGMA_BIAS is an uncertainty parameter, sampled by a MC algorithm in the bounded interval
68 : defined by SIGMA_MIN and SIGMA_MAX. The initial value is set at SIGMA0. The MC move is a
69 : random displacement of maximum value equal to DSIGMA. If the number of data point is
70 : too large and the acceptance rate drops it is possible to make the MC move over mutually
71 : exclusive, random subset of size MC_CHUNKSIZE and run more than one move setting MC_STEPS
72 : in such a way that MC_CHUNKSIZE*MC_STEPS will cover all the data points.
73 :
74 : Calculated and experimental data can be compared modulo a scaling factor and/or an offset
75 : using SCALEDATA and/or ADDOFFSET, the sampling is obtained by a MC algorithm either using
76 : a flat or a gaussian prior setting it with SCALE_PRIOR or OFFSET_PRIOR. With a flat prior the
77 : moves are reflected into the interval given by SCALE_MIN and SCALE_MAX (OFFSET_MIN and OFFSET_MAX),
78 : with a gaussian prior the prior is centered in SCALE0 (OFFSET0) and its width is DSCALE (DOFFSET),
79 : which also sets the size of the moves.
80 :
81 : The random numbers of the Monte Carlo are seeded from the clock, so that two runs are different;
82 : SEED makes them reproducible. The status of Metainference (uncertainties, history of sigma mean,
83 : scale, offset, averaged weights) is written every WRITE_STRIDE steps in the file MISTATUS followed
84 : by the label of the action (STATUS_FILE changes MISTATUS) and is read back when restarting.
85 :
86 : The same keywords can be used in the ISDB collective variables (e.g. [NOE](NOE.md), [RDC](RDC.md))
87 : together with DOSCORE: in this case the score is calculated directly by the collective variable,
88 : which is faster with many data.
89 :
90 : ## Examples
91 :
92 : In the following example we calculate a set of [RDC](RDC.md), take the replica-average of
93 : them and comparing them with a set of experimental values. RDCs are compared with
94 : the experimental data but for a multiplication factor SCALE that is also sampled by
95 : MC on-the-fly
96 :
97 : ```plumed
98 : RDC ...
99 : LABEL=rdc
100 : SCALE=0.0001
101 : GYROM=-72.5388
102 : ATOMS1=22,23
103 : ATOMS2=25,27
104 : ATOMS3=29,31
105 : ATOMS4=33,34
106 : ... RDC
107 :
108 : METAINFERENCE ...
109 : ARG=rdc.*
110 : NOISETYPE=MGAUSS
111 : PARAMETERS=1.9190,2.9190,3.9190,4.9190
112 : SCALEDATA SCALE0=1 SCALE_MIN=0.1 SCALE_MAX=3 DSCALE=0.01
113 : SIGMA0=0.01 SIGMA_MIN=0.00001 SIGMA_MAX=3 DSIGMA=0.01
114 : SIGMA_MEAN0=0.001
115 : LABEL=spe
116 : ... METAINFERENCE
117 :
118 : PRINT ARG=spe.bias FILE=BIAS STRIDE=1
119 : ```
120 :
121 : in the following example instead of using one uncertainty parameter per data point we use
122 : a single uncertainty value in a long-tailed gaussian to take into account for outliers, furthermore
123 : the data are weighted for the bias applied to other variables of the system.
124 :
125 : ```plumed
126 : RDC ...
127 : LABEL=rdc
128 : SCALE=0.0001
129 : GYROM=-72.5388
130 : ATOMS1=22,23
131 : ATOMS2=25,27
132 : ATOMS3=29,31
133 : ATOMS4=33,34
134 : ... RDC
135 :
136 : cv1: TORSION ATOMS=1,2,3,4
137 : cv2: TORSION ATOMS=2,3,4,5
138 : mm: METAD ARG=cv1,cv2 HEIGHT=0.5 SIGMA=0.3,0.3 PACE=200 BIASFACTOR=8 WALKERS_MPI
139 :
140 : METAINFERENCE ...
141 : #SETTINGS NREPLICAS=2
142 : ARG=rdc.*,mm.bias
143 : REWEIGHT
144 : NOISETYPE=OUTLIERS
145 : PARAMETERS=1.9190,2.9190,3.9190,4.9190
146 : SCALEDATA SCALE0=1 SCALE_MIN=0.1 SCALE_MAX=3 DSCALE=0.01
147 : SIGMA0=0.01 SIGMA_MIN=0.00001 SIGMA_MAX=3 DSIGMA=0.01
148 : SIGMA_MEAN0=0.001
149 : LABEL=spe
150 : ... METAINFERENCE
151 : ```
152 :
153 : (See also [RDC](RDC.md), [PBMETAD](PBMETAD.md)).
154 :
155 : */
156 : //+ENDPLUMEDOC
157 :
158 : class Metainference : public bias::Bias {
159 : // number of experimental data (the bias used with REWEIGHT is not counted)
160 : unsigned narg;
161 : // the Metainference algorithm
162 : MetainferenceCore core_;
163 :
164 : public:
165 : explicit Metainference(const ActionOptions&);
166 : void calculate() override;
167 : void update() override;
168 : static void registerKeywords(Keywords& keys);
169 : };
170 :
171 :
172 : PLUMED_REGISTER_ACTION(Metainference,"METAINFERENCE")
173 :
174 93 : void Metainference::registerKeywords(Keywords& keys) {
175 93 : Bias::registerKeywords(keys);
176 186 : keys.addInputKeyword("optional","PARARG","scalar","reference values for the experimental data, these can be provided as arguments without derivatives");
177 93 : keys.add("optional","PARAMETERS","reference values for the experimental data");
178 93 : MetainferenceCore::registerKeywords(keys);
179 93 : keys.addDOI("10.1126/sciadv.1501177");
180 93 : keys.addDOI("10.1038/srep31232");
181 93 : keys.addDOI("10.1063/1.4981211");
182 93 : }
183 :
184 89 : Metainference::Metainference(const ActionOptions&ao):
185 : PLUMED_BIAS_INIT(ao),
186 89 : narg(0),
187 89 : core_(*this,true) {
188 89 : narg = getNumberOfArguments();
189 89 : if(core_.doReweight()) {
190 16 : narg--;
191 : }
192 :
193 89 : std::vector<double> parameters;
194 178 : parseVector("PARAMETERS",parameters);
195 89 : if(parameters.size()!=static_cast<unsigned>(narg)&&!parameters.empty()) {
196 0 : error("Size of PARAMETERS array should be either 0 or the same as of the number of arguments in ARG1");
197 : }
198 :
199 89 : std::vector<Value*> arg2;
200 178 : parseArgumentList("PARARG",arg2);
201 89 : if(!arg2.empty()) {
202 4 : if(parameters.size()>0) {
203 0 : error("It is not possible to use PARARG and PARAMETERS together");
204 : }
205 4 : if(arg2.size()!=narg) {
206 0 : error("Size of PARARG array should be the same as number for arguments in ARG");
207 : }
208 2360 : for(unsigned i=0; i<arg2.size(); i++) {
209 2356 : parameters.push_back(arg2[i]->get());
210 2356 : if(arg2[i]->hasDerivatives()==true) {
211 0 : error("PARARG can only accept arguments without derivatives");
212 : }
213 : }
214 : }
215 :
216 89 : if(parameters.size()!=narg) {
217 0 : error("PARARG or PARAMETERS arrays should include the same number of elements as the arguments in ARG");
218 : }
219 :
220 89 : checkRead();
221 :
222 178 : core_.setup(parameters, *this, [this](const std::string& name) {
223 5942 : return getPntrToComponent(name);
224 : });
225 89 : }
226 :
227 512 : void Metainference::calculate() {
228 : std::vector<double>& data=core_.data();
229 5702 : for(unsigned i=0; i<narg; ++i) {
230 5190 : data[i]=getArgument(i);
231 : }
232 :
233 512 : const double ene = core_.calculate(core_.getSelector(), core_.doReweight() ? getArgument(narg) : 0.);
234 :
235 5702 : for(unsigned i=0; i<narg; ++i) {
236 5190 : setOutputForce(i, -core_.getMetaDer(i));
237 : }
238 512 : if(core_.doReweight()) {
239 148 : setOutputForce(narg, -core_.getBiasDer());
240 : }
241 :
242 512 : setBias(ene);
243 512 : }
244 :
245 512 : void Metainference::update() {
246 512 : core_.update();
247 512 : }
248 :
249 : }
250 : }
251 :
|
