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