Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2021 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 Gaussian kernels (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 Gaussian kernels](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 Gaussian kernels as these distributions properly
108 : account for the periodicity of the CVs.
109 :
110 :
111 : \par Examples
112 :
113 : One single Gaussian kernel in one-dimension.
114 : \plumedfile
115 : td: TD_GAUSSIAN CENTER1=-1.5 SIGMA1=0.8
116 : \endplumedfile
117 :
118 : Sum of three Gaussian kernels 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 Gaussian kernels 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 Gaussian kernels 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 : 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 override;
173 : };
174 :
175 :
176 10610 : PLUMED_REGISTER_ACTION(TD_Gaussian,"TD_GAUSSIAN")
177 :
178 :
179 192 : void TD_Gaussian::registerKeywords(Keywords& keys) {
180 192 : TargetDistribution::registerKeywords(keys);
181 384 : keys.add("numbered","CENTER","The centers of the Gaussian distributions.");
182 384 : keys.add("numbered","SIGMA","The standard deviations of the Gaussian distributions.");
183 384 : 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 Gaussian kernels is considered as un-correlated (i.e. value of 0.0).");
184 384 : 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 192 : keys.use("WELLTEMPERED_FACTOR");
186 192 : keys.use("SHIFT_TO_ZERO");
187 192 : keys.use("NORMALIZE");
188 192 : }
189 :
190 :
191 191 : TD_Gaussian::TD_Gaussian(const ActionOptions& ao):
192 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
193 382 : centers_(0),
194 191 : sigmas_(0),
195 191 : correlation_(0),
196 191 : weights_(0),
197 191 : diagonal_(true),
198 382 : ncenters_(0)
199 : {
200 211 : for(unsigned int i=1;; i++) {
201 : std::vector<double> tmp_center;
202 804 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {break;}
203 211 : centers_.push_back(tmp_center);
204 211 : }
205 211 : for(unsigned int i=1;; i++) {
206 : std::vector<double> tmp_sigma;
207 804 : if(!parseNumberedVector("SIGMA",i,tmp_sigma) ) {break;}
208 211 : sigmas_.push_back(tmp_sigma);
209 211 : }
210 :
211 191 : 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 191 : 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 191 : setDimension(centers_[0].size());
220 191 : ncenters_ = centers_.size();
221 : // check centers and sigmas
222 402 : for(unsigned int i=0; i<ncenters_; i++) {
223 211 : if(centers_[i].size()!=getDimension()) {
224 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
225 : }
226 211 : 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 191 : correlation_.resize(ncenters_);
232 :
233 402 : for(unsigned int i=0; i<ncenters_; i++) {
234 : std::vector<double> corr;
235 422 : parseNumberedVector("CORRELATION",(i+1),corr);
236 211 : if(corr.size()>0) {
237 3 : diagonal_ = false;
238 : }
239 : else {
240 208 : corr.assign(1,0.0);
241 : }
242 211 : correlation_[i] = corr;
243 : }
244 :
245 191 : if(!diagonal_ && getDimension()!=2) {
246 0 : plumed_merror(getName()+": CORRELATION is only defined for two-dimensional Gaussians for now.");
247 : }
248 402 : for(unsigned int i=0; i<correlation_.size(); i++) {
249 211 : if(correlation_[i].size()!=1) {
250 0 : plumed_merror(getName()+": only one value should be given in CORRELATION");
251 : }
252 422 : for(unsigned int k=0; k<correlation_[i].size(); k++) {
253 211 : 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 382 : parseVector("WEIGHTS",weights_);
260 191 : if(weights_.size()==0) {weights_.assign(centers_.size(),1.0);}
261 191 : 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 402 : for(unsigned int i=0; i<weights_.size(); i++) {sum_weights+=weights_[i];}
267 402 : for(unsigned int i=0; i<weights_.size(); i++) {weights_[i]/=sum_weights;}
268 : //
269 191 : checkRead();
270 191 : }
271 :
272 :
273 229739 : double TD_Gaussian::getValue(const std::vector<double>& argument) const {
274 : double value=0.0;
275 229739 : if(diagonal_) {
276 501589 : for(unsigned int i=0; i<ncenters_; i++) {
277 292252 : value+=weights_[i]*GaussianDiagonal(argument, centers_[i], sigmas_[i]);
278 : }
279 : }
280 20402 : else if(!diagonal_ && getDimension()==2) {
281 61206 : for(unsigned int i=0; i<ncenters_; i++) {
282 40804 : value+=weights_[i]*Gaussian2D(argument, centers_[i], sigmas_[i],correlation_[i]);
283 : }
284 : }
285 229739 : return value;
286 : }
287 :
288 :
289 292252 : 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 828323 : for(unsigned int k=0; k<argument.size(); k++) {
292 536071 : double arg=(argument[k]-center[k])/sigma[k];
293 536071 : double tmp_exp = exp(-0.5*arg*arg);
294 536071 : if(normalize) {tmp_exp/=(sigma[k]*sqrt(2.0*pi));}
295 536071 : value*=tmp_exp;
296 : }
297 292252 : 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 40804 : 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 : }
|
