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 : #include "GridIntegrationWeights.h"
25 :
26 : #include "core/ActionRegister.h"
27 : #include "tools/Tools.h"
28 :
29 : #include <iostream>
30 :
31 :
32 :
33 : namespace PLMD {
34 : namespace ves {
35 :
36 : //+PLUMEDOC VES_TARGETDIST TD_VONMISES
37 : /*
38 : Target distribution given by a sum of Von Mises distributions (static).
39 :
40 : Employ a target distribution that is given by a sum where each
41 : term is a product of one-dimensional
42 : [Von Mises distributions](https://en.wikipedia.org/wiki/Von_Mises_distribution),
43 :
44 : $$
45 : p(\mathbf{s}) = \sum_{i} \, w_{i}
46 : \prod_{k}^{d}
47 : \frac{\exp\left(\kappa_{k,i} \, \cos (s_{k}-\mu_{k,i}) \right)}
48 : {2\pi I_{0}(\kappa_{k,i})}
49 : $$
50 :
51 : where $(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})$
52 : are the centers of the distributions,
53 : $(\kappa_{1,i},\kappa_{2,i},\ldots,\kappa_{d,i})$
54 : are parameters that determine the extend of each distribution,
55 : and $I_{0}(x)$ is the modified Bessel function of order 0.
56 : The weights $w_{i}$ are normalized to 1, $\sum_{i}w_{i}=1$.
57 :
58 : The Von Mises distribution is defined for periodic variables with a
59 : periodicity of $2\pi$ and is analogous to the Gaussian distribution.
60 : The parameter $ \sqrt{1/\kappa}$ is comparable to the standard deviation
61 : $\sigma$ for the Gaussian distribution.
62 :
63 : To use this target distribution you need to give the centers
64 : $(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})$ by
65 : using the numbered CENTER keywords and the "standard deviations"
66 : $(\sqrt{1/\kappa_{1,i}},\sqrt{1/\kappa_{2,i}},\ldots,\sqrt{1/\kappa_{d,i}})$ using the numbered SIGMA keywords.
67 :
68 : ## Examples
69 :
70 : Sum of two Von Mises distribution in one dimension that have equal weights
71 : as no weights are given.
72 :
73 : ```plumed
74 : TD_VONMISES ...
75 : CENTER1=+2.0 SIGMA1=0.6
76 : CENTER2=-2.0 SIGMA2=0.7
77 : LABEL=td
78 : ... TD_VONMISES
79 : ```
80 :
81 : Sum of two Von Mises distribution in two dimensions that have different weights.
82 : Note that the weights are automatically normalized to 1 such that
83 : specifying WEIGHTS=1.0,2.0 is equal to specifying WEIGHTS=0.33333,0.66667.
84 :
85 : ```plumed
86 : TD_VONMISES ...
87 : CENTER1=+2.0,+2.0 SIGMA1=0.6,0.7
88 : CENTER2=-2.0,+2.0 SIGMA2=0.7,0.6
89 : WEIGHTS=1.0,2.0
90 : LABEL=td
91 : ... TD_VONMISES
92 : ```
93 :
94 : */
95 : //+ENDPLUMEDOC
96 :
97 : class TD_VonMises: public TargetDistribution {
98 : // properties of the Gaussians
99 : std::vector< std::vector<double> > sigmas_;
100 : std::vector< std::vector<double> > kappas_;
101 : std::vector< std::vector<double> > centers_;
102 : std::vector< std::vector<double> > normalization_;
103 : std::vector<double> weights_;
104 : std::vector<double> periods_;
105 : unsigned int ncenters_;
106 : double VonMisesDiagonal(const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&) const;
107 : double getNormalization(const double, const double) const;
108 : public:
109 : static void registerKeywords(Keywords&);
110 : explicit TD_VonMises(const ActionOptions& ao);
111 : double getValue(const std::vector<double>&) const override;
112 : };
113 :
114 :
115 : PLUMED_REGISTER_ACTION(TD_VonMises,"TD_VONMISES")
116 :
117 :
118 11 : void TD_VonMises::registerKeywords(Keywords& keys) {
119 11 : TargetDistribution::registerKeywords(keys);
120 11 : keys.add("numbered","CENTER","The centers of the Von Mises distributions.");
121 11 : keys.add("numbered","SIGMA","The standard deviations of the Von Mises distributions.");
122 11 : keys.add("optional","WEIGHTS","The weights of the Von Mises 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.");
123 11 : keys.add("hidden","PERIODS","The periods for each of the dimensions. By default they are 2*pi for each dimension.");
124 11 : keys.use("WELLTEMPERED_FACTOR");
125 11 : keys.use("SHIFT_TO_ZERO");
126 : //keys.use("NORMALIZE");
127 11 : }
128 :
129 :
130 9 : TD_VonMises::TD_VonMises(const ActionOptions& ao):
131 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
132 18 : sigmas_(0),
133 9 : centers_(0),
134 9 : normalization_(0),
135 9 : weights_(0),
136 9 : periods_(0),
137 18 : ncenters_(0) {
138 13 : for(unsigned int i=1;; i++) {
139 : std::vector<double> tmp_center;
140 44 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {
141 : break;
142 : }
143 13 : centers_.push_back(tmp_center);
144 13 : }
145 13 : for(unsigned int i=1;; i++) {
146 : std::vector<double> tmp_sigma;
147 44 : if(!parseNumberedVector("SIGMA",i,tmp_sigma) ) {
148 : break;
149 : }
150 13 : sigmas_.push_back(tmp_sigma);
151 13 : }
152 : //
153 9 : plumed_massert(centers_.size()==sigmas_.size(),"there has to be an equal amount of CENTER and SIGMA keywords");
154 9 : if(centers_.size()==0) {
155 0 : plumed_merror(getName()+": CENTER and SIGMA keywords seem to be missing. Note that numbered keywords start at CENTER1 and SIGMA1.");
156 : }
157 : //
158 9 : setDimension(centers_[0].size());
159 9 : ncenters_ = centers_.size();
160 : //
161 : // check centers and sigmas
162 22 : for(unsigned int i=0; i<ncenters_; i++) {
163 13 : if(centers_[i].size()!=getDimension()) {
164 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
165 : }
166 13 : if(sigmas_[i].size()!=getDimension()) {
167 0 : plumed_merror(getName()+": one of the SIGMA keyword does not match the given dimension");
168 : }
169 : }
170 : //
171 9 : kappas_.resize(sigmas_.size());
172 22 : for(unsigned int i=0; i<sigmas_.size(); i++) {
173 13 : kappas_[i].resize(sigmas_[i].size());
174 32 : for(unsigned int k=0; k<kappas_[i].size(); k++) {
175 19 : kappas_[i][k] = 1.0/(sigmas_[i][k]*sigmas_[i][k]);
176 : }
177 : }
178 : //
179 18 : parseVector("WEIGHTS",weights_);
180 9 : if(weights_.size()==0) {
181 9 : weights_.assign(centers_.size(),1.0);
182 : }
183 9 : if(centers_.size()!=weights_.size()) {
184 0 : plumed_merror(getName() + ": there has to be as many weights given in WEIGHTS as numbered CENTER keywords");
185 : }
186 : //
187 9 : if(periods_.size()==0) {
188 9 : periods_.assign(getDimension(),2*pi);
189 : }
190 18 : parseVector("PERIODS",periods_);
191 9 : if(periods_.size()!=getDimension()) {
192 0 : plumed_merror(getName() + ": the number of values given in PERIODS does not match the dimension of the distribution");
193 : }
194 : //
195 : double sum_weights=0.0;
196 22 : for(unsigned int i=0; i<weights_.size(); i++) {
197 13 : sum_weights+=weights_[i];
198 : }
199 22 : for(unsigned int i=0; i<weights_.size(); i++) {
200 13 : weights_[i]/=sum_weights;
201 : }
202 : //
203 9 : normalization_.resize(ncenters_);
204 22 : for(unsigned int i=0; i<ncenters_; i++) {
205 13 : normalization_[i].resize(getDimension());
206 32 : for(unsigned int k=0; k<getDimension(); k++) {
207 19 : normalization_[i][k] = getNormalization(kappas_[i][k],periods_[k]);
208 : }
209 : }
210 9 : checkRead();
211 9 : }
212 :
213 :
214 31100 : double TD_VonMises::getValue(const std::vector<double>& argument) const {
215 : double value=0.0;
216 92400 : for(unsigned int i=0; i<ncenters_; i++) {
217 61300 : value+=weights_[i]*VonMisesDiagonal(argument, centers_[i], kappas_[i],periods_,normalization_[i]);
218 : }
219 31100 : return value;
220 : }
221 :
222 :
223 80319 : double TD_VonMises::VonMisesDiagonal(const std::vector<double>& argument, const std::vector<double>& center, const std::vector<double>& kappa, const std::vector<double>& periods, const std::vector<double>& normalization) const {
224 : double value = 1.0;
225 220638 : for(unsigned int k=0; k<argument.size(); k++) {
226 140319 : double arg = kappa[k]*cos( ((2*pi)/periods[k])*(argument[k]-center[k]) );
227 140319 : value*=normalization[k]*exp(arg);
228 : }
229 80319 : return value;
230 : }
231 :
232 :
233 19 : double TD_VonMises::getNormalization(const double kappa, const double period) const {
234 : //
235 19 : std::vector<double> centers(1);
236 19 : centers[0] = 0.0;
237 19 : std::vector<double> kappas(1);
238 19 : kappas[0] = kappa;
239 19 : std::vector<double> periods(1);
240 19 : periods[0] = period;
241 19 : std::vector<double> norm(1);
242 19 : norm[0] = 1.0;
243 : //
244 : const unsigned int nbins = 1001;
245 : std::vector<double> points;
246 : std::vector<double> weights;
247 : double min = 0.0;
248 : double max = period;
249 19 : GridIntegrationWeights::getOneDimensionalIntegrationPointsAndWeights(points,weights,nbins,min,max);
250 : //
251 : double sum = 0.0;
252 19038 : for(unsigned int l=0; l<nbins; l++) {
253 19019 : std::vector<double> arg(1);
254 19019 : arg[0]= points[l];
255 19019 : sum += weights[l] * VonMisesDiagonal(arg,centers,kappas,periods,norm);
256 : }
257 38 : return 1.0/sum;
258 : }
259 :
260 :
261 : }
262 : }
|
