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