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 : PLUMED_REGISTER_ACTION(TD_Gaussian,"TD_GAUSSIAN")
177 :
178 :
179 193 : void TD_Gaussian::registerKeywords(Keywords& keys) {
180 193 : TargetDistribution::registerKeywords(keys);
181 193 : keys.add("numbered","CENTER","The centers of the Gaussian distributions.");
182 193 : keys.add("numbered","SIGMA","The standard deviations of the Gaussian distributions.");
183 193 : 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 193 : 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 193 : keys.use("WELLTEMPERED_FACTOR");
186 193 : keys.use("SHIFT_TO_ZERO");
187 193 : keys.use("NORMALIZE");
188 193 : }
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 211 : for(unsigned int i=1;; i++) {
200 : std::vector<double> tmp_center;
201 804 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {
202 : break;
203 : }
204 211 : centers_.push_back(tmp_center);
205 211 : }
206 211 : for(unsigned int i=1;; i++) {
207 : std::vector<double> tmp_sigma;
208 804 : if(!parseNumberedVector("SIGMA",i,tmp_sigma) ) {
209 : break;
210 : }
211 211 : sigmas_.push_back(tmp_sigma);
212 211 : }
213 :
214 191 : if(centers_.size()==0) {
215 0 : plumed_merror(getName()+": CENTER keywords seem to be missing. Note that numbered keywords start at CENTER1.");
216 : }
217 : //
218 191 : if(centers_.size()!=sigmas_.size()) {
219 0 : plumed_merror(getName()+": there has to be an equal amount of CENTER and SIGMA keywords");
220 : }
221 : //
222 191 : setDimension(centers_[0].size());
223 191 : ncenters_ = centers_.size();
224 : // check centers and sigmas
225 402 : for(unsigned int i=0; i<ncenters_; i++) {
226 211 : if(centers_[i].size()!=getDimension()) {
227 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
228 : }
229 211 : if(sigmas_[i].size()!=getDimension()) {
230 0 : plumed_merror(getName()+": one of the SIGMA keyword does not match the given dimension");
231 : }
232 : }
233 : //
234 191 : correlation_.resize(ncenters_);
235 :
236 402 : for(unsigned int i=0; i<ncenters_; i++) {
237 : std::vector<double> corr;
238 422 : parseNumberedVector("CORRELATION",(i+1),corr);
239 211 : if(corr.size()>0) {
240 3 : diagonal_ = false;
241 : } else {
242 208 : corr.assign(1,0.0);
243 : }
244 211 : correlation_[i] = corr;
245 : }
246 :
247 191 : if(!diagonal_ && getDimension()!=2) {
248 0 : plumed_merror(getName()+": CORRELATION is only defined for two-dimensional Gaussians for now.");
249 : }
250 402 : for(unsigned int i=0; i<correlation_.size(); i++) {
251 211 : if(correlation_[i].size()!=1) {
252 0 : plumed_merror(getName()+": only one value should be given in CORRELATION");
253 : }
254 422 : for(unsigned int k=0; k<correlation_[i].size(); k++) {
255 211 : if(correlation_[i][k] <= -1.0 || correlation_[i][k] >= 1.0) {
256 0 : plumed_merror(getName()+": values given in CORRELATION should be between -1.0 and 1.0" );
257 : }
258 : }
259 : }
260 : //
261 382 : parseVector("WEIGHTS",weights_);
262 191 : if(weights_.size()==0) {
263 188 : weights_.assign(centers_.size(),1.0);
264 : }
265 191 : if(centers_.size()!=weights_.size()) {
266 0 : plumed_merror(getName()+": there has to be as many weights given in WEIGHTS as numbered CENTER keywords");
267 : }
268 : //
269 : double sum_weights=0.0;
270 402 : for(unsigned int i=0; i<weights_.size(); i++) {
271 211 : sum_weights+=weights_[i];
272 : }
273 402 : for(unsigned int i=0; i<weights_.size(); i++) {
274 211 : weights_[i]/=sum_weights;
275 : }
276 : //
277 191 : checkRead();
278 191 : }
279 :
280 :
281 229739 : double TD_Gaussian::getValue(const std::vector<double>& argument) const {
282 : double value=0.0;
283 229739 : if(diagonal_) {
284 501589 : for(unsigned int i=0; i<ncenters_; i++) {
285 292252 : value+=weights_[i]*GaussianDiagonal(argument, centers_[i], sigmas_[i]);
286 : }
287 20402 : } else if(!diagonal_ && getDimension()==2) {
288 61206 : for(unsigned int i=0; i<ncenters_; i++) {
289 40804 : value+=weights_[i]*Gaussian2D(argument, centers_[i], sigmas_[i],correlation_[i]);
290 : }
291 : }
292 229739 : return value;
293 : }
294 :
295 :
296 292252 : double TD_Gaussian::GaussianDiagonal(const std::vector<double>& argument, const std::vector<double>& center, const std::vector<double>& sigma, bool normalize) const {
297 : double value = 1.0;
298 828323 : for(unsigned int k=0; k<argument.size(); k++) {
299 536071 : double arg=(argument[k]-center[k])/sigma[k];
300 536071 : double tmp_exp = exp(-0.5*arg*arg);
301 536071 : if(normalize) {
302 536071 : tmp_exp/=(sigma[k]*sqrt(2.0*pi));
303 : }
304 536071 : value*=tmp_exp;
305 : }
306 292252 : return value;
307 : }
308 :
309 :
310 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 {
311 40804 : double arg1 = (argument[0]-center[0])/sigma[0];
312 40804 : double arg2 = (argument[1]-center[1])/sigma[1];
313 40804 : double corr = correlation[0];
314 40804 : double value = (arg1*arg1 + arg2*arg2 - 2.0*corr*arg1*arg2);
315 40804 : value *= -1.0 / ( 2.0*(1.0-corr*corr) );
316 40804 : value = exp(value);
317 40804 : if(normalize) {
318 40804 : value /= 2*pi*sigma[0]*sigma[1]*sqrt(1.0-corr*corr);
319 : }
320 40804 : return value;
321 : }
322 :
323 : }
324 : }
|
