Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2020 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 "bias/ActionRegister.h"
25 : #include "core/PlumedMain.h"
26 : #include "core/Atoms.h"
27 : #include "core/Value.h"
28 : #include "tools/File.h"
29 : #include "tools/OpenMP.h"
30 : #include "tools/Random.h"
31 : #include <cmath>
32 : #include <ctime>
33 : #include <numeric>
34 : using namespace std;
35 :
36 : #ifndef M_PI
37 : #define M_PI 3.14159265358979323846
38 : #endif
39 :
40 : namespace PLMD {
41 : namespace isdb {
42 :
43 : //+PLUMEDOC ISDB_BIAS METAINFERENCE
44 : /*
45 : Calculates the Metainference energy for a set of experimental data.
46 :
47 : Metainference \cite Bonomi:2016ip is a Bayesian framework
48 : to model heterogeneous systems by integrating prior information with noisy, ensemble-averaged data.
49 : Metainference models a system and quantifies the level of noise in the data by considering a set of replicas of the system.
50 :
51 : Calculated experimental data are given in input as ARG while reference experimental values
52 : can be given either from fixed components of other actions using PARARG or as numbers using
53 : PARAMETERS. The default behavior is that of averaging the data over the available replicas,
54 : if this is not wanted the keyword NOENSEMBLE prevent this averaging.
55 :
56 : Metadynamics Metainference \cite Bonomi:2016ge or more in general biased Metainference requires the knowledge of
57 : biasing potential in order to calculate the weighted average. In this case the value of the bias
58 : can be provided as the last argument in ARG and adding the keyword REWEIGHT. To avoid the noise
59 : resulting from the instantaneous value of the bias the weight of each replica can be averaged
60 : over a give time using the keyword AVERAGING.
61 :
62 : The data can be averaged by using multiple replicas and weighted for a bias if present.
63 : The functional form of Metainference can be chosen among four variants selected
64 : with NOISE=GAUSS,MGAUSS,OUTLIERS,MOUTLIERS,GENERIC which correspond to modelling the noise for
65 : the arguments as a single gaussian common to all the data points, a gaussian per data
66 : point, a single long-tailed gaussian common to all the data points, a log-tailed
67 : gaussian per data point or using two distinct noises as for the most general formulation of Metainference.
68 : In this latter case the noise of the replica-averaging is gaussian (one per data point) and the noise for
69 : the comparison with the experimental data can chosen using the keyword LIKELIHOOD
70 : between gaussian or log-normal (one per data point), furthermore the evolution of the estimated average
71 : over an infinite number of replicas is driven by DFTILDE.
72 :
73 : As for Metainference theory there are two sigma values: SIGMA_MEAN represent the
74 : error of calculating an average quantity using a finite set of replica and should
75 : be set as small as possible following the guidelines for replica-averaged simulations
76 : in the framework of the Maximum Entropy Principle. Alternatively, this can be obtained
77 : automatically using the internal sigma mean optimization as introduced in \cite Lohr:2017gc
78 : (OPTSIGMAMEAN=SEM), in this second case sigma_mean is estimated from the maximum standard error
79 : of the mean either over the simulation or over a defined time using the keyword AVERAGING.
80 : SIGMA_BIAS is an uncertainty parameter, sampled by a MC algorithm in the bounded interval
81 : defined by SIGMA_MIN and SIGMA_MAX. The initial value is set at SIGMA0. The MC move is a
82 : random displacement of maximum value equal to DSIGMA. If the number of data point is
83 : too large and the acceptance rate drops it is possible to make the MC move over mutually
84 : exclusive, random subset of size MC_CHUNKSIZE and run more than one move setting MC_STEPS
85 : in such a way that MC_CHUNKSIZE*MC_STEPS will cover all the data points.
86 :
87 : Calculated and experimental data can be compared modulo a scaling factor and/or an offset
88 : using SCALEDATA and/or ADDOFFSET, the sampling is obtained by a MC algorithm either using
89 : a flat or a gaussian prior setting it with SCALE_PRIOR or OFFSET_PRIOR.
90 :
91 : \par Examples
92 :
93 : In the following example we calculate a set of \ref RDC, take the replica-average of
94 : them and comparing them with a set of experimental values. RDCs are compared with
95 : the experimental data but for a multiplication factor SCALE that is also sampled by
96 : MC on-the-fly
97 :
98 : \plumedfile
99 : RDC ...
100 : LABEL=rdc
101 : SCALE=0.0001
102 : GYROM=-72.5388
103 : ATOMS1=22,23
104 : ATOMS2=25,27
105 : ATOMS3=29,31
106 : ATOMS4=33,34
107 : ... RDC
108 :
109 : METAINFERENCE ...
110 : ARG=rdc.*
111 : NOISETYPE=MGAUSS
112 : PARAMETERS=1.9190,2.9190,3.9190,4.9190
113 : SCALEDATA SCALE0=1 SCALE_MIN=0.1 SCALE_MAX=3 DSCALE=0.01
114 : SIGMA0=0.01 SIGMA_MIN=0.00001 SIGMA_MAX=3 DSIGMA=0.01
115 : SIGMA_MEAN0=0.001
116 : LABEL=spe
117 : ... METAINFERENCE
118 :
119 : PRINT ARG=spe.bias FILE=BIAS STRIDE=1
120 : \endplumedfile
121 :
122 : in the following example instead of using one uncertainty parameter per data point we use
123 : a single uncertainty value in a long-tailed gaussian to take into account for outliers, furthermore
124 : the data are weighted for the bias applied to other variables of the system.
125 :
126 : \plumedfile
127 : RDC ...
128 : LABEL=rdc
129 : SCALE=0.0001
130 : GYROM=-72.5388
131 : ATOMS1=22,23
132 : ATOMS2=25,27
133 : ATOMS3=29,31
134 : ATOMS4=33,34
135 : ... RDC
136 :
137 : cv1: TORSION ATOMS=1,2,3,4
138 : cv2: TORSION ATOMS=2,3,4,5
139 : mm: METAD ARG=cv1,cv2 HEIGHT=0.5 SIGMA=0.3,0.3 PACE=200 BIASFACTOR=8 WALKERS_MPI
140 :
141 : METAINFERENCE ...
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_MEAN=0.001
149 : LABEL=spe
150 : ... METAINFERENCE
151 : \endplumedfile
152 :
153 : (See also \ref RDC, \ref PBMETAD).
154 :
155 : */
156 : //+ENDPLUMEDOC
157 :
158 : class Metainference : public bias::Bias
159 : {
160 : // experimental values
161 : vector<double> parameters;
162 : // noise type
163 : unsigned noise_type_;
164 : enum { GAUSS, MGAUSS, OUTLIERS, MOUTLIERS, GENERIC };
165 : unsigned gen_likelihood_;
166 : enum { LIKE_GAUSS, LIKE_LOGN };
167 : // scale is data scaling factor
168 : // noise type
169 : unsigned scale_prior_;
170 : enum { SC_GAUSS, SC_FLAT };
171 : bool doscale_;
172 : double scale_;
173 : double scale_mu_;
174 : double scale_min_;
175 : double scale_max_;
176 : double Dscale_;
177 : // scale is data scaling factor
178 : // noise type
179 : unsigned offset_prior_;
180 : bool dooffset_;
181 : double offset_;
182 : double offset_mu_;
183 : double offset_min_;
184 : double offset_max_;
185 : double Doffset_;
186 : // scale and offset regression
187 : bool doregres_zero_;
188 : int nregres_zero_;
189 : // sigma is data uncertainty
190 : vector<double> sigma_;
191 : vector<double> sigma_min_;
192 : vector<double> sigma_max_;
193 : vector<double> Dsigma_;
194 : // sigma_mean is uncertainty in the mean estimate
195 : vector<double> sigma_mean2_;
196 : // this is the estimator of the mean value per replica for generic metainference
197 : vector<double> ftilde_;
198 : double Dftilde_;
199 :
200 : // temperature in kbt
201 : double kbt_;
202 :
203 : // Monte Carlo stuff
204 : vector<Random> random;
205 : unsigned MCsteps_;
206 : unsigned MCstride_;
207 : long unsigned MCaccept_;
208 : long unsigned MCacceptScale_;
209 : long unsigned MCacceptFT_;
210 : long unsigned MCtrial_;
211 : unsigned MCchunksize_;
212 :
213 : // output
214 : Value* valueScale;
215 : Value* valueOffset;
216 : Value* valueAccept;
217 : Value* valueAcceptScale;
218 : Value* valueAcceptFT;
219 : vector<Value*> valueSigma;
220 : vector<Value*> valueSigmaMean;
221 : vector<Value*> valueFtilde;
222 :
223 : // restart
224 : unsigned write_stride_;
225 : OFile sfile_;
226 :
227 : // others
228 : bool firstTime;
229 : vector<bool> firstTimeW;
230 : bool master;
231 : bool do_reweight_;
232 : unsigned do_optsigmamean_;
233 : unsigned nrep_;
234 : unsigned replica_;
235 : unsigned narg;
236 :
237 : // selector
238 : string selector_;
239 :
240 : // optimize sigma mean
241 : vector< vector < vector <double> > > sigma_mean2_last_;
242 : unsigned optsigmamean_stride_;
243 :
244 : // average weights
245 : unsigned average_weights_stride_;
246 : vector< vector <double> > average_weights_;
247 :
248 : double getEnergyMIGEN(const vector<double> &mean, const vector<double> &ftilde, const vector<double> &sigma,
249 : const double scale, const double offset);
250 : double getEnergySP(const vector<double> &mean, const vector<double> &sigma,
251 : const double scale, const double offset);
252 : double getEnergySPE(const vector<double> &mean, const vector<double> &sigma,
253 : const double scale, const double offset);
254 : double getEnergyGJ(const vector<double> &mean, const vector<double> &sigma,
255 : const double scale, const double offset);
256 : double getEnergyGJE(const vector<double> &mean, const vector<double> &sigma,
257 : const double scale, const double offset);
258 : void doMonteCarlo(const vector<double> &mean);
259 : double getEnergyForceMIGEN(const vector<double> &mean, const vector<double> &dmean_x, const vector<double> &dmean_b);
260 : double getEnergyForceSP(const vector<double> &mean, const vector<double> &dmean_x, const vector<double> &dmean_b);
261 : double getEnergyForceSPE(const vector<double> &mean, const vector<double> &dmean_x, const vector<double> &dmean_b);
262 : double getEnergyForceGJ(const vector<double> &mean, const vector<double> &dmean_x, const vector<double> &dmean_b);
263 : double getEnergyForceGJE(const vector<double> &mean, const vector<double> &dmean_x, const vector<double> &dmean_b);
264 : void get_weights(const unsigned iselect, double &fact, double &var_fact);
265 : void replica_averaging(const double fact, std::vector<double> &mean, std::vector<double> &dmean_b);
266 : void get_sigma_mean(const unsigned iselect, const double fact, const double var_fact, const vector<double> &mean);
267 : void writeStatus();
268 : void do_regression_zero(const vector<double> &mean);
269 :
270 : public:
271 : explicit Metainference(const ActionOptions&);
272 : ~Metainference();
273 : void calculate();
274 : void update();
275 : static void registerKeywords(Keywords& keys);
276 : };
277 :
278 :
279 7394 : PLUMED_REGISTER_ACTION(Metainference,"METAINFERENCE")
280 :
281 20 : void Metainference::registerKeywords(Keywords& keys) {
282 20 : Bias::registerKeywords(keys);
283 40 : keys.use("ARG");
284 80 : keys.add("optional","PARARG","reference values for the experimental data, these can be provided as arguments without derivatives");
285 80 : keys.add("optional","PARAMETERS","reference values for the experimental data");
286 60 : keys.addFlag("NOENSEMBLE",false,"don't perform any replica-averaging");
287 60 : keys.addFlag("REWEIGHT",false,"simple REWEIGHT using the latest ARG as energy");
288 80 : keys.add("optional","AVERAGING", "Stride for calculation of averaged weights and sigma_mean");
289 100 : keys.add("compulsory","NOISETYPE","MGAUSS","functional form of the noise (GAUSS,MGAUSS,OUTLIERS,MOUTLIERS,GENERIC)");
290 100 : keys.add("compulsory","LIKELIHOOD","GAUSS","the likelihood for the GENERIC metainference model, GAUSS or LOGN");
291 100 : keys.add("compulsory","DFTILDE","0.1","fraction of sigma_mean used to evolve ftilde");
292 60 : keys.addFlag("SCALEDATA",false,"Set to TRUE if you want to sample a scaling factor common to all values and replicas");
293 100 : keys.add("compulsory","SCALE0","1.0","initial value of the scaling factor");
294 100 : keys.add("compulsory","SCALE_PRIOR","FLAT","either FLAT or GAUSSIAN");
295 80 : keys.add("optional","SCALE_MIN","minimum value of the scaling factor");
296 80 : keys.add("optional","SCALE_MAX","maximum value of the scaling factor");
297 80 : keys.add("optional","DSCALE","maximum MC move of the scaling factor");
298 60 : keys.addFlag("ADDOFFSET",false,"Set to TRUE if you want to sample an offset common to all values and replicas");
299 100 : keys.add("compulsory","OFFSET0","0.0","initial value of the offset");
300 100 : keys.add("compulsory","OFFSET_PRIOR","FLAT","either FLAT or GAUSSIAN");
301 80 : keys.add("optional","OFFSET_MIN","minimum value of the offset");
302 80 : keys.add("optional","OFFSET_MAX","maximum value of the offset");
303 80 : keys.add("optional","DOFFSET","maximum MC move of the offset");
304 80 : keys.add("optional","REGRES_ZERO","stride for regression with zero offset");
305 100 : keys.add("compulsory","SIGMA0","1.0","initial value of the uncertainty parameter");
306 100 : keys.add("compulsory","SIGMA_MIN","0.0","minimum value of the uncertainty parameter");
307 100 : keys.add("compulsory","SIGMA_MAX","10.","maximum value of the uncertainty parameter");
308 80 : keys.add("optional","DSIGMA","maximum MC move of the uncertainty parameter");
309 100 : keys.add("compulsory","OPTSIGMAMEAN","NONE","Set to NONE/SEM to manually set sigma mean, or to estimate it on the fly");
310 80 : keys.add("optional","SIGMA_MEAN0","starting value for the uncertainty in the mean estimate");
311 80 : keys.add("optional","TEMP","the system temperature - this is only needed if code doesn't pass the temperature to plumed");
312 80 : keys.add("optional","MC_STEPS","number of MC steps");
313 80 : keys.add("optional","MC_STRIDE","MC stride");
314 80 : keys.add("optional","MC_CHUNKSIZE","MC chunksize");
315 80 : keys.add("optional","STATUS_FILE","write a file with all the data useful for restart/continuation of Metainference");
316 100 : keys.add("compulsory","WRITE_STRIDE","10000","write the status to a file every N steps, this can be used for restart/continuation");
317 80 : keys.add("optional","SELECTOR","name of selector");
318 80 : keys.add("optional","NSELECT","range of values for selector [0, N-1]");
319 40 : keys.use("RESTART");
320 20 : useCustomisableComponents(keys);
321 80 : keys.addOutputComponent("sigma", "default", "uncertainty parameter");
322 80 : keys.addOutputComponent("sigmaMean", "default", "uncertainty in the mean estimate");
323 80 : keys.addOutputComponent("acceptSigma", "default", "MC acceptance");
324 80 : keys.addOutputComponent("acceptScale", "SCALEDATA", "MC acceptance");
325 80 : keys.addOutputComponent("weight", "REWEIGHT", "weights of the weighted average");
326 80 : keys.addOutputComponent("biasDer", "REWEIGHT", "derivatives with respect to the bias");
327 80 : keys.addOutputComponent("scale", "SCALEDATA", "scale parameter");
328 80 : keys.addOutputComponent("offset", "ADDOFFSET", "offset parameter");
329 80 : keys.addOutputComponent("ftilde", "GENERIC", "ensemble average estimator");
330 20 : }
331 :
332 19 : Metainference::Metainference(const ActionOptions&ao):
333 : PLUMED_BIAS_INIT(ao),
334 : doscale_(false),
335 : scale_(1.),
336 : scale_mu_(0),
337 : scale_min_(1),
338 : scale_max_(-1),
339 : Dscale_(-1),
340 : dooffset_(false),
341 : offset_(0.),
342 : offset_mu_(0),
343 : offset_min_(1),
344 : offset_max_(-1),
345 : Doffset_(-1),
346 : doregres_zero_(false),
347 : nregres_zero_(0),
348 : Dftilde_(0.1),
349 : random(3),
350 : MCsteps_(1),
351 : MCstride_(1),
352 : MCaccept_(0),
353 : MCacceptScale_(0),
354 : MCacceptFT_(0),
355 : MCtrial_(0),
356 : MCchunksize_(0),
357 : write_stride_(0),
358 : firstTime(true),
359 : do_reweight_(false),
360 : do_optsigmamean_(0),
361 : optsigmamean_stride_(0),
362 247 : average_weights_stride_(1)
363 : {
364 19 : bool noensemble = false;
365 38 : parseFlag("NOENSEMBLE", noensemble);
366 :
367 : // set up replica stuff
368 19 : master = (comm.Get_rank()==0);
369 19 : if(master) {
370 11 : nrep_ = multi_sim_comm.Get_size();
371 11 : replica_ = multi_sim_comm.Get_rank();
372 11 : if(noensemble) nrep_ = 1;
373 : } else {
374 8 : nrep_ = 0;
375 8 : replica_ = 0;
376 : }
377 19 : comm.Sum(&nrep_,1);
378 19 : comm.Sum(&replica_,1);
379 :
380 19 : unsigned nsel = 1;
381 38 : parse("SELECTOR", selector_);
382 38 : parse("NSELECT", nsel);
383 : // do checks
384 19 : if(selector_.length()>0 && nsel<=1) error("With SELECTOR active, NSELECT must be greater than 1");
385 19 : if(selector_.length()==0 && nsel>1) error("With NSELECT greater than 1, you must specify SELECTOR");
386 :
387 : // initialise firstTimeW
388 19 : firstTimeW.resize(nsel, true);
389 :
390 : // reweight implies a different number of arguments (the latest one must always be the bias)
391 38 : parseFlag("REWEIGHT", do_reweight_);
392 19 : if(do_reweight_&&nrep_<2) error("REWEIGHT can only be used in parallel with 2 or more replicas");
393 38 : if(!getRestart()) average_weights_.resize(nsel, vector<double> (nrep_, 1./static_cast<double>(nrep_)));
394 0 : else average_weights_.resize(nsel, vector<double> (nrep_, 0.));
395 19 : narg = getNumberOfArguments();
396 19 : if(do_reweight_) narg--;
397 :
398 19 : unsigned averaging=0;
399 38 : parse("AVERAGING", averaging);
400 19 : if(averaging>0) {
401 0 : average_weights_stride_ = averaging;
402 0 : optsigmamean_stride_ = averaging;
403 : }
404 :
405 38 : parseVector("PARAMETERS",parameters);
406 19 : if(parameters.size()!=static_cast<unsigned>(narg)&&!parameters.empty())
407 0 : error("Size of PARAMETERS array should be either 0 or the same as of the number of arguments in ARG1");
408 :
409 : vector<Value*> arg2;
410 38 : parseArgumentList("PARARG",arg2);
411 19 : if(!arg2.empty()) {
412 4 : if(parameters.size()>0) error("It is not possible to use PARARG and PARAMETERS together");
413 4 : if(arg2.size()!=narg) error("Size of PARARG array should be the same as number for arguments in ARG");
414 7076 : for(unsigned i=0; i<arg2.size(); i++) {
415 7068 : parameters.push_back(arg2[i]->get());
416 4712 : if(arg2[i]->hasDerivatives()==true) error("PARARG can only accept arguments without derivatives");
417 : }
418 : }
419 :
420 19 : if(parameters.size()!=narg)
421 0 : error("PARARG or PARAMETERS arrays should include the same number of elements as the arguments in ARG");
422 :
423 : string stringa_noise;
424 38 : parse("NOISETYPE",stringa_noise);
425 19 : if(stringa_noise=="GAUSS") noise_type_ = GAUSS;
426 18 : else if(stringa_noise=="MGAUSS") noise_type_ = MGAUSS;
427 10 : else if(stringa_noise=="OUTLIERS") noise_type_ = OUTLIERS;
428 5 : else if(stringa_noise=="MOUTLIERS") noise_type_ = MOUTLIERS;
429 1 : else if(stringa_noise=="GENERIC") noise_type_ = GENERIC;
430 0 : else error("Unknown noise type!");
431 :
432 19 : if(noise_type_== GENERIC) {
433 : string stringa_like;
434 2 : parse("LIKELIHOOD",stringa_like);
435 1 : if(stringa_like=="GAUSS") gen_likelihood_ = LIKE_GAUSS;
436 0 : else if(stringa_like=="LOGN") gen_likelihood_ = LIKE_LOGN;
437 0 : else error("Unknown likelihood type!");
438 :
439 2 : parse("DFTILDE",Dftilde_);
440 : }
441 :
442 38 : parse("WRITE_STRIDE",write_stride_);
443 : string status_file_name_;
444 38 : parse("STATUS_FILE",status_file_name_);
445 57 : if(status_file_name_=="") status_file_name_ = "MISTATUS"+getLabel();
446 0 : else status_file_name_ = status_file_name_+getLabel();
447 :
448 : string stringa_optsigma;
449 38 : parse("OPTSIGMAMEAN", stringa_optsigma);
450 19 : if(stringa_optsigma=="NONE") do_optsigmamean_=0;
451 0 : else if(stringa_optsigma=="SEM") do_optsigmamean_=1;
452 :
453 : // resize vector for sigma_mean history
454 19 : sigma_mean2_last_.resize(nsel);
455 57 : for(unsigned i=0; i<nsel; i++) sigma_mean2_last_[i].resize(narg);
456 :
457 : vector<double> read_sigma_mean_;
458 38 : parseVector("SIGMA_MEAN0",read_sigma_mean_);
459 38 : if(!do_optsigmamean_ && read_sigma_mean_.size()==0 && !getRestart())
460 0 : error("If you don't use OPTSIGMAMEAN and you are not RESTARTING then you MUST SET SIGMA_MEAN0");
461 :
462 19 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
463 13 : if(read_sigma_mean_.size()==narg) {
464 0 : sigma_mean2_.resize(narg);
465 0 : for(unsigned i=0; i<narg; i++) sigma_mean2_[i]=read_sigma_mean_[i]*read_sigma_mean_[i];
466 13 : } else if(read_sigma_mean_.size()==1) {
467 13 : sigma_mean2_.resize(narg,read_sigma_mean_[0]*read_sigma_mean_[0]);
468 0 : } else if(read_sigma_mean_.size()==0) {
469 0 : sigma_mean2_.resize(narg,0.000001);
470 : } else {
471 0 : error("SIGMA_MEAN0 can accept either one single value or as many values as the arguments (with NOISETYPE=MGAUSS|MOUTLIERS)");
472 : }
473 : // set the initial value for the history
474 2416 : for(unsigned i=0; i<nsel; i++) for(unsigned j=0; j<narg; j++) sigma_mean2_last_[i][j].push_back(sigma_mean2_[j]);
475 : } else {
476 6 : if(read_sigma_mean_.size()==1) {
477 6 : sigma_mean2_.resize(1, read_sigma_mean_[0]*read_sigma_mean_[0]);
478 0 : } else if(read_sigma_mean_.size()==0) {
479 0 : sigma_mean2_.resize(1, 0.000001);
480 : } else {
481 0 : error("If you want to use more than one SIGMA_MEAN0 you should use NOISETYPE=MGAUSS|MOUTLIERS");
482 : }
483 : // set the initial value for the history
484 37 : for(unsigned i=0; i<nsel; i++) for(unsigned j=0; j<narg; j++) sigma_mean2_last_[i][j].push_back(sigma_mean2_[0]);
485 : }
486 :
487 38 : parseFlag("SCALEDATA", doscale_);
488 19 : if(doscale_) {
489 : string stringa_noise;
490 24 : parse("SCALE_PRIOR",stringa_noise);
491 12 : if(stringa_noise=="GAUSSIAN") scale_prior_ = SC_GAUSS;
492 12 : else if(stringa_noise=="FLAT") scale_prior_ = SC_FLAT;
493 0 : else error("Unknown SCALE_PRIOR type!");
494 24 : parse("SCALE0",scale_);
495 24 : parse("DSCALE",Dscale_);
496 12 : if(Dscale_<0.) error("DSCALE must be set when using SCALEDATA");
497 12 : if(scale_prior_==SC_GAUSS) {
498 0 : scale_mu_=scale_;
499 : } else {
500 24 : parse("SCALE_MIN",scale_min_);
501 24 : parse("SCALE_MAX",scale_max_);
502 12 : if(scale_max_<scale_min_) error("SCALE_MAX and SCALE_MIN must be set when using SCALE_PRIOR=FLAT");
503 : }
504 : }
505 :
506 38 : parseFlag("ADDOFFSET", dooffset_);
507 19 : if(dooffset_) {
508 : string stringa_noise;
509 4 : parse("OFFSET_PRIOR",stringa_noise);
510 2 : if(stringa_noise=="GAUSSIAN") offset_prior_ = SC_GAUSS;
511 2 : else if(stringa_noise=="FLAT") offset_prior_ = SC_FLAT;
512 0 : else error("Unknown OFFSET_PRIOR type!");
513 4 : parse("OFFSET0",offset_);
514 4 : parse("DOFFSET",Doffset_);
515 2 : if(offset_prior_==SC_GAUSS) {
516 0 : offset_mu_=offset_;
517 0 : if(Doffset_<0.) error("DOFFSET must be set when using OFFSET_PRIOR=GAUSS");
518 : } else {
519 4 : parse("OFFSET_MIN",offset_min_);
520 4 : parse("OFFSET_MAX",offset_max_);
521 2 : if(Doffset_<0) Doffset_ = 0.05*(offset_max_ - offset_min_);
522 2 : if(offset_max_<offset_min_) error("OFFSET_MAX and OFFSET_MIN must be set when using OFFSET_PRIOR=FLAT");
523 : }
524 : }
525 :
526 : // regression with zero intercept
527 38 : parse("REGRES_ZERO", nregres_zero_);
528 19 : if(nregres_zero_>0) {
529 : // set flag
530 0 : doregres_zero_=true;
531 : // check if already sampling scale and offset
532 0 : if(doscale_) error("REGRES_ZERO and SCALEDATA are mutually exclusive");
533 0 : if(dooffset_) error("REGRES_ZERO and ADDOFFSET are mutually exclusive");
534 : }
535 :
536 : vector<double> readsigma;
537 38 : parseVector("SIGMA0",readsigma);
538 25 : if((noise_type_!=MGAUSS&&noise_type_!=MOUTLIERS&&noise_type_!=GENERIC)&&readsigma.size()>1)
539 0 : error("If you want to use more than one SIGMA you should use NOISETYPE=MGAUSS|MOUTLIERS|GENERIC");
540 19 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
541 13 : sigma_.resize(readsigma.size());
542 13 : sigma_=readsigma;
543 6 : } else sigma_.resize(1, readsigma[0]);
544 :
545 : vector<double> readsigma_min;
546 38 : parseVector("SIGMA_MIN",readsigma_min);
547 25 : if((noise_type_!=MGAUSS&&noise_type_!=MOUTLIERS&&noise_type_!=GENERIC)&&readsigma_min.size()>1)
548 0 : error("If you want to use more than one SIGMA you should use NOISETYPE=MGAUSS|MOUTLIERS|GENERIC");
549 19 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
550 13 : sigma_min_.resize(readsigma_min.size());
551 13 : sigma_min_=readsigma_min;
552 6 : } else sigma_min_.resize(1, readsigma_min[0]);
553 :
554 : vector<double> readsigma_max;
555 38 : parseVector("SIGMA_MAX",readsigma_max);
556 25 : if((noise_type_!=MGAUSS&&noise_type_!=MOUTLIERS&&noise_type_!=GENERIC)&&readsigma_max.size()>1)
557 0 : error("If you want to use more than one SIGMA you should use NOISETYPE=MGAUSS|MOUTLIERS|GENERIC");
558 19 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
559 13 : sigma_max_.resize(readsigma_max.size());
560 13 : sigma_max_=readsigma_max;
561 6 : } else sigma_max_.resize(1, readsigma_max[0]);
562 :
563 19 : if(sigma_max_.size()!=sigma_min_.size()) error("The number of values for SIGMA_MIN and SIGMA_MAX must be the same");
564 :
565 : vector<double> read_dsigma;
566 38 : parseVector("DSIGMA",read_dsigma);
567 25 : if((noise_type_!=MGAUSS&&noise_type_!=MOUTLIERS&&noise_type_!=GENERIC)&&readsigma_max.size()>1)
568 0 : error("If you want to use more than one SIGMA you should use NOISETYPE=MGAUSS|MOUTLIERS|GENERIC");
569 19 : if(read_dsigma.size()>0) {
570 19 : Dsigma_.resize(read_dsigma.size());
571 19 : Dsigma_=read_dsigma;
572 : } else {
573 0 : Dsigma_.resize(sigma_max_.size());
574 0 : for(unsigned i=0; i<sigma_max_.size(); i++) Dsigma_[i] = 0.05*(sigma_max_[i] - sigma_min_[i]);
575 : }
576 :
577 : // monte carlo stuff
578 38 : parse("MC_STEPS",MCsteps_);
579 38 : parse("MC_STRIDE",MCstride_);
580 38 : parse("MC_CHUNKSIZE", MCchunksize_);
581 : // adjust for multiple-time steps
582 19 : MCstride_ *= getStride();
583 : // get temperature
584 19 : double temp=0.0;
585 38 : parse("TEMP",temp);
586 38 : if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
587 0 : else kbt_=plumed.getAtoms().getKbT();
588 19 : if(kbt_==0.0) error("Unless the MD engine passes the temperature to plumed, you must specify it using TEMP");
589 :
590 19 : checkRead();
591 :
592 : // set sigma_bias
593 19 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
594 13 : if(sigma_.size()==1) {
595 13 : double tmp = sigma_[0];
596 13 : sigma_.resize(narg, tmp);
597 0 : } else if(sigma_.size()>1&&sigma_.size()!=narg) {
598 0 : error("SIGMA0 can accept either one single value or as many values as the number of arguments (with NOISETYPE=MGAUSS|MOUTLIERS|GENERIC)");
599 : }
600 13 : if(sigma_min_.size()==1) {
601 13 : double tmp = sigma_min_[0];
602 13 : sigma_min_.resize(narg, tmp);
603 0 : } else if(sigma_min_.size()>1&&sigma_min_.size()!=narg) {
604 0 : error("SIGMA_MIN can accept either one single value or as many values as the number of arguments (with NOISETYPE=MGAUSS|MOUTLIERS|GENERIC)");
605 : }
606 13 : if(sigma_max_.size()==1) {
607 13 : double tmp = sigma_max_[0];
608 13 : sigma_max_.resize(narg, tmp);
609 0 : } else if(sigma_max_.size()>1&&sigma_max_.size()!=narg) {
610 0 : error("SIGMA_MAX can accept either one single value or as many values as the number of arguments (with NOISETYPE=MGAUSS|MOUTLIERS|GENERIC)");
611 : }
612 13 : if(Dsigma_.size()==1) {
613 13 : double tmp = Dsigma_[0];
614 13 : Dsigma_.resize(narg, tmp);
615 0 : } else if(Dsigma_.size()>1&&Dsigma_.size()!=narg) {
616 0 : error("DSIGMA can accept either one single value or as many values as the number of arguments (with NOISETYPE=MGAUSS|MOUTLIERS|GENERIC)");
617 : }
618 : }
619 :
620 38 : IFile restart_sfile;
621 19 : restart_sfile.link(*this);
622 19 : if(getRestart()&&restart_sfile.FileExist(status_file_name_)) {
623 0 : firstTime = false;
624 0 : for(unsigned i=0; i<nsel; i++) firstTimeW[i] = false;
625 0 : restart_sfile.open(status_file_name_);
626 0 : log.printf(" Restarting from %s\n", status_file_name_.c_str());
627 : double dummy;
628 0 : if(restart_sfile.scanField("time",dummy)) {
629 : // nsel
630 0 : for(unsigned i=0; i<sigma_mean2_last_.size(); i++) {
631 : std::string msg_i;
632 0 : Tools::convert(i,msg_i);
633 : // narg
634 0 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
635 0 : for(unsigned j=0; j<narg; ++j) {
636 : std::string msg_j;
637 0 : Tools::convert(j,msg_j);
638 0 : std::string msg = msg_i+"_"+msg_j;
639 : double read_sm;
640 0 : restart_sfile.scanField("sigmaMean_"+msg,read_sm);
641 0 : sigma_mean2_last_[i][j][0] = read_sm*read_sm;
642 : }
643 : }
644 0 : if(noise_type_==GAUSS||noise_type_==OUTLIERS) {
645 : double read_sm;
646 : std::string msg_j;
647 0 : Tools::convert(0,msg_j);
648 0 : std::string msg = msg_i+"_"+msg_j;
649 0 : restart_sfile.scanField("sigmaMean_"+msg,read_sm);
650 0 : for(unsigned j=0; j<narg; j++) sigma_mean2_last_[i][j][0] = read_sm*read_sm;
651 : }
652 : }
653 :
654 0 : for(unsigned i=0; i<sigma_.size(); ++i) {
655 : std::string msg;
656 0 : Tools::convert(i,msg);
657 0 : restart_sfile.scanField("sigma_"+msg,sigma_[i]);
658 : }
659 0 : if(noise_type_==GENERIC) {
660 0 : for(unsigned i=0; i<ftilde_.size(); ++i) {
661 : std::string msg;
662 0 : Tools::convert(i,msg);
663 0 : restart_sfile.scanField("ftilde_"+msg,ftilde_[i]);
664 : }
665 : }
666 0 : restart_sfile.scanField("scale0_",scale_);
667 0 : restart_sfile.scanField("offset0_",offset_);
668 :
669 0 : for(unsigned i=0; i<nsel; i++) {
670 : std::string msg;
671 0 : Tools::convert(i,msg);
672 : double tmp_w;
673 0 : restart_sfile.scanField("weight_"+msg,tmp_w);
674 0 : if(master) {
675 0 : average_weights_[i][replica_] = tmp_w;
676 0 : if(nrep_>1) multi_sim_comm.Sum(&average_weights_[i][0], nrep_);
677 : }
678 0 : comm.Sum(&average_weights_[i][0], nrep_);
679 : }
680 :
681 : }
682 0 : restart_sfile.scanField();
683 0 : restart_sfile.close();
684 : }
685 :
686 19 : switch(noise_type_) {
687 1 : case GENERIC:
688 1 : log.printf(" with general metainference ");
689 1 : if(gen_likelihood_==LIKE_GAUSS) log.printf(" and a gaussian likelihood\n");
690 0 : else if(gen_likelihood_==LIKE_LOGN) log.printf(" and a log-normal likelihood\n");
691 1 : log.printf(" ensemble average parameter sampled with a step %lf of sigma_mean\n", Dftilde_);
692 : break;
693 1 : case GAUSS:
694 1 : log.printf(" with gaussian noise and a single noise parameter for all the data\n");
695 : break;
696 8 : case MGAUSS:
697 8 : log.printf(" with gaussian noise and a noise parameter for each data point\n");
698 : break;
699 5 : case OUTLIERS:
700 5 : log.printf(" with long tailed gaussian noise and a single noise parameter for all the data\n");
701 : break;
702 4 : case MOUTLIERS:
703 4 : log.printf(" with long tailed gaussian noise and a noise parameter for each data point\n");
704 : break;
705 : }
706 :
707 19 : if(doscale_) {
708 12 : log.printf(" sampling a common scaling factor with:\n");
709 12 : log.printf(" initial scale parameter %f\n",scale_);
710 12 : if(scale_prior_==SC_GAUSS) {
711 0 : log.printf(" gaussian prior with mean %f and width %f\n",scale_mu_,Dscale_);
712 : }
713 12 : if(scale_prior_==SC_FLAT) {
714 12 : log.printf(" flat prior between %f - %f\n",scale_min_,scale_max_);
715 12 : log.printf(" maximum MC move of scale parameter %f\n",Dscale_);
716 : }
717 : }
718 :
719 19 : if(dooffset_) {
720 2 : log.printf(" sampling a common offset with:\n");
721 2 : log.printf(" initial offset parameter %f\n",offset_);
722 2 : if(offset_prior_==SC_GAUSS) {
723 0 : log.printf(" gaussian prior with mean %f and width %f\n",offset_mu_,Doffset_);
724 : }
725 2 : if(offset_prior_==SC_FLAT) {
726 2 : log.printf(" flat prior between %f - %f\n",offset_min_,offset_max_);
727 2 : log.printf(" maximum MC move of offset parameter %f\n",Doffset_);
728 : }
729 : }
730 :
731 19 : if(doregres_zero_)
732 0 : log.printf(" doing regression with zero intercept with stride: %d\n", nregres_zero_);
733 :
734 19 : log.printf(" number of experimental data points %u\n",narg);
735 19 : log.printf(" number of replicas %u\n",nrep_);
736 19 : log.printf(" initial data uncertainties");
737 7226 : for(unsigned i=0; i<sigma_.size(); ++i) log.printf(" %f", sigma_[i]);
738 19 : log.printf("\n");
739 19 : log.printf(" minimum data uncertainties");
740 7226 : for(unsigned i=0; i<sigma_.size(); ++i) log.printf(" %f",sigma_min_[i]);
741 19 : log.printf("\n");
742 19 : log.printf(" maximum data uncertainties");
743 7226 : for(unsigned i=0; i<sigma_.size(); ++i) log.printf(" %f",sigma_max_[i]);
744 19 : log.printf("\n");
745 19 : log.printf(" maximum MC move of data uncertainties");
746 7226 : for(unsigned i=0; i<sigma_.size(); ++i) log.printf(" %f",Dsigma_[i]);
747 19 : log.printf("\n");
748 19 : log.printf(" temperature of the system %f\n",kbt_);
749 19 : log.printf(" MC steps %u\n",MCsteps_);
750 19 : log.printf(" MC stride %u\n",MCstride_);
751 19 : log.printf(" initial standard errors of the mean");
752 7226 : for(unsigned i=0; i<sigma_mean2_.size(); ++i) log.printf(" %f", sqrt(sigma_mean2_[i]));
753 19 : log.printf("\n");
754 :
755 19 : if(do_reweight_) {
756 32 : addComponent("biasDer");
757 32 : componentIsNotPeriodic("biasDer");
758 32 : addComponent("weight");
759 32 : componentIsNotPeriodic("weight");
760 : }
761 :
762 19 : if(doscale_ || doregres_zero_) {
763 24 : addComponent("scale");
764 24 : componentIsNotPeriodic("scale");
765 24 : valueScale=getPntrToComponent("scale");
766 : }
767 :
768 19 : if(dooffset_) {
769 4 : addComponent("offset");
770 4 : componentIsNotPeriodic("offset");
771 4 : valueOffset=getPntrToComponent("offset");
772 : }
773 :
774 19 : if(dooffset_||doscale_) {
775 28 : addComponent("acceptScale");
776 28 : componentIsNotPeriodic("acceptScale");
777 28 : valueAcceptScale=getPntrToComponent("acceptScale");
778 : }
779 :
780 19 : if(noise_type_==GENERIC) {
781 2 : addComponent("acceptFT");
782 2 : componentIsNotPeriodic("acceptFT");
783 2 : valueAcceptFT=getPntrToComponent("acceptFT");
784 : }
785 :
786 38 : addComponent("acceptSigma");
787 38 : componentIsNotPeriodic("acceptSigma");
788 38 : valueAccept=getPntrToComponent("acceptSigma");
789 :
790 19 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
791 7196 : for(unsigned i=0; i<sigma_mean2_.size(); ++i) {
792 2390 : std::string num; Tools::convert(i,num);
793 7170 : addComponent("sigmaMean_"+num); componentIsNotPeriodic("sigmaMean_"+num);
794 7170 : valueSigmaMean.push_back(getPntrToComponent("sigmaMean_"+num));
795 7170 : getPntrToComponent("sigmaMean_"+num)->set(sqrt(sigma_mean2_[i]));
796 7170 : addComponent("sigma_"+num); componentIsNotPeriodic("sigma_"+num);
797 7170 : valueSigma.push_back(getPntrToComponent("sigma_"+num));
798 4780 : getPntrToComponent("sigma_"+num)->set(sigma_[i]);
799 2390 : if(noise_type_==GENERIC) {
800 6 : addComponent("ftilde_"+num); componentIsNotPeriodic("ftilde_"+num);
801 6 : valueFtilde.push_back(getPntrToComponent("ftilde_"+num));
802 : }
803 : }
804 : } else {
805 18 : addComponent("sigmaMean"); componentIsNotPeriodic("sigmaMean");
806 18 : valueSigmaMean.push_back(getPntrToComponent("sigmaMean"));
807 18 : getPntrToComponent("sigmaMean")->set(sqrt(sigma_mean2_[0]));
808 18 : addComponent("sigma"); componentIsNotPeriodic("sigma");
809 18 : valueSigma.push_back(getPntrToComponent("sigma"));
810 18 : getPntrToComponent("sigma")->set(sigma_[0]);
811 : }
812 :
813 : // initialize random seed
814 : unsigned iseed;
815 19 : if(master) iseed = time(NULL)+replica_;
816 8 : else iseed = 0;
817 19 : comm.Sum(&iseed, 1);
818 19 : random[0].setSeed(-iseed);
819 : // Random chunk
820 19 : if(master) iseed = time(NULL)+replica_;
821 8 : else iseed = 0;
822 19 : comm.Sum(&iseed, 1);
823 19 : random[2].setSeed(-iseed);
824 19 : if(doscale_||dooffset_) {
825 : // in this case we want the same seed everywhere
826 14 : iseed = time(NULL);
827 14 : if(master&&nrep_>1) multi_sim_comm.Bcast(iseed,0);
828 14 : comm.Bcast(iseed,0);
829 14 : random[1].setSeed(-iseed);
830 : }
831 :
832 : // outfile stuff
833 19 : if(write_stride_>0) {
834 19 : sfile_.link(*this);
835 19 : sfile_.open(status_file_name_);
836 : }
837 :
838 57 : log<<" Bibliography "<<plumed.cite("Bonomi, Camilloni, Cavalli, Vendruscolo, Sci. Adv. 2, e150117 (2016)");
839 51 : if(do_reweight_) log<<plumed.cite("Bonomi, Camilloni, Vendruscolo, Sci. Rep. 6, 31232 (2016)");
840 19 : if(do_optsigmamean_>0) log<<plumed.cite("Loehr, Jussupow, Camilloni, J. Chem. Phys. 146, 165102 (2017)");
841 57 : log<<plumed.cite("Bonomi, Camilloni, Bioinformatics, 33, 3999 (2017)");
842 19 : log<<"\n";
843 19 : }
844 :
845 114 : Metainference::~Metainference()
846 : {
847 19 : if(sfile_.isOpen()) sfile_.close();
848 38 : }
849 :
850 156 : double Metainference::getEnergySP(const vector<double> &mean, const vector<double> &sigma,
851 : const double scale, const double offset)
852 : {
853 156 : const double scale2 = scale*scale;
854 156 : const double sm2 = sigma_mean2_[0];
855 156 : const double ss2 = sigma[0]*sigma[0] + scale2*sm2;
856 156 : const double sss = sigma[0]*sigma[0] + sm2;
857 :
858 : double ene = 0.0;
859 468 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(ene)
860 : {
861 312 : #pragma omp for reduction( + : ene)
862 312 : for(unsigned i=0; i<narg; ++i) {
863 1980 : const double dev = scale*mean[i]-parameters[i]+offset;
864 660 : const double a2 = 0.5*dev*dev + ss2;
865 660 : ene += std::log(2.0*a2/(1.0-exp(-a2/sm2)));
866 : }
867 : }
868 : // add one single Jeffrey's prior and one normalisation per data point
869 156 : ene += 0.5*std::log(sss) + static_cast<double>(narg)*0.5*std::log(0.5*M_PI*M_PI/ss2);
870 156 : if(doscale_ || doregres_zero_) ene += 0.5*std::log(sss);
871 156 : if(dooffset_) ene += 0.5*std::log(sss);
872 156 : return kbt_ * ene;
873 : }
874 :
875 144 : double Metainference::getEnergySPE(const vector<double> &mean, const vector<double> &sigma,
876 : const double scale, const double offset)
877 : {
878 144 : const double scale2 = scale*scale;
879 : double ene = 0.0;
880 432 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(ene)
881 : {
882 288 : #pragma omp for reduction( + : ene)
883 288 : for(unsigned i=0; i<narg; ++i) {
884 1152 : const double sm2 = sigma_mean2_[i];
885 1152 : const double ss2 = sigma[i]*sigma[i] + scale2*sm2;
886 576 : const double sss = sigma[i]*sigma[i] + sm2;
887 1728 : const double dev = scale*mean[i]-parameters[i]+offset;
888 576 : const double a2 = 0.5*dev*dev + ss2;
889 576 : ene += 0.5*std::log(sss) + 0.5*std::log(0.5*M_PI*M_PI/ss2) + std::log(2.0*a2/(1.0-exp(-a2/sm2)));
890 1152 : if(doscale_ || doregres_zero_) ene += 0.5*std::log(sss);
891 576 : if(dooffset_) ene += 0.5*std::log(sss);
892 : }
893 : }
894 144 : return kbt_ * ene;
895 : }
896 :
897 48 : double Metainference::getEnergyMIGEN(const vector<double> &mean, const vector<double> &ftilde, const vector<double> &sigma,
898 : const double scale, const double offset)
899 : {
900 : double ene = 0.0;
901 144 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(ene)
902 : {
903 96 : #pragma omp for reduction( + : ene)
904 96 : for(unsigned i=0; i<narg; ++i) {
905 192 : const double inv_sb2 = 1./(sigma[i]*sigma[i]);
906 96 : const double inv_sm2 = 1./sigma_mean2_[i];
907 : double devb = 0;
908 384 : if(gen_likelihood_==LIKE_GAUSS) devb = scale*ftilde[i]-parameters[i]+offset;
909 0 : else if(gen_likelihood_==LIKE_LOGN) devb = std::log(scale*ftilde[i]/parameters[i]);
910 288 : double devm = mean[i] - ftilde[i];
911 : // deviation + normalisation + jeffrey
912 : double normb = 0.;
913 192 : if(gen_likelihood_==LIKE_GAUSS) normb = -0.5*std::log(0.5/M_PI*inv_sb2);
914 0 : else if(gen_likelihood_==LIKE_LOGN) normb = -0.5*std::log(0.5/M_PI*inv_sb2/(parameters[i]*parameters[i]));
915 96 : const double normm = -0.5*std::log(0.5/M_PI*inv_sm2);
916 96 : const double jeffreys = -0.5*std::log(2.*inv_sb2);
917 96 : ene += 0.5*devb*devb*inv_sb2 + 0.5*devm*devm*inv_sm2 + normb + normm + jeffreys;
918 96 : if(doscale_ || doregres_zero_) ene += jeffreys;
919 96 : if(dooffset_) ene += jeffreys;
920 : }
921 : }
922 48 : return kbt_ * ene;
923 : }
924 :
925 36 : double Metainference::getEnergyGJ(const vector<double> &mean, const vector<double> &sigma,
926 : const double scale, const double offset)
927 : {
928 36 : const double scale2 = scale*scale;
929 72 : const double inv_s2 = 1./(sigma[0]*sigma[0] + scale2*sigma_mean2_[0]);
930 36 : const double inv_sss = 1./(sigma[0]*sigma[0] + sigma_mean2_[0]);
931 :
932 : double ene = 0.0;
933 108 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(ene)
934 : {
935 72 : #pragma omp for reduction( + : ene)
936 72 : for(unsigned i=0; i<narg; ++i) {
937 216 : double dev = scale*mean[i]-parameters[i]+offset;
938 72 : ene += 0.5*dev*dev*inv_s2;
939 : }
940 : }
941 36 : const double normalisation = -0.5*std::log(0.5/M_PI*inv_s2);
942 36 : const double jeffreys = -0.5*std::log(2.*inv_sss);
943 : // add Jeffrey's prior in case one sigma for all data points + one normalisation per datapoint
944 36 : ene += jeffreys + static_cast<double>(narg)*normalisation;
945 36 : if(doscale_ || doregres_zero_) ene += jeffreys;
946 36 : if(dooffset_) ene += jeffreys;
947 :
948 36 : return kbt_ * ene;
949 : }
950 :
951 152 : double Metainference::getEnergyGJE(const vector<double> &mean, const vector<double> &sigma,
952 : const double scale, const double offset)
953 : {
954 152 : const double scale2 = scale*scale;
955 :
956 : double ene = 0.0;
957 456 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(ene)
958 : {
959 304 : #pragma omp for reduction( + : ene)
960 304 : for(unsigned i=0; i<narg; ++i) {
961 15864 : const double inv_s2 = 1./(sigma[i]*sigma[i] + scale2*sigma_mean2_[i]);
962 5288 : const double inv_sss = 1./(sigma[i]*sigma[i] + sigma_mean2_[i]);
963 15864 : double dev = scale*mean[i]-parameters[i]+offset;
964 : // deviation + normalisation + jeffrey
965 5288 : const double normalisation = -0.5*std::log(0.5/M_PI*inv_s2);
966 5288 : const double jeffreys = -0.5*std::log(2.*inv_sss);
967 5288 : ene += 0.5*dev*dev*inv_s2 + normalisation + jeffreys;
968 5864 : if(doscale_ || doregres_zero_) ene += jeffreys;
969 5288 : if(dooffset_) ene += jeffreys;
970 : }
971 : }
972 152 : return kbt_ * ene;
973 : }
974 :
975 178 : void Metainference::doMonteCarlo(const vector<double> &mean_)
976 : {
977 : // calculate old energy with the updated coordinates
978 : double old_energy=0.;
979 :
980 178 : switch(noise_type_) {
981 12 : case GAUSS:
982 12 : old_energy = getEnergyGJ(mean_,sigma_,scale_,offset_);
983 : break;
984 52 : case MGAUSS:
985 52 : old_energy = getEnergyGJE(mean_,sigma_,scale_,offset_);
986 : break;
987 54 : case OUTLIERS:
988 54 : old_energy = getEnergySP(mean_,sigma_,scale_,offset_);
989 : break;
990 48 : case MOUTLIERS:
991 48 : old_energy = getEnergySPE(mean_,sigma_,scale_,offset_);
992 : break;
993 12 : case GENERIC:
994 12 : old_energy = getEnergyMIGEN(mean_,ftilde_,sigma_,scale_,offset_);
995 : break;
996 : }
997 :
998 : // Create vector of random sigma indices
999 : vector<unsigned> indices;
1000 178 : if (MCchunksize_ > 0) {
1001 0 : for (unsigned j=0; j<sigma_.size(); j++) {
1002 0 : indices.push_back(j);
1003 : }
1004 0 : random[2].Shuffle(indices);
1005 : }
1006 : bool breaknow = false;
1007 :
1008 : // cycle on MC steps
1009 534 : for(unsigned i=0; i<MCsteps_; ++i) {
1010 :
1011 178 : MCtrial_++;
1012 :
1013 : // propose move for ftilde
1014 178 : vector<double> new_ftilde(sigma_.size());
1015 178 : new_ftilde = ftilde_;
1016 :
1017 178 : if(noise_type_==GENERIC) {
1018 : // change all sigmas
1019 96 : for(unsigned j=0; j<sigma_.size(); j++) {
1020 24 : const double r3 = random[0].Gaussian();
1021 48 : const double ds3 = Dftilde_*sqrt(sigma_mean2_[j])*r3;
1022 48 : new_ftilde[j] = ftilde_[j] + ds3;
1023 : }
1024 : // calculate new energy
1025 12 : double new_energy = getEnergyMIGEN(mean_,new_ftilde,sigma_,scale_,offset_);
1026 :
1027 : // accept or reject
1028 12 : const double delta = ( new_energy - old_energy ) / kbt_;
1029 : // if delta is negative always accept move
1030 12 : if( delta <= 0.0 ) {
1031 : old_energy = new_energy;
1032 0 : ftilde_ = new_ftilde;
1033 0 : MCacceptFT_++;
1034 : // otherwise extract random number
1035 : } else {
1036 12 : const double s = random[0].RandU01();
1037 12 : if( s < exp(-delta) ) {
1038 : old_energy = new_energy;
1039 0 : ftilde_ = new_ftilde;
1040 0 : MCacceptFT_++;
1041 : }
1042 : }
1043 : }
1044 :
1045 : // propose move for scale and/or offset
1046 178 : double new_scale = scale_;
1047 178 : double new_offset = offset_;
1048 178 : if(doscale_||dooffset_) {
1049 168 : if(doscale_) {
1050 144 : if(scale_prior_==SC_FLAT) {
1051 144 : const double r1 = random[1].Gaussian();
1052 144 : const double ds1 = Dscale_*r1;
1053 144 : new_scale += ds1;
1054 : // check boundaries
1055 144 : if(new_scale > scale_max_) {new_scale = 2.0 * scale_max_ - new_scale;}
1056 144 : if(new_scale < scale_min_) {new_scale = 2.0 * scale_min_ - new_scale;}
1057 : } else {
1058 0 : const double r1 = random[1].Gaussian();
1059 0 : const double ds1 = 0.5*(scale_mu_-new_scale)+Dscale_*exp(1)/M_PI*r1;
1060 0 : new_scale += ds1;
1061 : }
1062 : }
1063 :
1064 168 : if(dooffset_) {
1065 24 : if(offset_prior_==SC_FLAT) {
1066 24 : const double r1 = random[1].Gaussian();
1067 24 : const double ds1 = Doffset_*r1;
1068 24 : new_offset += ds1;
1069 : // check boundaries
1070 24 : if(new_offset > offset_max_) {new_offset = 2.0 * offset_max_ - new_offset;}
1071 24 : if(new_offset < offset_min_) {new_offset = 2.0 * offset_min_ - new_offset;}
1072 : } else {
1073 0 : const double r1 = random[1].Gaussian();
1074 0 : const double ds1 = 0.5*(offset_mu_-new_offset)+Doffset_*exp(1)/M_PI*r1;
1075 0 : new_offset += ds1;
1076 : }
1077 : }
1078 :
1079 : // calculate new energy
1080 : double new_energy = 0.;
1081 :
1082 168 : switch(noise_type_) {
1083 12 : case GAUSS:
1084 12 : new_energy = getEnergyGJ(mean_,sigma_,new_scale,new_offset);
1085 : break;
1086 48 : case MGAUSS:
1087 48 : new_energy = getEnergyGJE(mean_,sigma_,new_scale,new_offset);
1088 : break;
1089 48 : case OUTLIERS:
1090 48 : new_energy = getEnergySP(mean_,sigma_,new_scale,new_offset);
1091 : break;
1092 48 : case MOUTLIERS:
1093 48 : new_energy = getEnergySPE(mean_,sigma_,new_scale,new_offset);
1094 : break;
1095 12 : case GENERIC:
1096 12 : new_energy = getEnergyMIGEN(mean_,ftilde_,sigma_,new_scale,new_offset);
1097 : break;
1098 : }
1099 : // for the scale we need to consider the total energy
1100 168 : vector<double> totenergies(2);
1101 168 : if(master) {
1102 96 : totenergies[0] = old_energy;
1103 96 : totenergies[1] = new_energy;
1104 96 : if(nrep_>1) multi_sim_comm.Sum(totenergies);
1105 : } else {
1106 72 : totenergies[0] = 0;
1107 72 : totenergies[1] = 0;
1108 : }
1109 168 : comm.Sum(totenergies);
1110 :
1111 : // accept or reject
1112 168 : const double delta = ( totenergies[1] - totenergies[0] ) / kbt_;
1113 : // if delta is negative always accept move
1114 168 : if( delta <= 0.0 ) {
1115 : old_energy = new_energy;
1116 156 : scale_ = new_scale;
1117 156 : offset_ = new_offset;
1118 156 : MCacceptScale_++;
1119 : // otherwise extract random number
1120 : } else {
1121 12 : double s = random[1].RandU01();
1122 12 : if( s < exp(-delta) ) {
1123 : old_energy = new_energy;
1124 0 : scale_ = new_scale;
1125 0 : offset_ = new_offset;
1126 0 : MCacceptScale_++;
1127 : }
1128 : }
1129 : }
1130 :
1131 : // propose move for sigma
1132 178 : vector<double> new_sigma(sigma_.size());
1133 178 : new_sigma = sigma_;
1134 :
1135 : // change MCchunksize_ sigmas
1136 178 : if (MCchunksize_ > 0) {
1137 0 : if ((MCchunksize_ * i) >= sigma_.size()) {
1138 : // This means we are not moving any sigma, so we should break immediately
1139 : breaknow = true;
1140 : }
1141 :
1142 : // change random sigmas
1143 0 : for(unsigned j=0; j<MCchunksize_; j++) {
1144 0 : const unsigned shuffle_index = j + MCchunksize_ * i;
1145 0 : if (shuffle_index >= sigma_.size()) {
1146 : // Going any further will segfault but we should still evaluate the sigmas we changed
1147 : break;
1148 : }
1149 0 : const unsigned index = indices[shuffle_index];
1150 0 : const double r2 = random[0].Gaussian();
1151 0 : const double ds2 = Dsigma_[index]*r2;
1152 0 : new_sigma[index] = sigma_[index] + ds2;
1153 : // check boundaries
1154 0 : if(new_sigma[index] > sigma_max_[index]) {new_sigma[index] = 2.0 * sigma_max_[index] - new_sigma[index];}
1155 0 : if(new_sigma[index] < sigma_min_[index]) {new_sigma[index] = 2.0 * sigma_min_[index] - new_sigma[index];}
1156 : }
1157 : } else {
1158 : // change all sigmas
1159 8846 : for(unsigned j=0; j<sigma_.size(); j++) {
1160 2830 : const double r2 = random[0].Gaussian();
1161 2830 : const double ds2 = Dsigma_[j]*r2;
1162 5660 : new_sigma[j] = sigma_[j] + ds2;
1163 : // check boundaries
1164 5660 : if(new_sigma[j] > sigma_max_[j]) {new_sigma[j] = 2.0 * sigma_max_[j] - new_sigma[j];}
1165 5660 : if(new_sigma[j] < sigma_min_[j]) {new_sigma[j] = 2.0 * sigma_min_[j] - new_sigma[j];}
1166 : }
1167 : }
1168 :
1169 178 : if (breaknow) {
1170 : // We didnt move any sigmas, so no sense in evaluating anything
1171 : break;
1172 : }
1173 :
1174 : // calculate new energy
1175 : double new_energy = 0.;
1176 178 : switch(noise_type_) {
1177 12 : case GAUSS:
1178 12 : new_energy = getEnergyGJ(mean_,new_sigma,scale_,offset_);
1179 : break;
1180 52 : case MGAUSS:
1181 52 : new_energy = getEnergyGJE(mean_,new_sigma,scale_,offset_);
1182 : break;
1183 54 : case OUTLIERS:
1184 54 : new_energy = getEnergySP(mean_,new_sigma,scale_,offset_);
1185 : break;
1186 48 : case MOUTLIERS:
1187 48 : new_energy = getEnergySPE(mean_,new_sigma,scale_,offset_);
1188 : break;
1189 12 : case GENERIC:
1190 12 : new_energy = getEnergyMIGEN(mean_,ftilde_,new_sigma,scale_,offset_);
1191 : break;
1192 : }
1193 :
1194 : // accept or reject
1195 178 : const double delta = ( new_energy - old_energy ) / kbt_;
1196 : // if delta is negative always accept move
1197 178 : if( delta <= 0.0 ) {
1198 : old_energy = new_energy;
1199 166 : sigma_ = new_sigma;
1200 166 : MCaccept_++;
1201 : // otherwise extract random number
1202 : } else {
1203 12 : const double s = random[0].RandU01();
1204 12 : if( s < exp(-delta) ) {
1205 : old_energy = new_energy;
1206 0 : sigma_ = new_sigma;
1207 0 : MCaccept_++;
1208 : }
1209 : }
1210 :
1211 : }
1212 : /* save the result of the sampling */
1213 178 : double accept = static_cast<double>(MCaccept_) / static_cast<double>(MCtrial_);
1214 178 : valueAccept->set(accept);
1215 178 : if(doscale_ || doregres_zero_) valueScale->set(scale_);
1216 178 : if(dooffset_) valueOffset->set(offset_);
1217 178 : if(doscale_||dooffset_) {
1218 168 : accept = static_cast<double>(MCacceptScale_) / static_cast<double>(MCtrial_);
1219 168 : valueAcceptScale->set(accept);
1220 : }
1221 14506 : for(unsigned i=0; i<sigma_.size(); i++) valueSigma[i]->set(sigma_[i]);
1222 178 : if(noise_type_==GENERIC) {
1223 12 : accept = static_cast<double>(MCacceptFT_) / static_cast<double>(MCtrial_);
1224 12 : valueAcceptFT->set(accept);
1225 144 : for(unsigned i=0; i<sigma_.size(); i++) valueFtilde[i]->set(ftilde_[i]);
1226 : }
1227 178 : }
1228 :
1229 : /*
1230 : In the following energy-force functions we don't add the normalisation and the jeffreys priors
1231 : because they are not needed for the forces, the correct MetaInference energy is the one calculated
1232 : in the Monte-Carlo
1233 : */
1234 :
1235 54 : double Metainference::getEnergyForceSP(const vector<double> &mean, const vector<double> &dmean_x,
1236 : const vector<double> &dmean_b)
1237 : {
1238 54 : const double scale2 = scale_*scale_;
1239 54 : const double sm2 = sigma_mean2_[0];
1240 54 : const double ss2 = sigma_[0]*sigma_[0] + scale2*sm2;
1241 54 : vector<double> f(narg+1,0);
1242 :
1243 54 : if(master) {
1244 : double omp_ene=0.;
1245 90 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(omp_ene)
1246 : {
1247 60 : #pragma omp for reduction( + : omp_ene)
1248 60 : for(unsigned i=0; i<narg; ++i) {
1249 414 : const double dev = scale_*mean[i]-parameters[i]+offset_;
1250 138 : const double a2 = 0.5*dev*dev + ss2;
1251 138 : const double t = exp(-a2/sm2);
1252 138 : const double dt = 1./t;
1253 138 : const double it = 1./(1.-t);
1254 138 : const double dit = 1./(1.-dt);
1255 138 : omp_ene += std::log(2.*a2*it);
1256 276 : f[i] = -scale_*dev*(dit/sm2 + 1./a2);
1257 : }
1258 : }
1259 60 : f[narg] = omp_ene;
1260 : // collect contribution to forces and energy from other replicas
1261 54 : if(nrep_>1) multi_sim_comm.Sum(&f[0],narg+1);
1262 : }
1263 : // intra-replica summation
1264 108 : comm.Sum(&f[0],narg+1);
1265 :
1266 108 : const double ene = f[narg];
1267 : double w_tmp = 0.;
1268 522 : for(unsigned i=0; i<narg; ++i) {
1269 702 : setOutputForce(i, kbt_*f[i]*dmean_x[i]);
1270 702 : w_tmp += kbt_*f[i]*dmean_b[i];
1271 : }
1272 :
1273 54 : if(do_reweight_) {
1274 48 : setOutputForce(narg, w_tmp);
1275 96 : getPntrToComponent("biasDer")->set(-w_tmp);
1276 : }
1277 :
1278 108 : return kbt_*ene;
1279 : }
1280 :
1281 48 : double Metainference::getEnergyForceSPE(const vector<double> &mean, const vector<double> &dmean_x,
1282 : const vector<double> &dmean_b)
1283 : {
1284 48 : const double scale2 = scale_*scale_;
1285 48 : vector<double> f(narg+1,0);
1286 :
1287 48 : if(master) {
1288 : double omp_ene = 0;
1289 72 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(omp_ene)
1290 : {
1291 48 : #pragma omp for reduction( + : omp_ene)
1292 48 : for(unsigned i=0; i<narg; ++i) {
1293 192 : const double sm2 = sigma_mean2_[i];
1294 96 : const double ss2 = sigma_[i]*sigma_[i] + scale2*sm2;
1295 288 : const double dev = scale_*mean[i]-parameters[i]+offset_;
1296 96 : const double a2 = 0.5*dev*dev + ss2;
1297 96 : const double t = exp(-a2/sm2);
1298 96 : const double dt = 1./t;
1299 96 : const double it = 1./(1.-t);
1300 96 : const double dit = 1./(1.-dt);
1301 96 : omp_ene += std::log(2.*a2*it);
1302 192 : f[i] = -scale_*dev*(dit/sm2 + 1./a2);
1303 : }
1304 : }
1305 48 : f[narg] = omp_ene;
1306 : // collect contribution to forces and energy from other replicas
1307 48 : if(nrep_>1) multi_sim_comm.Sum(&f[0],narg+1);
1308 : }
1309 96 : comm.Sum(&f[0],narg+1);
1310 :
1311 96 : const double ene = f[narg];
1312 : double w_tmp = 0.;
1313 432 : for(unsigned i=0; i<narg; ++i) {
1314 576 : setOutputForce(i, kbt_ * dmean_x[i] * f[i]);
1315 576 : w_tmp += kbt_ * dmean_b[i] *f[i];
1316 : }
1317 :
1318 48 : if(do_reweight_) {
1319 48 : setOutputForce(narg, w_tmp);
1320 96 : getPntrToComponent("biasDer")->set(-w_tmp);
1321 : }
1322 :
1323 96 : return kbt_*ene;
1324 : }
1325 :
1326 12 : double Metainference::getEnergyForceGJ(const vector<double> &mean, const vector<double> &dmean_x,
1327 : const vector<double> &dmean_b)
1328 : {
1329 12 : const double scale2 = scale_*scale_;
1330 12 : double inv_s2=0.;
1331 :
1332 12 : if(master) {
1333 24 : inv_s2 = 1./(sigma_[0]*sigma_[0] + scale2*sigma_mean2_[0]);
1334 12 : if(nrep_>1) multi_sim_comm.Sum(inv_s2);
1335 : }
1336 12 : comm.Sum(inv_s2);
1337 :
1338 : double ene = 0.;
1339 : double w_tmp = 0.;
1340 36 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(ene,w_tmp)
1341 : {
1342 24 : #pragma omp for reduction( + : ene,w_tmp)
1343 24 : for(unsigned i=0; i<narg; ++i) {
1344 72 : const double dev = scale_*mean[i]-parameters[i]+offset_;
1345 24 : const double mult = dev*scale_*inv_s2;
1346 24 : ene += 0.5*dev*dev*inv_s2;
1347 48 : setOutputForce(i, -kbt_*dmean_x[i]*mult);
1348 48 : w_tmp += kbt_*dmean_b[i]*mult;
1349 : }
1350 : }
1351 :
1352 12 : if(do_reweight_) {
1353 0 : setOutputForce(narg, -w_tmp);
1354 0 : getPntrToComponent("biasDer")->set(w_tmp);
1355 : }
1356 :
1357 12 : return kbt_*ene;
1358 : }
1359 :
1360 52 : double Metainference::getEnergyForceGJE(const vector<double> &mean, const vector<double> &dmean_x,
1361 : const vector<double> &dmean_b)
1362 : {
1363 52 : const double scale2 = scale_*scale_;
1364 104 : vector<double> inv_s2(sigma_.size(),0.);
1365 :
1366 52 : if(master) {
1367 6422 : for(unsigned i=0; i<sigma_.size(); ++i) inv_s2[i] = 1./(sigma_[i]*sigma_[i] + scale2*sigma_mean2_[i]);
1368 52 : if(nrep_>1) multi_sim_comm.Sum(&inv_s2[0],sigma_.size());
1369 : }
1370 156 : comm.Sum(&inv_s2[0],sigma_.size());
1371 :
1372 : double ene = 0.;
1373 : double w_tmp = 0.;
1374 156 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(ene,w_tmp)
1375 : {
1376 104 : #pragma omp for reduction( + : ene,w_tmp)
1377 104 : for(unsigned i=0; i<narg; ++i) {
1378 7644 : const double dev = scale_*mean[i]-parameters[i]+offset_;
1379 5096 : const double mult = dev*scale_*inv_s2[i];
1380 2548 : ene += 0.5*dev*dev*inv_s2[i];
1381 5096 : setOutputForce(i, -kbt_*dmean_x[i]*mult);
1382 5096 : w_tmp += kbt_*dmean_b[i]*mult;
1383 : }
1384 : }
1385 :
1386 52 : if(do_reweight_) {
1387 52 : setOutputForce(narg, -w_tmp);
1388 104 : getPntrToComponent("biasDer")->set(w_tmp);
1389 : }
1390 :
1391 104 : return kbt_*ene;
1392 : }
1393 :
1394 12 : double Metainference::getEnergyForceMIGEN(const vector<double> &mean, const vector<double> &dmean_x, const vector<double> &dmean_b)
1395 : {
1396 24 : vector<double> inv_s2(sigma_.size(),0.);
1397 24 : vector<double> dev(sigma_.size(),0.);
1398 24 : vector<double> dev2(sigma_.size(),0.);
1399 :
1400 96 : for(unsigned i=0; i<sigma_.size(); ++i) {
1401 48 : inv_s2[i] = 1./sigma_mean2_[i];
1402 24 : if(master) {
1403 72 : dev[i] = (mean[i]-ftilde_[i]);
1404 48 : dev2[i] = dev[i]*dev[i];
1405 : }
1406 : }
1407 12 : if(master&&nrep_>1) {
1408 0 : multi_sim_comm.Sum(&dev[0],dev.size());
1409 0 : multi_sim_comm.Sum(&dev2[0],dev2.size());
1410 : }
1411 24 : comm.Sum(&dev[0],dev.size());
1412 24 : comm.Sum(&dev2[0],dev2.size());
1413 :
1414 : double dene_b = 0.;
1415 : double ene = 0.;
1416 36 : #pragma omp parallel num_threads(OpenMP::getNumThreads()) shared(ene,dene_b)
1417 : {
1418 24 : #pragma omp for reduction( + : ene,dene_b)
1419 24 : for(unsigned i=0; i<narg; ++i) {
1420 96 : const double dene_x = kbt_*inv_s2[i]*dmean_x[i]*dev[i];
1421 48 : dene_b += kbt_*inv_s2[i]*dmean_b[i]*dev[i];
1422 48 : ene += 0.5*dev2[i]*inv_s2[i];
1423 24 : setOutputForce(i, -dene_x);
1424 : }
1425 : }
1426 :
1427 12 : if(do_reweight_) {
1428 0 : setOutputForce(narg, -dene_b);
1429 0 : getPntrToComponent("biasDer")->set(dene_b);
1430 : }
1431 :
1432 24 : return kbt_*ene;
1433 : }
1434 :
1435 178 : void Metainference::get_weights(const unsigned iselect, double &fact, double &var_fact)
1436 : {
1437 178 : const double dnrep = static_cast<double>(nrep_);
1438 178 : const double ave_fact = 1.0/dnrep;
1439 :
1440 : double norm = 0.0;
1441 :
1442 : // calculate the weights either from BIAS
1443 178 : if(do_reweight_) {
1444 148 : vector<double> bias(nrep_,0);
1445 148 : if(master) {
1446 148 : bias[replica_] = getArgument(narg);
1447 148 : if(nrep_>1) multi_sim_comm.Sum(&bias[0], nrep_);
1448 : }
1449 296 : comm.Sum(&bias[0], nrep_);
1450 :
1451 148 : const double maxbias = *(std::max_element(bias.begin(), bias.end()));
1452 740 : for(unsigned i=0; i<nrep_; ++i) {
1453 592 : bias[i] = exp((bias[i]-maxbias)/kbt_);
1454 296 : norm += bias[i];
1455 : }
1456 :
1457 : // accumulate weights
1458 148 : const double decay = 1./static_cast<double> (average_weights_stride_);
1459 296 : if(!firstTimeW[iselect]) {
1460 660 : for(unsigned i=0; i<nrep_; ++i) {
1461 792 : const double delta=bias[i]/norm-average_weights_[iselect][i];
1462 264 : average_weights_[iselect][i]+=decay*delta;
1463 : }
1464 : } else {
1465 : firstTimeW[iselect] = false;
1466 80 : for(unsigned i=0; i<nrep_; ++i) {
1467 96 : average_weights_[iselect][i] = bias[i]/norm;
1468 : }
1469 : }
1470 :
1471 : // set average back into bias and set norm to one
1472 740 : for(unsigned i=0; i<nrep_; ++i) bias[i] = average_weights_[iselect][i];
1473 : // set local weight, norm and weight variance
1474 296 : fact = bias[replica_];
1475 : norm = 1.0;
1476 740 : for(unsigned i=0; i<nrep_; ++i) var_fact += (bias[i]/norm-ave_fact)*(bias[i]/norm-ave_fact);
1477 296 : getPntrToComponent("weight")->set(fact);
1478 : } else {
1479 : // or arithmetic ones
1480 : norm = dnrep;
1481 30 : fact = 1.0/norm;
1482 : }
1483 178 : }
1484 :
1485 178 : void Metainference::get_sigma_mean(const unsigned iselect, const double fact, const double var_fact, const vector<double> &mean)
1486 : {
1487 178 : const double dnrep = static_cast<double>(nrep_);
1488 178 : const double ave_fact = 1.0/dnrep;
1489 :
1490 356 : vector<double> sigma_mean2_tmp(sigma_mean2_.size(), 0.);
1491 :
1492 178 : if(do_optsigmamean_>0) {
1493 : // remove first entry of the history vector
1494 0 : if(sigma_mean2_last_[iselect][0].size()==optsigmamean_stride_&&optsigmamean_stride_>0)
1495 0 : for(unsigned i=0; i<narg; ++i) sigma_mean2_last_[iselect][i].erase(sigma_mean2_last_[iselect][i].begin());
1496 : /* this is the current estimate of sigma mean for each argument
1497 : there is one of this per argument in any case because it is
1498 : the maximum among these to be used in case of GAUSS/OUTLIER */
1499 0 : vector<double> sigma_mean2_now(narg,0);
1500 0 : if(do_reweight_) {
1501 0 : if(master) {
1502 0 : for(unsigned i=0; i<narg; ++i) {
1503 0 : double tmp1 = (fact*getArgument(i)-ave_fact*mean[i])*(fact*getArgument(i)-ave_fact*mean[i]);
1504 0 : double tmp2 = -2.*mean[i]*(fact-ave_fact)*(fact*getArgument(i)-ave_fact*mean[i]);
1505 0 : sigma_mean2_now[i] = tmp1 + tmp2;
1506 : }
1507 0 : if(nrep_>1) multi_sim_comm.Sum(&sigma_mean2_now[0], narg);
1508 : }
1509 0 : comm.Sum(&sigma_mean2_now[0], narg);
1510 0 : for(unsigned i=0; i<narg; ++i) sigma_mean2_now[i] = dnrep/(dnrep-1.)*(sigma_mean2_now[i] + mean[i]*mean[i]*var_fact);
1511 : } else {
1512 0 : if(master) {
1513 0 : for(unsigned i=0; i<narg; ++i) {
1514 0 : double tmp = getArgument(i)-mean[i];
1515 0 : sigma_mean2_now[i] = fact*tmp*tmp;
1516 : }
1517 0 : if(nrep_>1) multi_sim_comm.Sum(&sigma_mean2_now[0], narg);
1518 : }
1519 0 : comm.Sum(&sigma_mean2_now[0], narg);
1520 0 : for(unsigned i=0; i<narg; ++i) sigma_mean2_now[i] /= dnrep;
1521 : }
1522 :
1523 : // add sigma_mean2 to history
1524 0 : if(optsigmamean_stride_>0) {
1525 0 : for(unsigned i=0; i<narg; ++i) sigma_mean2_last_[iselect][i].push_back(sigma_mean2_now[i]);
1526 : } else {
1527 0 : for(unsigned i=0; i<narg; ++i) if(sigma_mean2_now[i] > sigma_mean2_last_[iselect][i][0]) sigma_mean2_last_[iselect][i][0] = sigma_mean2_now[i];
1528 : }
1529 :
1530 0 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
1531 0 : for(unsigned i=0; i<narg; ++i) {
1532 : /* set to the maximum in history vector */
1533 0 : sigma_mean2_tmp[i] = *max_element(sigma_mean2_last_[iselect][i].begin(), sigma_mean2_last_[iselect][i].end());
1534 : /* the standard error of the mean */
1535 0 : valueSigmaMean[i]->set(sqrt(sigma_mean2_tmp[i]));
1536 0 : if(noise_type_==GENERIC) {
1537 0 : sigma_min_[i] = sqrt(sigma_mean2_tmp[i]);
1538 0 : if(sigma_[i] < sigma_min_[i]) sigma_[i] = sigma_min_[i];
1539 : }
1540 : }
1541 0 : } else if(noise_type_==GAUSS||noise_type_==OUTLIERS) {
1542 : // find maximum for each data point
1543 : vector <double> max_values;
1544 0 : for(unsigned i=0; i<narg; ++i) max_values.push_back(*max_element(sigma_mean2_last_[iselect][i].begin(), sigma_mean2_last_[iselect][i].end()));
1545 : // find maximum across data points
1546 0 : const double max_now = *max_element(max_values.begin(), max_values.end());
1547 : // set new value
1548 0 : sigma_mean2_tmp[0] = max_now;
1549 0 : valueSigmaMean[0]->set(sqrt(sigma_mean2_tmp[0]));
1550 : }
1551 : // endif sigma optimization
1552 : } else {
1553 178 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
1554 5640 : for(unsigned i=0; i<narg; ++i) {
1555 5528 : sigma_mean2_tmp[i] = sigma_mean2_last_[iselect][i][0];
1556 5528 : valueSigmaMean[i]->set(sqrt(sigma_mean2_tmp[i]));
1557 : }
1558 66 : } else if(noise_type_==GAUSS||noise_type_==OUTLIERS) {
1559 132 : sigma_mean2_tmp[0] = sigma_mean2_last_[iselect][0][0];
1560 132 : valueSigmaMean[0]->set(sqrt(sigma_mean2_tmp[0]));
1561 : }
1562 : }
1563 :
1564 178 : sigma_mean2_ = sigma_mean2_tmp;
1565 178 : }
1566 :
1567 178 : void Metainference::replica_averaging(const double fact, vector<double> &mean, vector<double> &dmean_b)
1568 : {
1569 178 : if(master) {
1570 3216 : for(unsigned i=0; i<narg; ++i) mean[i] = fact*getArgument(i);
1571 178 : if(nrep_>1) multi_sim_comm.Sum(&mean[0], narg);
1572 : }
1573 356 : comm.Sum(&mean[0], narg);
1574 : // set the derivative of the mean with respect to the bias
1575 9244 : for(unsigned i=0; i<narg; ++i) dmean_b[i] = fact/kbt_*(getArgument(i)-mean[i])/static_cast<double>(average_weights_stride_);
1576 :
1577 : // this is only for generic metainference
1578 178 : if(firstTime) {ftilde_ = mean; firstTime = false;}
1579 178 : }
1580 :
1581 0 : void Metainference::do_regression_zero(const vector<double> &mean)
1582 : {
1583 : // parameters[i] = scale_ * mean[i]: find scale_ with linear regression
1584 : double num = 0.0;
1585 : double den = 0.0;
1586 0 : for(unsigned i=0; i<parameters.size(); ++i) {
1587 0 : num += mean[i] * parameters[i];
1588 0 : den += mean[i] * mean[i];
1589 : }
1590 0 : if(den>0) {
1591 0 : scale_ = num / den;
1592 : } else {
1593 0 : scale_ = 1.0;
1594 : }
1595 0 : }
1596 :
1597 178 : void Metainference::calculate()
1598 : {
1599 : // get step
1600 178 : const long int step = getStep();
1601 :
1602 : unsigned iselect = 0;
1603 : // set the value of selector for REM-like stuff
1604 178 : if(selector_.length()>0) iselect = static_cast<unsigned>(plumed.passMap[selector_]);
1605 :
1606 178 : double fact = 0.0;
1607 178 : double var_fact = 0.0;
1608 : // get weights for ensemble average
1609 178 : get_weights(iselect, fact, var_fact);
1610 : // calculate the mean
1611 178 : vector<double> mean(narg,0);
1612 : // this is the derivative of the mean with respect to the argument
1613 178 : vector<double> dmean_x(narg,fact);
1614 : // this is the derivative of the mean with respect to the bias
1615 178 : vector<double> dmean_b(narg,0);
1616 : // calculate it
1617 178 : replica_averaging(fact, mean, dmean_b);
1618 : // calculate sigma mean
1619 178 : get_sigma_mean(iselect, fact, var_fact, mean);
1620 : // in case of regression with zero intercept, calculate scale
1621 178 : if(doregres_zero_ && step%nregres_zero_==0) do_regression_zero(mean);
1622 :
1623 : /* MONTE CARLO */
1624 178 : if(step%MCstride_==0&&!getExchangeStep()) doMonteCarlo(mean);
1625 :
1626 : // calculate bias and forces
1627 : double ene = 0;
1628 178 : switch(noise_type_) {
1629 12 : case GAUSS:
1630 12 : ene = getEnergyForceGJ(mean, dmean_x, dmean_b);
1631 : break;
1632 52 : case MGAUSS:
1633 52 : ene = getEnergyForceGJE(mean, dmean_x, dmean_b);
1634 : break;
1635 54 : case OUTLIERS:
1636 54 : ene = getEnergyForceSP(mean, dmean_x, dmean_b);
1637 : break;
1638 48 : case MOUTLIERS:
1639 48 : ene = getEnergyForceSPE(mean, dmean_x, dmean_b);
1640 : break;
1641 12 : case GENERIC:
1642 12 : ene = getEnergyForceMIGEN(mean, dmean_x, dmean_b);
1643 : break;
1644 : }
1645 :
1646 : // set value of the bias
1647 : setBias(ene);
1648 178 : }
1649 :
1650 19 : void Metainference::writeStatus()
1651 : {
1652 19 : sfile_.rewind();
1653 38 : sfile_.printField("time",getTimeStep()*getStep());
1654 : //nsel
1655 95 : for(unsigned i=0; i<sigma_mean2_last_.size(); i++) {
1656 : std::string msg_i,msg_j;
1657 19 : Tools::convert(i,msg_i);
1658 : vector <double> max_values;
1659 : //narg
1660 4849 : for(unsigned j=0; j<narg; ++j) {
1661 2415 : Tools::convert(j,msg_j);
1662 4830 : std::string msg = msg_i+"_"+msg_j;
1663 2415 : if(noise_type_==MGAUSS||noise_type_==MOUTLIERS||noise_type_==GENERIC) {
1664 7170 : sfile_.printField("sigmaMean_"+msg,sqrt(*max_element(sigma_mean2_last_[i][j].begin(), sigma_mean2_last_[i][j].end())));
1665 : } else {
1666 : // find maximum for each data point
1667 50 : max_values.push_back(*max_element(sigma_mean2_last_[i][j].begin(), sigma_mean2_last_[i][j].end()));
1668 : }
1669 : }
1670 19 : if(noise_type_==GAUSS||noise_type_==OUTLIERS) {
1671 : // find maximum across data points
1672 6 : const double max_now = sqrt(*max_element(max_values.begin(), max_values.end()));
1673 6 : Tools::convert(0,msg_j);
1674 12 : std::string msg = msg_i+"_"+msg_j;
1675 12 : sfile_.printField("sigmaMean_"+msg, max_now);
1676 : }
1677 : }
1678 7226 : for(unsigned i=0; i<sigma_.size(); ++i) {
1679 : std::string msg;
1680 2396 : Tools::convert(i,msg);
1681 4792 : sfile_.printField("sigma_"+msg,sigma_[i]);
1682 : }
1683 19 : if(noise_type_==GENERIC) {
1684 8 : for(unsigned i=0; i<ftilde_.size(); ++i) {
1685 : std::string msg;
1686 2 : Tools::convert(i,msg);
1687 4 : sfile_.printField("ftilde_"+msg,ftilde_[i]);
1688 : }
1689 : }
1690 38 : sfile_.printField("scale0_",scale_);
1691 38 : sfile_.printField("offset0_",offset_);
1692 95 : for(unsigned i=0; i<average_weights_.size(); i++) {
1693 : std::string msg_i;
1694 19 : Tools::convert(i,msg_i);
1695 57 : sfile_.printField("weight_"+msg_i,average_weights_[i][replica_]);
1696 : }
1697 19 : sfile_.printField();
1698 19 : sfile_.flush();
1699 19 : }
1700 :
1701 178 : void Metainference::update() {
1702 : // write status file
1703 178 : if(write_stride_>0&& (getStep()%write_stride_==0 || getCPT()) ) writeStatus();
1704 178 : }
1705 :
1706 : }
1707 5517 : }
1708 :
|
