Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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 "Bias.h"
23 : #include "ActionRegister.h"
24 : #include "core/PlumedMain.h"
25 : #include "core/Atoms.h"
26 : #include "core/Value.h"
27 : #include <cmath>
28 : #include <ctime>
29 :
30 : using namespace std;
31 :
32 : namespace PLMD {
33 : namespace bias {
34 :
35 : //+PLUMEDOC BIAS METAINFERENCE
36 : /*
37 : Calculate the Metainference Score for a set of back calculated experimental data.
38 :
39 :
40 : The back calculated data, that are expected to be averages over replicas (NR=1,2,..,N)
41 : The functional form of this bias can be chosen between three variants selected
42 : with NOISE=GAUSS,MGAUSS,OUTLIERS which correspond to modelling the noise for
43 : the arguments as a single gaussian common to all the data points, a gaussian per data
44 : point or a single long-tailed gaussian common to all the data points.
45 :
46 : As from Metainference theory there are two sigma values: SIGMA_MEAN represent the
47 : error of calculating an average quanity using a finite set of replica and should
48 : be set as small as possible following the guidelines for replica-averaged simulations
49 : in the framework of the Maximum Entropy Principle.
50 : SIGMA_BIAS is an uncertainty parameter, sampled by a MC algorithm in the bounded interval
51 : defined by SIGMA_MIN and SIGMA_MAX. The initial value is set at SIGMA0. The MC move is a
52 : random displacement of maximum value equal to DSIGMA.
53 :
54 : \par Examples
55 :
56 : In the following example we calculate a set of \ref RDC, take the replica-average of
57 : them and comparing them with a set of experimental values. RDCs are compared with
58 : the experimental data but for a multiplication factor SCALE that is also sampled by
59 : MC on-the-fly
60 :
61 : \verbatim
62 : RDC ...
63 : LABEL=rdc
64 : SCALE=0.0001
65 : GYROM=-72.5388
66 : ATOMS1=22,23
67 : ATOMS2=25,27
68 : ATOMS3=29,31
69 : ATOMS4=33,34
70 : ... RDC
71 :
72 : ardc: ENSEMBLE ARG=rdc.*
73 :
74 : METAINFERENCE ...
75 : ARG=ardc.*
76 : NOISETYPE=MGAUSS
77 : PARAMETERS=1.9190,2.9190,3.9190,4.9190
78 : SCALEDATA SCALE0=1 SCALE_MIN=0.00001 SCALE_MAX=3 DSCALE=0.00
79 : SIGMA0=0.01 SIGMA_MIN=0.00001 SIGMA_MAX=3 DSIGMA=0.00
80 : SIGMA_MEAN=0.001
81 : TEMP=300
82 : LABEL=spe
83 : ... METAINFERENCE
84 :
85 : PRINT ARG=spe.bias FILE=BIAS STRIDE=1
86 : \endverbatim
87 :
88 : in the following example instead of using one uncertainty parameter per data point we use
89 : a single uncertainty value in a long-tailed gaussian to take into account for outliers.
90 :
91 : \verbatim
92 : METAINFERENCE ...
93 : ARG=ardc.*
94 : NOISETYPE=OUTLIERS
95 : PARAMETERS=1.9190,2.9190,3.9190,4.9190
96 : SCALEDATA SCALE0=1 SCALE_MIN=0.00001 SCALE_MAX=3 DSCALE=0.00
97 : SIGMA0=0.01 SIGMA_MIN=0.00001 SIGMA_MAX=3 DSIGMA=0.00
98 : SIGMA_MEAN=0.001
99 : TEMP=300
100 : LABEL=spe
101 : ... METAINFERENCE
102 : \endverbatim
103 :
104 : (See also \ref RDC and \ref ENSEMBLE).
105 :
106 : */
107 : //+ENDPLUMEDOC
108 :
109 16 : class Metainference : public Bias
110 : {
111 : const double sqrt2_div_pi;
112 : const double sqrt2_pi;
113 : // experimental values
114 : vector<double> parameters;
115 : // noise type
116 : unsigned noise_type_;
117 : enum { GAUSS, MGAUSS, OUTLIERS };
118 : // scale is data scaling factor
119 : bool doscale_;
120 : double scale_;
121 : double scale_min_;
122 : double scale_max_;
123 : double Dscale_;
124 : // sigma is data uncertainty
125 : vector<double> sigma_;
126 : double sigma_min_;
127 : double sigma_max_;
128 : double Dsigma_;
129 : // sigma_mean is uncertainty in the mean estimate
130 : double sigma_mean_;
131 : // temperature in kbt
132 : double kbt_;
133 : // number of data points
134 : unsigned ndata_;
135 : // Monte Carlo stuff
136 : unsigned MCsteps_;
137 : unsigned MCstride_;
138 : unsigned MCaccept_;
139 : long int MCfirst_;
140 : // output
141 : Value* valueScale;
142 : Value* valueAccept;
143 : vector<Value*> valueSigma;
144 :
145 : unsigned nrep_;
146 : unsigned replica_;
147 :
148 : double getEnergySPE(const vector<double> &sigma, const double scale);
149 : double getEnergyGJE(const vector<double> &sigma, const double scale);
150 : void doMonteCarlo();
151 : double getEnergyForceSPE();
152 : double getEnergyForceGJE();
153 :
154 : public:
155 : explicit Metainference(const ActionOptions&);
156 : void calculate();
157 : static void registerKeywords(Keywords& keys);
158 : };
159 :
160 :
161 2531 : PLUMED_REGISTER_ACTION(Metainference,"METAINFERENCE")
162 :
163 9 : void Metainference::registerKeywords(Keywords& keys) {
164 9 : Bias::registerKeywords(keys);
165 9 : keys.use("ARG");
166 9 : keys.add("optional","PARARG","reference values for the experimental data, these can be provided as arguments without derivatives");
167 9 : keys.add("optional","PARAMETERS","reference values for the experimental data");
168 9 : keys.add("compulsory","NOISETYPE","functional form of the noise (GAUSS,MGAUSS,OUTLIERS)");
169 9 : keys.addFlag("SCALEDATA",false,"Set to TRUE if you want to sample a scaling factor common to all values and replicas");
170 9 : keys.add("compulsory","SCALE0","initial value of the uncertainty parameter");
171 9 : keys.add("compulsory","SCALE_MIN","minimum value of the uncertainty parameter");
172 9 : keys.add("compulsory","SCALE_MAX","maximum value of the uncertainty parameter");
173 9 : keys.add("compulsory","DSCALE","maximum MC move of the uncertainty parameter");
174 9 : keys.add("compulsory","SIGMA0","initial value of the uncertainty parameter");
175 9 : keys.add("compulsory","SIGMA_MIN","minimum value of the uncertainty parameter");
176 9 : keys.add("compulsory","SIGMA_MAX","maximum value of the uncertainty parameter");
177 9 : keys.add("compulsory","DSIGMA","maximum MC move of the uncertainty parameter");
178 9 : keys.add("compulsory","SIGMA_MEAN","starting value for the uncertainty in the mean estimate");
179 9 : keys.add("optional","TEMP","the system temperature - this is only needed if code doesnt' pass the temperature to plumed");
180 9 : keys.add("optional","MC_STEPS","number of MC steps");
181 9 : keys.add("optional","MC_STRIDE","MC stride");
182 9 : useCustomisableComponents(keys);
183 9 : keys.addOutputComponent("sigma", "default","uncertainty parameter");
184 9 : keys.addOutputComponent("scale", "default","scale parameter");
185 9 : keys.addOutputComponent("accept","default","MC acceptance");
186 9 : }
187 :
188 8 : Metainference::Metainference(const ActionOptions&ao):
189 : PLUMED_BIAS_INIT(ao),
190 : sqrt2_div_pi(0.45015815807855),
191 : sqrt2_pi(2.50662827463100050240),
192 : doscale_(false),
193 8 : ndata_(getNumberOfArguments()),
194 : MCsteps_(1),
195 : MCstride_(1),
196 : MCaccept_(0),
197 16 : MCfirst_(-1)
198 : {
199 8 : parseVector("PARAMETERS",parameters);
200 8 : if(parameters.size()!=static_cast<unsigned>(getNumberOfArguments())&&!parameters.empty())
201 0 : error("Size of PARAMETERS array should be either 0 or the same as of the number of arguments in ARG1");
202 :
203 8 : vector<Value*> arg2;
204 8 : parseArgumentList("PARARG",arg2);
205 :
206 8 : if(!arg2.empty()) {
207 0 : if(parameters.size()>0) error("It is not possible to use PARARG and PARAMETERS together");
208 0 : if(arg2.size()!=getNumberOfArguments()) error("Size of PARARG array should be the same as number for arguments in ARG");
209 0 : for(unsigned i=0; i<arg2.size(); i++) {
210 0 : parameters.push_back(arg2[i]->get());
211 0 : if(arg2[i]->hasDerivatives()==true) error("PARARG can only accept arguments without derivatives");
212 : }
213 : }
214 :
215 8 : if(parameters.size()!=getNumberOfArguments())
216 0 : error("PARARG or PARAMETERS arrays should include the same number of elements as the arguments in ARG");
217 :
218 16 : string stringa_noise;
219 8 : parse("NOISETYPE",stringa_noise);
220 8 : if(stringa_noise=="GAUSS") noise_type_ = GAUSS;
221 8 : else if(stringa_noise=="MGAUSS") noise_type_ = MGAUSS;
222 4 : else if(stringa_noise=="OUTLIERS") noise_type_ = OUTLIERS;
223 0 : else error("Unkwnow noise type");
224 :
225 8 : parseFlag("SCALEDATA", doscale_);
226 8 : if(doscale_) {
227 8 : parse("SCALE0",scale_);
228 8 : parse("SCALE_MIN",scale_min_);
229 8 : parse("SCALE_MAX",scale_max_);
230 8 : parse("DSCALE",Dscale_);
231 : } else {
232 0 : scale_=1.0;
233 : }
234 :
235 16 : vector<double> readsigma;
236 8 : parseVector("SIGMA0",readsigma);
237 8 : if(noise_type_!=MGAUSS&&readsigma.size()>1) error("If you want to use more than one sigma you should use NOISETYPE=MGAUSS");
238 8 : if(noise_type_==MGAUSS) {
239 4 : if(readsigma.size()==getNumberOfArguments()) {
240 0 : sigma_.resize(getNumberOfArguments());
241 0 : sigma_=readsigma;
242 4 : } else if(readsigma.size()==1) {
243 4 : sigma_.resize(getNumberOfArguments(),readsigma[0]);
244 : } else {
245 0 : error("SIGMA0 can accept either one single value or as many values as the number of arguments (with NOISETYPE=MGAUSS)");
246 : }
247 4 : } else sigma_.resize(1, readsigma[0]);
248 :
249 8 : parse("SIGMA_MIN",sigma_min_);
250 8 : parse("SIGMA_MAX",sigma_max_);
251 8 : parse("DSIGMA",Dsigma_);
252 8 : parse("MC_STEPS",MCsteps_);
253 8 : parse("MC_STRIDE",MCstride_);
254 : // get temperature
255 8 : double temp=0.0;
256 8 : parse("TEMP",temp);
257 8 : parse("SIGMA_MEAN",sigma_mean_);
258 :
259 8 : checkRead();
260 :
261 8 : if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
262 0 : else kbt_=plumed.getAtoms().getKbT();
263 :
264 : // get number of replicas
265 8 : if(comm.Get_rank()==0) {
266 4 : nrep_ = multi_sim_comm.Get_size();
267 4 : replica_ = multi_sim_comm.Get_rank();
268 : } else {
269 4 : nrep_ = 0;
270 4 : replica_ = 0;
271 : }
272 8 : comm.Sum(&nrep_,1);
273 8 : comm.Sum(&replica_,1);
274 :
275 : // divide sigma_mean by the square root of the number of replicas
276 8 : sigma_mean_ /= sqrt(static_cast<double>(nrep_));
277 :
278 : // adjust for multiple-time steps
279 8 : MCstride_ *= getStride();
280 :
281 8 : switch(noise_type_) {
282 : case GAUSS:
283 0 : log.printf(" with gaussian noise and a single noise parameter for all the data\n");
284 0 : break;
285 : case MGAUSS:
286 4 : log.printf(" with gaussian noise and a noise parameter for each data point\n");
287 4 : break;
288 : case OUTLIERS:
289 4 : log.printf(" with long tailed gaussian noise and a single noise parameter for all the data\n");
290 4 : break;
291 : }
292 :
293 8 : if(doscale_) {
294 8 : log.printf(" sampling a common scaling factor with:\n");
295 8 : log.printf(" initial scale parameter %f\n",scale_);
296 8 : log.printf(" minimum scale parameter %f\n",scale_min_);
297 8 : log.printf(" maximum scale parameter %f\n",scale_max_);
298 8 : log.printf(" maximum MC move of scale parameter %f\n",Dscale_);
299 : }
300 :
301 8 : if(readsigma.size()==1) log.printf(" initial data uncertainty %f\n",sigma_[0]);
302 : else {
303 0 : log.printf(" initial data uncertainties");
304 0 : for(unsigned i=0; i<sigma_.size(); ++i) log.printf(" %f", sigma_[i]);
305 0 : log.printf("\n");
306 : }
307 8 : log.printf(" minimum data uncertainty %f\n",sigma_min_);
308 8 : log.printf(" maximum data uncertainty %f\n",sigma_max_);
309 8 : log.printf(" maximum MC move of data uncertainty %f\n",Dsigma_);
310 8 : log.printf(" uncertainty in the mean estimate %f\n",sigma_mean_);
311 8 : log.printf(" temperature of the system %f\n",kbt_);
312 8 : log.printf(" number of experimental data points %u\n",getNumberOfArguments());
313 8 : log.printf(" number of replicas %u\n",nrep_);
314 8 : log.printf(" MC steps %u\n",MCsteps_);
315 8 : log.printf(" MC stride %u\n",MCstride_);
316 :
317 8 : if(doscale_) {
318 8 : addComponent("scale");
319 8 : componentIsNotPeriodic("scale");
320 8 : valueScale=getPntrToComponent("scale");
321 : }
322 8 : addComponent("accept");
323 8 : componentIsNotPeriodic("accept");
324 8 : valueAccept=getPntrToComponent("accept");
325 :
326 8 : if(noise_type_==MGAUSS) {
327 20 : for(unsigned i=0; i<sigma_.size(); ++i) {
328 16 : std::string num; Tools::convert(i,num);
329 16 : addComponent("sigma_"+num); componentIsNotPeriodic("sigma_"+num);
330 16 : valueSigma.push_back(getPntrToComponent("sigma_"+num));
331 16 : }
332 : } else {
333 4 : addComponent("sigma"); componentIsNotPeriodic("sigma");
334 4 : valueSigma.push_back(getPntrToComponent("sigma"));
335 : }
336 : // initialize random seed
337 : unsigned iseed;
338 8 : if(comm.Get_rank()==0) iseed = time(NULL)+replica_;
339 4 : else iseed = 0;
340 8 : comm.Sum(&iseed, 1);
341 8 : srand(iseed);
342 :
343 16 : log<<" Bibliography "<<plumed.cite("Bonomi, Camilloni, Cavalli, Vendruscolo, Sci. Adv. 2, e150117 (2016)");
344 8 : }
345 :
346 96 : double Metainference::getEnergySPE(const vector<double> &sigma, const double scale) {
347 : // calculate effective sigma
348 96 : const double smean2 = sigma_mean_*sigma_mean_;
349 96 : const double s = sqrt( sigma[0]*sigma[0] + smean2 );
350 : // cycle on arguments
351 96 : double ene = 0.0;
352 480 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
353 384 : const double dev = scale*getArgument(i)-parameters[i];
354 : // argument
355 384 : const double a2 = 0.5*dev*dev + s*s;
356 : // increment energy
357 384 : ene += std::log( 2.0 * a2 / ( 1.0 - exp(- a2 / smean2) ) );
358 : }
359 : // add normalization and Jeffrey's prior
360 96 : ene += std::log(s) - static_cast<double>(ndata_)*std::log(sqrt2_div_pi*s);
361 96 : return kbt_ * ene;
362 : }
363 :
364 96 : double Metainference::getEnergyGJE(const vector<double> &sigma, const double scale) {
365 : // cycle on arguments
366 96 : double ene = 0.0;
367 96 : const double smean2 = sigma_mean_*sigma_mean_;
368 96 : double ss = sigma[0]*sigma[0] + smean2;
369 480 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
370 384 : if(noise_type_==MGAUSS) {
371 384 : ss = sigma[i]*sigma[i] + smean2;
372 : // add Jeffrey's prior - one per sigma
373 384 : ene += 0.5*std::log(ss);
374 : }
375 384 : const double dev = scale*getArgument(i)-parameters[i];
376 384 : ene += 0.5*dev*dev/ss + 0.5*std::log(ss*sqrt2_pi);
377 : }
378 : // add Jeffrey's prior in case one sigma for all data points
379 96 : if(noise_type_==GAUSS) ene += 0.5*std::log(ss);
380 96 : return kbt_ * ene;
381 : }
382 :
383 96 : void Metainference::doMonteCarlo() {
384 : double old_energy;
385 96 : switch(noise_type_) {
386 : case GAUSS:
387 : case MGAUSS:
388 48 : old_energy = getEnergyGJE(sigma_,scale_);
389 48 : break;
390 : case OUTLIERS:
391 48 : old_energy = getEnergySPE(sigma_,scale_);
392 48 : break;
393 : }
394 :
395 : // cycle on MC steps
396 192 : for(unsigned i=0; i<MCsteps_; ++i) {
397 :
398 : // propose move for scale
399 96 : double new_scale = scale_;
400 96 : if(doscale_) {
401 96 : const double r1 = static_cast<double>(rand()) / RAND_MAX;
402 96 : const double ds1 = -Dscale_ + r1 * 2.0 * Dscale_;
403 96 : new_scale += ds1;
404 : // check boundaries
405 96 : if(new_scale > scale_max_) {new_scale = 2.0 * scale_max_ - new_scale;}
406 96 : if(new_scale < scale_min_) {new_scale = 2.0 * scale_min_ - new_scale;}
407 : // the scaling factor should be the same for all the replicas
408 96 : if(comm.Get_rank()==0) multi_sim_comm.Bcast(new_scale,0);
409 96 : comm.Bcast(new_scale,0);
410 : }
411 :
412 : // propose move for sigma
413 96 : vector<double> new_sigma(sigma_.size());
414 336 : for(unsigned j=0; j<sigma_.size(); j++) {
415 240 : const double r2 = static_cast<double>(rand()) / RAND_MAX;
416 240 : const double ds2 = -Dsigma_ + r2 * 2.0 * Dsigma_;
417 240 : new_sigma[j] = sigma_[j] + ds2;
418 : // check boundaries
419 240 : if(new_sigma[j] > sigma_max_) {new_sigma[j] = 2.0 * sigma_max_ - new_sigma[j];}
420 240 : if(new_sigma[j] < sigma_min_) {new_sigma[j] = 2.0 * sigma_min_ - new_sigma[j];}
421 : }
422 :
423 : // calculate new energy
424 96 : double new_energy=0;
425 96 : switch(noise_type_) {
426 : case GAUSS:
427 : case MGAUSS:
428 48 : new_energy = getEnergyGJE(new_sigma,new_scale);
429 48 : break;
430 : case OUTLIERS:
431 48 : new_energy = getEnergySPE(new_sigma,new_scale);
432 48 : break;
433 : }
434 : // accept or reject
435 96 : const double delta = ( new_energy - old_energy ) / kbt_;
436 : // if delta is negative always accept move
437 96 : if( delta <= 0.0 ) {
438 96 : old_energy = new_energy;
439 96 : scale_ = new_scale;
440 96 : sigma_ = new_sigma;
441 96 : MCaccept_++;
442 : // otherwise extract random number
443 : } else {
444 0 : const double s = static_cast<double>(rand()) / RAND_MAX;
445 0 : if( s < exp(-delta) ) {
446 0 : old_energy = new_energy;
447 0 : scale_ = new_scale;
448 0 : sigma_ = new_sigma;
449 0 : MCaccept_++;
450 : }
451 : }
452 :
453 96 : if(doscale_) {
454 : // the scaling factor should be the same for all the replicas
455 96 : if(comm.Get_rank()==0) multi_sim_comm.Bcast(scale_,0);
456 96 : comm.Bcast(scale_,0);
457 : }
458 96 : }
459 : /* save the result of the sampling */
460 96 : if(doscale_) valueScale->set(scale_);
461 96 : for(unsigned i=0; i<sigma_.size(); i++) valueSigma[i]->set(sigma_[i]);
462 96 : }
463 :
464 48 : double Metainference::getEnergyForceSPE()
465 : {
466 48 : double ene = 0.0;
467 48 : const unsigned narg=getNumberOfArguments();
468 :
469 48 : const double smean2 = sigma_mean_*sigma_mean_;
470 48 : const double s = sqrt( sigma_[0]*sigma_[0] + smean2 );
471 48 : vector<double> f(narg,0);
472 :
473 48 : if(comm.Get_rank()==0) {
474 120 : for(unsigned i=0; i<narg; ++i) {
475 96 : const double dev = scale_*getArgument(i)-parameters[i];
476 96 : const double a2 = 0.5*dev*dev + s*s;
477 96 : const double t = exp(-a2/smean2);
478 96 : const double dt = 1./t;
479 96 : const double it = 1./(1.-t);
480 96 : const double dit = 1./(1.-dt);
481 96 : ene += std::log(2.*a2*it);
482 96 : f[i] = -scale_*dev*(dit/smean2 + 1./a2);
483 : }
484 : // collect contribution to forces and energy from other replicas
485 24 : multi_sim_comm.Sum(&f[0],narg);
486 24 : multi_sim_comm.Sum(&ene,1);
487 : // add normalizations and priors of local replica
488 24 : ene += std::log(s) - static_cast<double>(ndata_)*std::log(sqrt2_div_pi*s);
489 : }
490 : // intra-replica summation
491 48 : comm.Sum(&f[0],narg);
492 48 : comm.Sum(&ene,1);
493 :
494 48 : for(unsigned i=0; i<narg; ++i) setOutputForce(i, kbt_ * f[i]);
495 48 : return ene;
496 : }
497 :
498 48 : double Metainference::getEnergyForceGJE()
499 : {
500 48 : double ene = 0.0;
501 :
502 48 : const unsigned ssize = sigma_.size();
503 48 : vector<double> ss(ssize);
504 96 : vector<double> inv_s2(ssize, 0.);
505 48 : const double smean2 = sigma_mean_*sigma_mean_;
506 :
507 240 : for(unsigned i=0; i<ssize; ++i) {
508 192 : ss[i] = sigma_[i]*sigma_[i] + smean2;
509 192 : if(comm.Get_rank()==0) inv_s2[i] = 1.0/ss[i];
510 : }
511 :
512 48 : if(comm.Get_rank()==0) multi_sim_comm.Sum(&inv_s2[0],ssize);
513 48 : comm.Sum(&inv_s2[0],ssize);
514 :
515 48 : const unsigned narg=getNumberOfArguments();
516 240 : for(unsigned i=0; i<narg; ++i) {
517 192 : const double dev = scale_*getArgument(i)-parameters[i];
518 192 : unsigned sel_sigma=0;
519 192 : if(noise_type_==MGAUSS) {
520 192 : sel_sigma=i;
521 : // add Jeffrey's prior - one per sigma
522 192 : ene += 0.5*std::log(ss[sel_sigma]);
523 : }
524 192 : ene += 0.5*dev*dev*inv_s2[sel_sigma] + 0.5*std::log(ss[sel_sigma]*sqrt2_pi);
525 192 : setOutputForce(i, -kbt_*dev*scale_*inv_s2[sel_sigma]);
526 : }
527 : // add Jeffrey's prior in case one sigma for all data points
528 48 : if(noise_type_==GAUSS) ene += 0.5*std::log(ss[0]);
529 96 : return ene;
530 : }
531 :
532 96 : void Metainference::calculate() {
533 : /* MONTE CARLO */
534 96 : const long int step = getStep();
535 96 : if(step%MCstride_==0&&!getExchangeStep()) doMonteCarlo();
536 : // this is needed when restarting simulations
537 96 : if(MCfirst_==-1) MCfirst_=step;
538 96 : const double MCtrials = floor(static_cast<double>(step-MCfirst_) / static_cast<double>(MCstride_))+1.0;
539 96 : const double accept = static_cast<double>(MCaccept_) / static_cast<double>(MCsteps_) / MCtrials;
540 96 : valueAccept->set(accept);
541 :
542 : // calculate bias and forces
543 96 : double ene = 0;
544 96 : switch(noise_type_) {
545 : case GAUSS:
546 : case MGAUSS:
547 48 : ene = getEnergyForceGJE();
548 48 : break;
549 : case OUTLIERS:
550 48 : ene = getEnergyForceSPE();
551 48 : break;
552 : }
553 : // set value of the bias
554 96 : setBias(kbt_*ene);
555 96 : }
556 :
557 : }
558 2523 : }
559 :
560 :
|
