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