Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 The VES code team
3 : (see the PEOPLE-VES file at the root of this folder for a list of names)
4 :
5 : See http://www.ves-code.org for more information.
6 :
7 : This file is part of VES code module.
8 :
9 : The VES code module 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 : The VES code module 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 the VES code module. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 :
23 : #include "TargetDistribution.h"
24 :
25 : #include "core/ActionRegister.h"
26 :
27 :
28 : namespace PLMD {
29 : namespace ves {
30 :
31 : //+PLUMEDOC VES_TARGETDIST TD_GAUSSIAN
32 : /*
33 : Target distribution given by a sum of Gaussians (static).
34 :
35 : Employ a target distribution that is given by a sum of multivariate Gaussian (or normal)
36 : distributions, defined as
37 : \f[
38 : p(\mathbf{s}) = \sum_{i} \, w_{i} \, N(\mathbf{s};\mathbf{\mu}_{i},\mathbf{\Sigma}_{i})
39 : \f]
40 : where \f$\mathbf{\mu}_{i}=(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})\f$
41 : and \f$\mathbf{\Sigma}_{i}\f$ are
42 : the center and the covariance matrix for the \f$i\f$-th Gaussian.
43 : The weights \f$w_{i}\f$ are normalized to 1, \f$\sum_{i}w_{i}=1\f$.
44 :
45 : By default the Gaussian distributions are considered as separable into
46 : independent one-dimensional Gaussian distributions. In other words,
47 : the covariance matrix is taken as diagonal
48 : \f$\mathbf{\Sigma}_{i}=(\sigma^2_{1,i},\sigma^2_{2,i},\ldots,\sigma^{2}_{d,i})\f$.
49 : The Gaussian distribution is then written as
50 : \f[
51 : N(\mathbf{s};\mathbf{\mu}_{i},\mathbf{\sigma}_{i}) =
52 : \prod^{d}_{k} \, \frac{1}{\sqrt{2\pi\sigma^2_{d,i}}} \,
53 : \exp\left(
54 : -\frac{(s_{d}-\mu_{d,i})^2}{2\sigma^2_{d,i}}
55 : \right)
56 : \f]
57 : where
58 : \f$\mathbf{\sigma}_{i}=(\sigma_{1,i},\sigma_{2,i},\ldots,\sigma_{d,i})\f$
59 : is the standard deviation.
60 : In this case you need to specify the centers \f$\mathbf{\mu}_{i}\f$ using the
61 : numbered CENTER keywords and the standard deviations \f$\mathbf{\sigma}_{i}\f$
62 : using the numbered SIGMA keywords.
63 :
64 : For two arguments it is possible to employ
65 : [bivariate Gaussians](https://en.wikipedia.org/wiki/Multivariate_normal_distribution)
66 : with correlation between arguments, defined as
67 : \f[
68 : N(\mathbf{s};\mathbf{\mu}_{i},\mathbf{\sigma}_{i},\rho_i) =
69 : \frac{1}{2 \pi \sigma_{1,i} \sigma_{2,i} \sqrt{1-\rho_i^2}}
70 : \,
71 : \exp\left(
72 : -\frac{1}{2(1-\rho_i^2)}
73 : \left[
74 : \frac{(s_{1}-\mu_{1,i})^2}{\sigma_{1,i}^2}+
75 : \frac{(s_{2}-\mu_{2,i})^2}{\sigma_{2,i}^2}-
76 : \frac{2 \rho_i (s_{1}-\mu_{1,i})(s_{2}-\mu_{2,i})}{\sigma_{1,i}\sigma_{2,i}}
77 : \right]
78 : \right)
79 : \f]
80 : where \f$\rho_i\f$ is the correlation between \f$s_{1}\f$ and \f$s_{2}\f$
81 : that goes from -1 to 1. In this case the covariance matrix is given as
82 : \f[
83 : \mathbf{\Sigma}=
84 : \left[
85 : \begin{array}{cc}
86 : \sigma^2_{1,i} & \rho_i \sigma_{1,i} \sigma_{2,i} \\
87 : \rho_i \sigma_{1,i} \sigma_{2,i} & \sigma^2_{2,i}
88 : \end{array}
89 : \right]
90 : \f]
91 : The correlation \f$\rho\f$ is given using
92 : the numbered CORRELATION keywords. A value of \f$\rho=0\f$ means
93 : that the arguments are considered as
94 : un-correlated, which is the default behavior.
95 :
96 : The Gaussian distributions are always defined with the conventional
97 : normalization factor such that they are normalized to 1 over an unbounded
98 : region. However, in calculation within VES we normally consider bounded
99 : region on which the target distribution is defined. Thus, if the center of
100 : a Gaussian is close to the boundary of the region it can happen that the
101 : tails go outside the region. In that case it might be needed to use the
102 : NORMALIZE keyword to make sure that the target distribution is properly
103 : normalized to 1 over the bounded region. The code will issue a warning
104 : if that is needed.
105 :
106 : For periodic CVs it is generally better to use \ref TD_VONMISES "Von Mises"
107 : distributions instead of Gaussians as these distributions properly
108 : account for the periodicity of the CVs.
109 :
110 :
111 : \par Examples
112 :
113 : One single Gaussians in one-dimension.
114 : \plumedfile
115 : td: TD_GAUSSIAN CENTER1=-1.5 SIGMA1=0.8
116 : \endplumedfile
117 :
118 : Sum of three Gaussians in two-dimensions with equal weights as
119 : no weights are given.
120 : \plumedfile
121 : TD_GAUSSIAN ...
122 : CENTER1=-1.5,+1.5 SIGMA1=0.8,0.3
123 : CENTER2=+1.5,-1.5 SIGMA2=0.3,0.8
124 : CENTER3=+1.5,+1.5 SIGMA3=0.4,0.4
125 : LABEL=td
126 : ... TD_GAUSSIAN
127 : \endplumedfile
128 :
129 : Sum of three Gaussians in two-dimensions which
130 : are weighted unequally. Note that weights are automatically
131 : normalized to 1 so that WEIGHTS=1.0,2.0,1.0 is equal to
132 : specifying WEIGHTS=0.25,0.50,0.25.
133 : \plumedfile
134 : TD_GAUSSIAN ...
135 : CENTER1=-1.5,+1.5 SIGMA1=0.8,0.3
136 : CENTER2=+1.5,-1.5 SIGMA2=0.3,0.8
137 : CENTER3=+1.5,+1.5 SIGMA3=0.4,0.4
138 : WEIGHTS=1.0,2.0,1.0
139 : LABEL=td
140 : ... TD_GAUSSIAN
141 : \endplumedfile
142 :
143 : Sum of two bivariate Gaussians where there is correlation of
144 : \f$\rho_{2}=0.75\f$ between the two arguments for the second Gaussian.
145 : \plumedfile
146 : TD_GAUSSIAN ...
147 : CENTER1=-1.5,+1.5 SIGMA1=0.8,0.3
148 : CENTER2=+1.5,-1.5 SIGMA2=0.3,0.8 CORRELATION2=0.75
149 : LABEL=td
150 : ... TD_GAUSSIAN
151 : \endplumedfile
152 :
153 :
154 :
155 :
156 :
157 : */
158 : //+ENDPLUMEDOC
159 :
160 429 : class TD_Gaussian: public TargetDistribution {
161 : std::vector< std::vector<double> > centers_;
162 : std::vector< std::vector<double> > sigmas_;
163 : std::vector< std::vector<double> > correlation_;
164 : std::vector<double> weights_;
165 : bool diagonal_;
166 : unsigned int ncenters_;
167 : double GaussianDiagonal(const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const bool normalize=true) const;
168 : double Gaussian2D(const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const bool normalize=true) const;
169 : public:
170 : static void registerKeywords(Keywords&);
171 : explicit TD_Gaussian(const ActionOptions& ao);
172 : double getValue(const std::vector<double>&) const;
173 : };
174 :
175 :
176 6595 : PLUMED_REGISTER_ACTION(TD_Gaussian,"TD_GAUSSIAN")
177 :
178 :
179 144 : void TD_Gaussian::registerKeywords(Keywords& keys) {
180 144 : TargetDistribution::registerKeywords(keys);
181 576 : keys.add("numbered","CENTER","The centers of the Gaussian distributions.");
182 576 : keys.add("numbered","SIGMA","The standard deviations of the Gaussian distributions.");
183 576 : keys.add("numbered","CORRELATION","The correlation for two-dimensional bivariate Gaussian distributions. Only works for two arguments. The value should be between -1 and 1. If no value is given the Gaussians is considered as un-correlated (i.e. value of 0.0).");
184 576 : keys.add("optional","WEIGHTS","The weights of the Gaussian distributions. Have to be as many as the number of centers given with the numbered CENTER keywords. If no weights are given the distributions are weighted equally. The weights are automatically normalized to 1.");
185 288 : keys.use("WELLTEMPERED_FACTOR");
186 288 : keys.use("SHIFT_TO_ZERO");
187 288 : keys.use("NORMALIZE");
188 144 : }
189 :
190 :
191 143 : TD_Gaussian::TD_Gaussian(const ActionOptions& ao):
192 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
193 : centers_(0),
194 : sigmas_(0),
195 : correlation_(0),
196 : weights_(0),
197 : diagonal_(true),
198 143 : ncenters_(0)
199 : {
200 163 : for(unsigned int i=1;; i++) {
201 : std::vector<double> tmp_center;
202 612 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {break;}
203 163 : centers_.push_back(tmp_center);
204 163 : }
205 163 : for(unsigned int i=1;; i++) {
206 : std::vector<double> tmp_sigma;
207 612 : if(!parseNumberedVector("SIGMA",i,tmp_sigma) ) {break;}
208 163 : sigmas_.push_back(tmp_sigma);
209 163 : }
210 :
211 143 : if(centers_.size()==0) {
212 0 : plumed_merror(getName()+": CENTER keywords seem to be missing. Note that numbered keywords start at CENTER1.");
213 : }
214 : //
215 143 : if(centers_.size()!=sigmas_.size()) {
216 0 : plumed_merror(getName()+": there has to be an equal amount of CENTER and SIGMA keywords");
217 : }
218 : //
219 143 : setDimension(centers_[0].size());
220 143 : ncenters_ = centers_.size();
221 : // check centers and sigmas
222 469 : for(unsigned int i=0; i<ncenters_; i++) {
223 326 : if(centers_[i].size()!=getDimension()) {
224 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
225 : }
226 163 : if(sigmas_[i].size()!=getDimension()) {
227 0 : plumed_merror(getName()+": one of the SIGMA keyword does not match the given dimension");
228 : }
229 : }
230 : //
231 143 : correlation_.resize(ncenters_);
232 :
233 469 : for(unsigned int i=0; i<ncenters_; i++) {
234 : std::vector<double> corr;
235 326 : parseNumberedVector("CORRELATION",(i+1),corr);
236 163 : if(corr.size()>0) {
237 3 : diagonal_ = false;
238 : }
239 : else {
240 320 : corr.assign(1,0.0);
241 : }
242 326 : correlation_[i] = corr;
243 : }
244 :
245 145 : if(!diagonal_ && getDimension()!=2) {
246 0 : plumed_merror(getName()+": CORRELATION is only defined for two-dimensional Gaussians for now.");
247 : }
248 775 : for(unsigned int i=0; i<correlation_.size(); i++) {
249 163 : if(correlation_[i].size()!=1) {
250 0 : plumed_merror(getName()+": only one value should be given in CORRELATION");
251 : }
252 489 : for(unsigned int k=0; k<correlation_[i].size(); k++) {
253 163 : if(correlation_[i][k] <= -1.0 || correlation_[i][k] >= 1.0) {
254 0 : plumed_merror(getName()+": values given in CORRELATION should be between -1.0 and 1.0" );
255 : }
256 : }
257 : }
258 : //
259 286 : parseVector("WEIGHTS",weights_);
260 283 : if(weights_.size()==0) {weights_.assign(centers_.size(),1.0);}
261 143 : if(centers_.size()!=weights_.size()) {
262 0 : plumed_merror(getName()+": there has to be as many weights given in WEIGHTS as numbered CENTER keywords");
263 : }
264 : //
265 : double sum_weights=0.0;
266 469 : for(unsigned int i=0; i<weights_.size(); i++) {sum_weights+=weights_[i];}
267 775 : for(unsigned int i=0; i<weights_.size(); i++) {weights_[i]/=sum_weights;}
268 : //
269 143 : checkRead();
270 143 : }
271 :
272 :
273 215307 : double TD_Gaussian::getValue(const std::vector<double>& argument) const {
274 : double value=0.0;
275 215307 : if(diagonal_) {
276 750545 : for(unsigned int i=0; i<ncenters_; i++) {
277 833460 : value+=weights_[i]*GaussianDiagonal(argument, centers_[i], sigmas_[i]);
278 : }
279 : }
280 20402 : else if(!diagonal_ && getDimension()==2) {
281 102010 : for(unsigned int i=0; i<ncenters_; i++) {
282 122412 : value+=weights_[i]*Gaussian2D(argument, centers_[i], sigmas_[i],correlation_[i]);
283 : }
284 : }
285 215307 : return value;
286 : }
287 :
288 :
289 277820 : double TD_Gaussian::GaussianDiagonal(const std::vector<double>& argument, const std::vector<double>& center, const std::vector<double>& sigma, bool normalize) const {
290 : double value = 1.0;
291 2120557 : for(unsigned int k=0; k<argument.size(); k++) {
292 1564917 : double arg=(argument[k]-center[k])/sigma[k];
293 521639 : double tmp_exp = exp(-0.5*arg*arg);
294 1043278 : if(normalize) {tmp_exp/=(sigma[k]*sqrt(2.0*pi));}
295 521639 : value*=tmp_exp;
296 : }
297 277820 : return value;
298 : }
299 :
300 :
301 40804 : double TD_Gaussian::Gaussian2D(const std::vector<double>& argument, const std::vector<double>& center, const std::vector<double>& sigma, const std::vector<double>& correlation, bool normalize) const {
302 122412 : double arg1 = (argument[0]-center[0])/sigma[0];
303 40804 : double arg2 = (argument[1]-center[1])/sigma[1];
304 40804 : double corr = correlation[0];
305 40804 : double value = (arg1*arg1 + arg2*arg2 - 2.0*corr*arg1*arg2);
306 40804 : value *= -1.0 / ( 2.0*(1.0-corr*corr) );
307 40804 : value = exp(value);
308 40804 : if(normalize) {
309 40804 : value /= 2*pi*sigma[0]*sigma[1]*sqrt(1.0-corr*corr);
310 : }
311 40804 : return value;
312 : }
313 :
314 : }
315 4839 : }
|
