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