From c9a71cfca3486834e548b7b88f64c8a7972c7f01 Mon Sep 17 00:00:00 2001
From: Carlo Camilloni <carlo.camilloni@gmail.com>
Date: Wed, 18 May 2016 16:11:37 +0200
Subject: [PATCH] METAINFERENCE: merged GJE and SPE in a new BIAS the old one
are still there for reference
---
regtest/.gitignore | 2 +-
.../rt-bayesgje/BIAS.0.reference | 26 +-
.../rt-bayesgje/BIAS.1.reference | 26 +-
.../rt-bayesgje/forces.0.reference | 24 +-
.../rt-bayesgje/forces.1.reference | 24 +-
regtest/metainference/rt-bayesgje/plumed.dat | 9 +-
regtest/metainference/rt-bayesspe/plumed.dat | 5 +-
src/bias/Metainference.cpp | 490 ++++++++++++++++++
8 files changed, 549 insertions(+), 57 deletions(-)
create mode 100644 src/bias/Metainference.cpp
diff --git a/regtest/.gitignore b/regtest/.gitignore
index d7b7fd89d..5874524a8 100644
--- a/regtest/.gitignore
+++ b/regtest/.gitignore
@@ -15,7 +15,7 @@
!/scripts
!/manyrestraints
!/molfile_plugin
-!/bayesbias
+!/metainference
# These files we just want to ignore completely
tmp
report.txt
diff --git a/regtest/metainference/rt-bayesgje/BIAS.0.reference b/regtest/metainference/rt-bayesgje/BIAS.0.reference
index b8a7514ab..bf05ee067 100644
--- a/regtest/metainference/rt-bayesgje/BIAS.0.reference
+++ b/regtest/metainference/rt-bayesgje/BIAS.0.reference
@@ -1,13 +1,13 @@
-#! FIELDS time spe.bias spe.scale spe.accept spe.sigma_0 spe.kappa_0 spe.sigma_1 spe.kappa_1 spe.sigma_2 spe.kappa_2 spe.sigma_3 spe.kappa_3
- 0.000000 1237852.869575 1.000000 1.000000 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.050000 1282471.107299 1.000000 0.181818 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.100000 1315165.380004 1.000000 0.142857 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.150000 1333063.481567 1.000000 0.129032 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.200000 1351555.666802 1.000000 0.121951 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.250000 1371772.655731 1.000000 0.117647 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.300000 1293440.279295 1.000000 0.114754 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.350000 1292004.343160 1.000000 0.112676 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.400000 1310371.270133 1.000000 0.111111 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.450000 1243985.739979 1.000000 0.109890 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.500000 1266194.635296 1.000000 0.108911 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.550000 1247194.127448 1.000000 0.108108 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
+#! FIELDS time spe.bias spe.scale spe.accept spe.sigma_0 spe.sigma_1 spe.sigma_2 spe.sigma_3
+ 0.000000 1244011.703127 1.000000 1.000000 0.010000 0.010000 0.010000 0.010000
+ 0.050000 1288851.922133 1.000000 0.090909 0.010000 0.010000 0.010000 0.010000
+ 0.100000 1321708.852911 1.000000 0.047619 0.010000 0.010000 0.010000 0.010000
+ 0.150000 1339695.999756 1.000000 0.032258 0.010000 0.010000 0.010000 0.010000
+ 0.200000 1358280.185913 1.000000 0.024390 0.010000 0.010000 0.010000 0.010000
+ 0.250000 1378597.756876 1.000000 0.019608 0.010000 0.010000 0.010000 0.010000
+ 0.300000 1299875.667124 1.000000 0.016393 0.010000 0.010000 0.010000 0.010000
+ 0.350000 1298432.587028 1.000000 0.014085 0.010000 0.010000 0.010000 0.010000
+ 0.400000 1316890.891747 1.000000 0.012346 0.010000 0.010000 0.010000 0.010000
+ 0.450000 1250175.085324 1.000000 0.010989 0.010000 0.010000 0.010000 0.010000
+ 0.500000 1272494.472658 1.000000 0.009901 0.010000 0.010000 0.010000 0.010000
+ 0.550000 1253399.434919 1.000000 0.009009 0.010000 0.010000 0.010000 0.010000
diff --git a/regtest/metainference/rt-bayesgje/BIAS.1.reference b/regtest/metainference/rt-bayesgje/BIAS.1.reference
index b8a7514ab..bf05ee067 100644
--- a/regtest/metainference/rt-bayesgje/BIAS.1.reference
+++ b/regtest/metainference/rt-bayesgje/BIAS.1.reference
@@ -1,13 +1,13 @@
-#! FIELDS time spe.bias spe.scale spe.accept spe.sigma_0 spe.kappa_0 spe.sigma_1 spe.kappa_1 spe.sigma_2 spe.kappa_2 spe.sigma_3 spe.kappa_3
- 0.000000 1237852.869575 1.000000 1.000000 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.050000 1282471.107299 1.000000 0.181818 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.100000 1315165.380004 1.000000 0.142857 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.150000 1333063.481567 1.000000 0.129032 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.200000 1351555.666802 1.000000 0.121951 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.250000 1371772.655731 1.000000 0.117647 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.300000 1293440.279295 1.000000 0.114754 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.350000 1292004.343160 1.000000 0.112676 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.400000 1310371.270133 1.000000 0.111111 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.450000 1243985.739979 1.000000 0.109890 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.500000 1266194.635296 1.000000 0.108911 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
- 0.550000 1247194.127448 1.000000 0.108108 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158 0.010000 49392.844158
+#! FIELDS time spe.bias spe.scale spe.accept spe.sigma_0 spe.sigma_1 spe.sigma_2 spe.sigma_3
+ 0.000000 1244011.703127 1.000000 1.000000 0.010000 0.010000 0.010000 0.010000
+ 0.050000 1288851.922133 1.000000 0.090909 0.010000 0.010000 0.010000 0.010000
+ 0.100000 1321708.852911 1.000000 0.047619 0.010000 0.010000 0.010000 0.010000
+ 0.150000 1339695.999756 1.000000 0.032258 0.010000 0.010000 0.010000 0.010000
+ 0.200000 1358280.185913 1.000000 0.024390 0.010000 0.010000 0.010000 0.010000
+ 0.250000 1378597.756876 1.000000 0.019608 0.010000 0.010000 0.010000 0.010000
+ 0.300000 1299875.667124 1.000000 0.016393 0.010000 0.010000 0.010000 0.010000
+ 0.350000 1298432.587028 1.000000 0.014085 0.010000 0.010000 0.010000 0.010000
+ 0.400000 1316890.891747 1.000000 0.012346 0.010000 0.010000 0.010000 0.010000
+ 0.450000 1250175.085324 1.000000 0.010989 0.010000 0.010000 0.010000 0.010000
+ 0.500000 1272494.472658 1.000000 0.009901 0.010000 0.010000 0.010000 0.010000
+ 0.550000 1253399.434919 1.000000 0.009009 0.010000 0.010000 0.010000 0.010000
diff --git a/regtest/metainference/rt-bayesgje/forces.0.reference b/regtest/metainference/rt-bayesgje/forces.0.reference
index faefc5998..0ed61d28d 100644
--- a/regtest/metainference/rt-bayesgje/forces.0.reference
+++ b/regtest/metainference/rt-bayesgje/forces.0.reference
@@ -1,13 +1,13 @@
#! FIELDS time ardc.rdc.rdc_0 ardc.rdc.rdc_1 ardc.rdc.rdc_2 ardc.rdc.rdc_3
- 0.000000 100642.38088002 141989.23823686 180293.36768282 243915.51628401
- 0.050000 96031.98232030 157899.37479933 182544.22046024 243354.07707703
- 0.100000 86716.33822558 159398.77609836 187600.18135685 248608.44413448
- 0.150000 81737.66067248 165952.55849053 203282.42805091 236961.41787945
- 0.200000 84584.72558079 163371.00532409 207785.88775295 237703.35057586
- 0.250000 93071.99457441 165199.43919893 194520.13904297 248451.87259207
- 0.300000 88502.40103501 166208.29851741 179909.50684797 244859.64655378
- 0.350000 88835.19479396 143704.44676296 195579.89949173 246668.89502228
- 0.400000 98275.89353361 137347.51224637 199311.50046416 247399.52636605
- 0.450000 98743.90812140 135101.70403599 196684.05229235 237083.74826703
- 0.500000 102074.94767658 127111.08569720 203939.31073060 238583.95889519
- 0.550000 102375.94644906 127568.39283256 186229.01549716 248550.43197152
+ 0.000000 101143.08924261 142695.65235745 181190.34961159 245129.02631528
+ 0.050000 96509.75337662 158684.94382818 183452.40066154 244564.79387841
+ 0.100000 87147.76279387 160191.80483517 188533.51559246 249845.30206550
+ 0.150000 82144.31570070 166778.19310989 204293.78341435 238140.33040621
+ 0.200000 85005.54511104 164183.79639535 208819.64838854 238885.95431007
+ 0.250000 93535.03932353 166021.32695614 195487.90092875 249687.95156019
+ 0.300000 88942.71148792 167035.20547521 180804.57902134 246077.85375057
+ 0.350000 89277.16093721 144419.39425930 196552.93381756 247896.10345523
+ 0.400000 98764.82832731 138030.83320282 200303.09996895 248630.36978081
+ 0.450000 99235.17134588 135773.85181726 197662.57991569 238263.26940268
+ 0.500000 102582.78323716 127743.47915838 204953.93416707 239770.94376532
+ 0.550000 102885.27951597 128203.06145362 187155.52801207 249787.00128481
diff --git a/regtest/metainference/rt-bayesgje/forces.1.reference b/regtest/metainference/rt-bayesgje/forces.1.reference
index faefc5998..0ed61d28d 100644
--- a/regtest/metainference/rt-bayesgje/forces.1.reference
+++ b/regtest/metainference/rt-bayesgje/forces.1.reference
@@ -1,13 +1,13 @@
#! FIELDS time ardc.rdc.rdc_0 ardc.rdc.rdc_1 ardc.rdc.rdc_2 ardc.rdc.rdc_3
- 0.000000 100642.38088002 141989.23823686 180293.36768282 243915.51628401
- 0.050000 96031.98232030 157899.37479933 182544.22046024 243354.07707703
- 0.100000 86716.33822558 159398.77609836 187600.18135685 248608.44413448
- 0.150000 81737.66067248 165952.55849053 203282.42805091 236961.41787945
- 0.200000 84584.72558079 163371.00532409 207785.88775295 237703.35057586
- 0.250000 93071.99457441 165199.43919893 194520.13904297 248451.87259207
- 0.300000 88502.40103501 166208.29851741 179909.50684797 244859.64655378
- 0.350000 88835.19479396 143704.44676296 195579.89949173 246668.89502228
- 0.400000 98275.89353361 137347.51224637 199311.50046416 247399.52636605
- 0.450000 98743.90812140 135101.70403599 196684.05229235 237083.74826703
- 0.500000 102074.94767658 127111.08569720 203939.31073060 238583.95889519
- 0.550000 102375.94644906 127568.39283256 186229.01549716 248550.43197152
+ 0.000000 101143.08924261 142695.65235745 181190.34961159 245129.02631528
+ 0.050000 96509.75337662 158684.94382818 183452.40066154 244564.79387841
+ 0.100000 87147.76279387 160191.80483517 188533.51559246 249845.30206550
+ 0.150000 82144.31570070 166778.19310989 204293.78341435 238140.33040621
+ 0.200000 85005.54511104 164183.79639535 208819.64838854 238885.95431007
+ 0.250000 93535.03932353 166021.32695614 195487.90092875 249687.95156019
+ 0.300000 88942.71148792 167035.20547521 180804.57902134 246077.85375057
+ 0.350000 89277.16093721 144419.39425930 196552.93381756 247896.10345523
+ 0.400000 98764.82832731 138030.83320282 200303.09996895 248630.36978081
+ 0.450000 99235.17134588 135773.85181726 197662.57991569 238263.26940268
+ 0.500000 102582.78323716 127743.47915838 204953.93416707 239770.94376532
+ 0.550000 102885.27951597 128203.06145362 187155.52801207 249787.00128481
diff --git a/regtest/metainference/rt-bayesgje/plumed.dat b/regtest/metainference/rt-bayesgje/plumed.dat
index c23798c09..70072bca8 100644
--- a/regtest/metainference/rt-bayesgje/plumed.dat
+++ b/regtest/metainference/rt-bayesgje/plumed.dat
@@ -11,15 +11,16 @@ ATOMS4=33,34
ardc: ENSEMBLE ARG=rdc.*
stat: STATS ARG=ardc.* PARAMETERS=1.9190,2.9190,3.9190,4.9190 SQDEV
-BAYESIANGJE ...
-ARG=ardc.*
+METAINFERENCE ...
+ARG=ardc.*
+NOISETYPE=MGAUSS
PARAMETERS=1.9190,2.9190,3.9190,4.9190
SCALEDATA SCALE0=1 SCALE_MIN=0.00001 SCALE_MAX=3 DSCALE=0.00
-MULTISIGMA SIGMA0=0.01 SIGMA_MIN=0.00001 SIGMA_MAX=3 DSIGMA=0.00
+SIGMA0=0.01 SIGMA_MIN=0.00001 SIGMA_MAX=3 DSIGMA=0.00
SIGMA_MEAN=0.001
TEMP=300
LABEL=spe
-... BAYESIANGJE
+... METAINFERENCE
DUMPFORCES ARG=ardc.* FILE=forces FMT=%15.8f
diff --git a/regtest/metainference/rt-bayesspe/plumed.dat b/regtest/metainference/rt-bayesspe/plumed.dat
index ac186bd9d..a070ff1ab 100644
--- a/regtest/metainference/rt-bayesspe/plumed.dat
+++ b/regtest/metainference/rt-bayesspe/plumed.dat
@@ -11,15 +11,16 @@ ATOMS4=33,34
ardc: ENSEMBLE ARG=rdc.*
stat: STATS ARG=ardc.* PARAMETERS=1.9190,2.9190,3.9190,4.9190 SQDEV
-BAYESIANSPE ...
+METAINFERENCE ...
ARG=ardc.*
PARAMETERS=1.9190,2.9190,3.9190,4.9190
+NOISETYPE=LTAIL
SCALEDATA SCALE0=1 SCALE_MIN=0.00001 SCALE_MAX=3 DSCALE=0.00
SIGMA0=0.01 SIGMA_MIN=0.00001 SIGMA_MAX=3 DSIGMA=0.00
SIGMA_MEAN=0.001
TEMP=300
LABEL=spe
-... BAYESIANSPE
+... METAINFERENCE
DUMPFORCES ARG=ardc.* FILE=forces
diff --git a/src/bias/Metainference.cpp b/src/bias/Metainference.cpp
new file mode 100644
index 000000000..911d07a50
--- /dev/null
+++ b/src/bias/Metainference.cpp
@@ -0,0 +1,490 @@
+/* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+ Copyright (c) 2013 The plumed team
+ (see the PEOPLE file at the root of the distribution for a list of names)
+
+ See http://www.plumed-code.org for more information.
+
+ This file is part of plumed, version 2.
+
+ plumed is free software: you can redistribute it and/or modify
+ it under the terms of the GNU Lesser General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
+
+ plumed is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU Lesser General Public License for more details.
+
+ You should have received a copy of the GNU Lesser General Public License
+ along with plumed. If not, see <http://www.gnu.org/licenses/>.
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
+#include "Bias.h"
+#include "ActionRegister.h"
+#include "core/PlumedMain.h"
+#include "core/Atoms.h"
+#include "core/Value.h"
+#include <cmath>
+
+using namespace std;
+
+namespace PLMD{
+namespace bias{
+
+//+PLUMEDOC BIAS METAINFERENCE
+/*
+Calculate the Metainference Score for a set of experimental Data.
+
+*/
+//+ENDPLUMEDOC
+
+class Metainference : public Bias
+{
+ const double sqrt2pi;
+ const double sqrt2_div_pi;
+ // experimental values
+ vector<double> parameters;
+ // noise type
+ unsigned noise_type_;
+ enum { GAUSS, MGAUSS, LTAIL};
+ // scale is data scaling factor
+ bool doscale_;
+ double scale_;
+ double scale_min_;
+ double scale_max_;
+ double Dscale_;
+ // sigma is data uncertainty
+ vector<double> sigma_;
+ double sigma_min_;
+ double sigma_max_;
+ double Dsigma_;
+ // sigma_mean is uncertainty in the mean estimate
+ double sigma_mean_;
+ // temperature in kbt
+ double kbt_;
+ // number of data points
+ unsigned ndata_;
+ // Monte Carlo stuff
+ double old_energy;
+ unsigned MCsteps_;
+ unsigned MCstride_;
+ unsigned MCaccept_;
+ long int MCfirst_;
+ // output
+ Value* valueBias;
+ Value* valueScale;
+ Value* valueAccept;
+ vector<Value*> valueSigma;
+
+ unsigned nrep_;
+ unsigned replica_;
+
+ double getEnergySPE(const double sigma, const double scale);
+ double getEnergyGJE(const vector<double> &sigma, const double scale);
+ void doMonteCarlo();
+ double getEnergyForceSPE();
+ double getEnergyForceGJE();
+
+public:
+ explicit Metainference(const ActionOptions&);
+ void calculate();
+ static void registerKeywords(Keywords& keys);
+};
+
+
+PLUMED_REGISTER_ACTION(Metainference,"METAINFERENCE")
+
+void Metainference::registerKeywords(Keywords& keys){
+ Bias::registerKeywords(keys);
+ keys.use("ARG");
+ keys.add("optional","PARARG","the input for this action is the scalar output from other actions without derivatives.");
+ keys.add("optional","PARAMETERS","the parameters of the arguments in your function");
+ keys.add("compulsory","NOISETYPE","functional form of the noise (GAUSS,MGAUSS,LTAIL)");
+ keys.addFlag("SCALEDATA",false,"Set to TRUE if you want to sample a scaling factor common to all values and replicas.");
+ keys.addFlag("OPTSIGMAMEAN", false, "Set to minimize sigma_mean on the fly");
+ keys.add("compulsory","SCALE0","initial value of the uncertainty parameter");
+ keys.add("compulsory","SCALE_MIN","minimum value of the uncertainty parameter");
+ keys.add("compulsory","SCALE_MAX","maximum value of the uncertainty parameter");
+ keys.add("compulsory","DSCALE","maximum MC move of the uncertainty parameter");
+ keys.add("compulsory","SIGMA0","initial value of the uncertainty parameter");
+ keys.add("compulsory","SIGMA_MIN","minimum value of the uncertainty parameter");
+ keys.add("compulsory","SIGMA_MAX","maximum value of the uncertainty parameter");
+ keys.add("compulsory","DSIGMA","maximum MC move of the uncertainty parameter");
+ keys.add("compulsory","SIGMA_MEAN","starting value for the uncertainty in the mean estimate");
+ keys.add("optional","TEMP","the system temperature - this is only needed if code doesnt' pass the temperature to plumed");
+ keys.add("optional","MC_STEPS","number of MC steps");
+ keys.add("optional","MC_STRIDE","MC stride");
+ componentsAreNotOptional(keys);
+ useCustomisableComponents(keys);
+ keys.addOutputComponent("bias","default","the instantaneous value of the bias potential");
+ keys.addOutputComponent("sigma", "default","uncertainty parameter");
+ keys.addOutputComponent("scale", "default","scale parameter");
+ keys.addOutputComponent("accept","default","MC acceptance");
+}
+
+Metainference::Metainference(const ActionOptions&ao):
+PLUMED_BIAS_INIT(ao),
+sqrt2pi(2.506628274631001),
+sqrt2_div_pi(0.45015815807855),
+doscale_(false),
+ndata_(getNumberOfArguments()),
+old_energy(0),
+MCsteps_(1),
+MCstride_(1),
+MCaccept_(0),
+MCfirst_(-1)
+{
+ parseVector("PARAMETERS",parameters);
+ if(parameters.size()!=static_cast<unsigned>(getNumberOfArguments())&&!parameters.empty())
+ error("Size of PARAMETERS array should be either 0 or the same as of the number of arguments in ARG1");
+
+ vector<Value*> arg2;
+ parseArgumentList("PARARG",arg2);
+
+ if(!arg2.empty()) {
+ if(parameters.size()>0) error("It is not possible to use PARARG and PARAMETERS together");
+ if(arg2.size()!=getNumberOfArguments()) error("Size of PARARG array should be the same as number for arguments in ARG");
+ for(unsigned i=0;i<arg2.size();i++){
+ parameters.push_back(arg2[i]->get());
+ if(arg2[i]->hasDerivatives()==true) error("PARARG can only accept arguments without derivatives");
+ }
+ }
+
+ if(parameters.size()!=getNumberOfArguments())
+ error("PARARG or PARAMETERS arrays should include the same number of elements as the arguments in ARG");
+
+ string stringa_noise;
+ parse("NOISETYPE",stringa_noise);
+ if(stringa_noise=="GAUSS") noise_type_ = GAUSS;
+ else if(stringa_noise=="MGAUSS") noise_type_ = MGAUSS;
+ else if(stringa_noise=="LTAIL") noise_type_ = LTAIL;
+ else error("Unkwnow noise type");
+
+ parseFlag("SCALEDATA", doscale_);
+ if(doscale_) {
+ parse("SCALE0",scale_);
+ parse("SCALE_MIN",scale_min_);
+ parse("SCALE_MAX",scale_max_);
+ parse("DSCALE",Dscale_);
+ } else {
+ scale_=1.0;
+ }
+
+ vector<double> readsigma;
+ parseVector("SIGMA0",readsigma);
+ if(noise_type_!=MGAUSS&&readsigma.size()>1) error("If you want to use more than one sigma you should add MULTISIGMA");
+
+ parse("SIGMA_MIN",sigma_min_);
+ parse("SIGMA_MAX",sigma_max_);
+ parse("DSIGMA",Dsigma_);
+ parse("MC_STEPS",MCsteps_);
+ parse("MC_STRIDE",MCstride_);
+ // get temperature
+ double temp=0.0;
+ parse("TEMP",temp);
+ parse("SIGMA_MEAN",sigma_mean_);
+
+ checkRead();
+
+ if(temp>0.0) kbt_=plumed.getAtoms().getKBoltzmann()*temp;
+ else kbt_=plumed.getAtoms().getKbT();
+
+ if(readsigma.size()==getNumberOfArguments()) {
+ sigma_.resize(getNumberOfArguments());
+ sigma_=readsigma;
+ } else if(readsigma.size()==1) {
+ if(noise_type_==MGAUSS) sigma_.resize(getNumberOfArguments(),readsigma[0]);
+ else sigma_.resize(1, readsigma[0]);
+ } else {
+ error("SIGMA0 can accept either one single value or as many values as the number of arguments (with MULTISIGMA)");
+ }
+
+ // get number of replicas
+ if(comm.Get_rank()==0) {
+ nrep_ = multi_sim_comm.Get_size();
+ replica_ = multi_sim_comm.Get_rank();
+ } else {
+ nrep_ = 0;
+ replica_ = 0;
+ }
+ comm.Sum(&nrep_,1);
+ comm.Sum(&replica_,1);
+
+ // divide sigma_mean by the square root of the number of replicas
+ sigma_mean_ /= sqrt(static_cast<double>(nrep_));
+
+ // adjust for multiple-time steps
+ MCstride_ *= getStride();
+
+ switch(noise_type_) {
+ case GAUSS:
+ log.printf(" with gaussian noise and a single noise parameter for all the data\n");
+ break;
+ case MGAUSS:
+ log.printf(" with gaussian noise and a noise parameter for each data point\n");
+ break;
+ case LTAIL:
+ log.printf(" with long tailed gaussian noise and a single noise parameter for all the data\n");
+ break;
+ }
+
+ if(doscale_) {
+ log.printf(" sampling a common scaling factor with:\n");
+ log.printf(" initial scale parameter %f\n",scale_);
+ log.printf(" minimum scale parameter %f\n",scale_min_);
+ log.printf(" maximum scale parameter %f\n",scale_max_);
+ log.printf(" maximum MC move of scale parameter %f\n",Dscale_);
+ }
+
+ if(readsigma.size()==1) log.printf(" initial data uncertainty %f\n",sigma_[0]);
+ else {
+ log.printf(" initial data uncertainties");
+ for(unsigned i=0;i<sigma_.size();++i) log.printf(" %f", sigma_[i]);
+ log.printf("\n");
+ }
+ log.printf(" minimum data uncertainty %f\n",sigma_min_);
+ log.printf(" maximum data uncertainty %f\n",sigma_max_);
+ log.printf(" maximum MC move of data uncertainty %f\n",Dsigma_);
+ log.printf(" uncertainty in the mean estimate %f\n",sigma_mean_);
+ log.printf(" temperature of the system %f\n",kbt_);
+ log.printf(" number of experimental data points %u\n",getNumberOfArguments());
+ log.printf(" number of replicas %u\n",nrep_);
+ log.printf(" MC steps %u\n",MCsteps_);
+ log.printf(" MC stride %u\n",MCstride_);
+
+ addComponent("bias");
+ componentIsNotPeriodic("bias");
+ valueBias=getPntrToComponent("bias");
+ if(doscale_) {
+ addComponent("scale");
+ componentIsNotPeriodic("scale");
+ valueScale=getPntrToComponent("scale");
+ }
+ addComponent("accept");
+ componentIsNotPeriodic("accept");
+ valueAccept=getPntrToComponent("accept");
+
+ if(noise_type_==MGAUSS) {
+ for(unsigned i=0;i<sigma_.size();++i){
+ std::string num; Tools::convert(i,num);
+ addComponent("sigma_"+num); componentIsNotPeriodic("sigma_"+num);
+ valueSigma.push_back(getPntrToComponent("sigma_"+num));
+ }
+ } else {
+ addComponent("sigma"); componentIsNotPeriodic("sigma");
+ valueSigma.push_back(getPntrToComponent("sigma"));
+ }
+ // initialize random seed
+ unsigned iseed;
+ if(comm.Get_rank()==0) iseed = time(NULL)+replica_;
+ else iseed = 0;
+ comm.Sum(&iseed, 1);
+ srand(iseed);
+}
+
+double Metainference::getEnergySPE(const double sigma, const double scale){
+ // calculate effective sigma
+ const double smean2 = sigma_mean_*sigma_mean_;
+ const double s = sqrt( sigma*sigma + smean2 );
+ // cycle on arguments
+ double ene = 0.0;
+ for(unsigned i=0;i<getNumberOfArguments();++i){
+ const double dev = scale*getArgument(i)-parameters[i];
+ // argument
+ const double a2 = 0.5*dev*dev + s*s;
+ // increment energy
+ ene += std::log( 2.0 * a2 / ( 1.0 - exp(- a2 / smean2) ) );
+ }
+ // add normalization and Jeffrey's prior
+ ene += std::log(s) - static_cast<double>(ndata_)*std::log(sqrt2_div_pi*s);
+ return kbt_ * ene;
+}
+
+double Metainference::getEnergyGJE(const vector<double> &sigma, const double scale){
+ // cycle on arguments
+ double ene = 0.0;
+ const double smean2 = sigma_mean_*sigma_mean_;
+ double ss = sigma[0]*sigma[0] + smean2;
+ for(unsigned i=0;i<getNumberOfArguments();++i){
+ if(noise_type_==MGAUSS) ss = sigma[i]*sigma[i] + smean2;
+ const double dev = scale*getArgument(i)-parameters[i];
+ ene += 0.5*dev*dev/ss + std::log(ss*sqrt2pi);
+ }
+ return kbt_ * ene;
+}
+
+void Metainference::doMonteCarlo(){
+ // calculate old energy (first time only)
+ if(MCfirst_==-1) {
+ switch(noise_type_) {
+ case GAUSS:
+ case MGAUSS:
+ old_energy = getEnergyGJE(sigma_,scale_);
+ break;
+ case LTAIL:
+ old_energy = getEnergySPE(sigma_[0],scale_);
+ break;
+ }
+ }
+
+ // cycle on MC steps
+ for(unsigned i=0;i<MCsteps_;++i){
+
+ // propose move for scale
+ double new_scale = scale_;
+ if(doscale_) {
+ const double r1 = static_cast<double>(rand()) / RAND_MAX;
+ const double ds1 = -Dscale_ + r1 * 2.0 * Dscale_;
+ new_scale += ds1;
+ // check boundaries
+ if(new_scale > scale_max_){new_scale = 2.0 * scale_max_ - new_scale;}
+ if(new_scale < scale_min_){new_scale = 2.0 * scale_min_ - new_scale;}
+ // the scaling factor should be the same for all the replicas
+ if(comm.Get_rank()==0) multi_sim_comm.Bcast(new_scale,0);
+ comm.Bcast(new_scale,0);
+ }
+
+ // propose move for sigma
+ vector<double> new_sigma(sigma_.size());
+ for(unsigned j=0;j<sigma_.size();++j) {
+ const double r2 = static_cast<double>(rand()) / RAND_MAX;
+ const double ds2 = -Dsigma_ + r2 * 2.0 * Dsigma_;
+ new_sigma[j] = sigma_[j] + ds2;
+ // check boundaries
+ if(new_sigma[j] > sigma_max_){new_sigma[j] = 2.0 * sigma_max_ - new_sigma[j];}
+ if(new_sigma[j] < sigma_min_){new_sigma[j] = 2.0 * sigma_min_ - new_sigma[j];}
+ }
+
+ // calculate new energy
+ double new_energy=0;
+ switch(noise_type_) {
+ case GAUSS:
+ case MGAUSS:
+ new_energy = getEnergyGJE(new_sigma,new_scale);
+ break;
+ case LTAIL:
+ new_energy = getEnergySPE(new_sigma[0],new_scale);
+ break;
+ }
+
+ // accept or reject
+ const double delta = ( new_energy - old_energy ) / kbt_;
+ // if delta is negative always accept move
+ if( delta <= 0.0 ){
+ old_energy = new_energy;
+ scale_ = new_scale;
+ sigma_ = new_sigma;
+ MCaccept_++;
+ // otherwise extract random number
+ } else {
+ const double s = static_cast<double>(rand()) / RAND_MAX;
+ if( s < exp(-delta) ){
+ old_energy = new_energy;
+ scale_ = new_scale;
+ sigma_ = new_sigma;
+ MCaccept_++;
+ }
+ }
+
+ if(doscale_) {
+ // the scaling factor should be the same for all the replicas
+ if(comm.Get_rank()==0) multi_sim_comm.Bcast(scale_,0);
+ comm.Bcast(scale_,0);
+ }
+ }
+ /* save the result of the sampling */
+ if(doscale_) valueScale->set(scale_);
+ for(unsigned i=0; i<sigma_.size(); i++) valueSigma[i]->set(sigma_[i]);
+}
+
+double Metainference::getEnergyForceSPE()
+{
+ double ene = 0.0;
+ const unsigned narg=getNumberOfArguments();
+
+ const double smean2 = sigma_mean_*sigma_mean_;
+ const double s = sqrt( sigma_[0]*sigma_[0] + smean2 );
+ vector<double> f(narg,0);
+
+ if(comm.Get_rank()==0){
+ for(unsigned i=0;i<narg;++i){
+ const double dev = scale_*getArgument(i)-parameters[i];
+ const double a2 = 0.5*dev*dev + s*s;
+ const double t = exp(-a2/smean2);
+ const double dt = 1./t;
+ const double it = 1./(1.-t);
+ const double dit = 1./(1.-dt);
+ ene += std::log(2.*a2*it);
+ f[i] = -scale_*dev*(dit/smean2 + 1./a2);
+ }
+ // collect contribution to forces and energy from other replicas
+ multi_sim_comm.Sum(&f[0],narg);
+ multi_sim_comm.Sum(&ene,1);
+ // add normalizations and priors of local replica
+ ene += std::log(s) - static_cast<double>(ndata_)*std::log(sqrt2_div_pi*s);
+ }
+ // intra-replica summation
+ comm.Sum(&f[0],narg);
+ comm.Sum(&ene,1);
+
+ for(unsigned i=0; i<narg; ++i) setOutputForce(i, kbt_ * f[i]);
+ return ene;
+}
+
+double Metainference::getEnergyForceGJE()
+{
+ double ene = 0.0;
+
+ const unsigned ssize = sigma_.size();
+ vector<double> ss(ssize);
+ vector<double> inv_s2(ssize, 0.);
+ const double smean2 = sigma_mean_*sigma_mean_;
+
+ for(unsigned i=0;i<ssize; ++i) {
+ ss[i] = sigma_[i]*sigma_[i] + smean2;
+ if(comm.Get_rank()==0) inv_s2[i] = 1.0/ss[i];
+ }
+
+ if(comm.Get_rank()==0) multi_sim_comm.Sum(&inv_s2[0],ssize);
+ comm.Sum(&inv_s2[0],ssize);
+
+ const unsigned narg=getNumberOfArguments();
+ for(unsigned i=0;i<narg;++i){
+ const double dev = scale_*getArgument(i)-parameters[i];
+ unsigned sel_sigma=0;
+ if(noise_type_==MGAUSS) sel_sigma=i;
+ ene += 0.5*dev*dev*inv_s2[sel_sigma] + std::log(ss[sel_sigma]*sqrt2pi);
+ setOutputForce(i, -kbt_*dev*scale_*inv_s2[sel_sigma]);
+ }
+ return ene;
+}
+
+void Metainference::calculate(){
+ /* MONTE CARLO */
+ const long int step = getStep();
+ if(step%MCstride_==0&&!getExchangeStep()) doMonteCarlo();
+ // this is needed when restarting simulations
+ if(MCfirst_==-1) MCfirst_=step;
+ const double MCtrials = floor(static_cast<double>(step-MCfirst_) / static_cast<double>(MCstride_))+1.0;
+ const double accept = static_cast<double>(MCaccept_) / static_cast<double>(MCsteps_) / MCtrials;
+ valueAccept->set(accept);
+
+ // calculate bias and forces
+ double ene = 0;
+ switch(noise_type_) {
+ case GAUSS:
+ case MGAUSS:
+ ene = getEnergyForceGJE();
+ break;
+ case LTAIL:
+ ene = getEnergyForceSPE();
+ break;
+ }
+ // set value of the bias
+ valueBias->set(kbt_*ene);
+}
+
+}
+}
+
+
--
GitLab