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