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