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