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 13 : void TD_VonMises::registerKeywords(Keywords& keys) {
119 13 : TargetDistribution::registerKeywords(keys);
120 13 : keys.add("numbered","CENTER","The centers of the Von Mises distributions.");
121 13 : keys.add("numbered","SIGMA","The standard deviations of the Von Mises distributions.");
122 13 : 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 13 : keys.add("hidden","PERIODS","The periods for each of the dimensions. By default they are 2*pi for each dimension.");
124 13 : keys.use("WELLTEMPERED_FACTOR");
125 13 : keys.use("SHIFT_TO_ZERO");
126 13 : }
127 :
128 :
129 9 : TD_VonMises::TD_VonMises(const ActionOptions& ao):
130 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
131 9 : sigmas_(0),
132 9 : centers_(0),
133 9 : normalization_(0),
134 9 : weights_(0),
135 9 : periods_(0),
136 18 : ncenters_(0) {
137 13 : for(unsigned int i=1;; i++) {
138 22 : std::vector<double> tmp_center;
139 44 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {
140 : break;
141 : }
142 13 : centers_.push_back(tmp_center);
143 13 : }
144 13 : for(unsigned int i=1;; i++) {
145 22 : std::vector<double> tmp_sigma;
146 44 : if(!parseNumberedVector("SIGMA",i,tmp_sigma) ) {
147 : break;
148 : }
149 13 : sigmas_.push_back(tmp_sigma);
150 13 : }
151 : //
152 9 : plumed_massert(centers_.size()==sigmas_.size(),"there has to be an equal amount of CENTER and SIGMA keywords");
153 9 : if(centers_.size()==0) {
154 0 : plumed_merror(getName()+": CENTER and SIGMA keywords seem to be missing. Note that numbered keywords start at CENTER1 and SIGMA1.");
155 : }
156 : //
157 9 : setDimension(centers_[0].size());
158 9 : ncenters_ = centers_.size();
159 : //
160 : // check centers and sigmas
161 22 : for(unsigned int i=0; i<ncenters_; i++) {
162 13 : if(centers_[i].size()!=getDimension()) {
163 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
164 : }
165 13 : if(sigmas_[i].size()!=getDimension()) {
166 0 : plumed_merror(getName()+": one of the SIGMA keyword does not match the given dimension");
167 : }
168 : }
169 : //
170 9 : kappas_.resize(sigmas_.size());
171 22 : for(unsigned int i=0; i<sigmas_.size(); i++) {
172 13 : kappas_[i].resize(sigmas_[i].size());
173 32 : for(unsigned int k=0; k<kappas_[i].size(); k++) {
174 19 : kappas_[i][k] = 1.0/(sigmas_[i][k]*sigmas_[i][k]);
175 : }
176 : }
177 : //
178 18 : parseVector("WEIGHTS",weights_);
179 9 : if(weights_.size()==0) {
180 9 : weights_.assign(centers_.size(),1.0);
181 : }
182 9 : if(centers_.size()!=weights_.size()) {
183 0 : plumed_merror(getName() + ": there has to be as many weights given in WEIGHTS as numbered CENTER keywords");
184 : }
185 : //
186 9 : if(periods_.size()==0) {
187 9 : periods_.assign(getDimension(),2*pi);
188 : }
189 18 : parseVector("PERIODS",periods_);
190 9 : if(periods_.size()!=getDimension()) {
191 0 : plumed_merror(getName() + ": the number of values given in PERIODS does not match the dimension of the distribution");
192 : }
193 : //
194 : double sum_weights=0.0;
195 22 : for(unsigned int i=0; i<weights_.size(); i++) {
196 13 : sum_weights+=weights_[i];
197 : }
198 22 : for(unsigned int i=0; i<weights_.size(); i++) {
199 13 : weights_[i]/=sum_weights;
200 : }
201 : //
202 9 : normalization_.resize(ncenters_);
203 22 : for(unsigned int i=0; i<ncenters_; i++) {
204 13 : normalization_[i].resize(getDimension());
205 32 : for(unsigned int k=0; k<getDimension(); k++) {
206 19 : normalization_[i][k] = getNormalization(kappas_[i][k],periods_[k]);
207 : }
208 : }
209 9 : checkRead();
210 9 : }
211 :
212 :
213 31100 : double TD_VonMises::getValue(const std::vector<double>& argument) const {
214 : double value=0.0;
215 92400 : for(unsigned int i=0; i<ncenters_; i++) {
216 61300 : value+=weights_[i]*VonMisesDiagonal(argument, centers_[i], kappas_[i],periods_,normalization_[i]);
217 : }
218 31100 : return value;
219 : }
220 :
221 :
222 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 {
223 : double value = 1.0;
224 220638 : for(unsigned int k=0; k<argument.size(); k++) {
225 140319 : double arg = kappa[k]*cos( ((2*pi)/periods[k])*(argument[k]-center[k]) );
226 140319 : value*=normalization[k]*exp(arg);
227 : }
228 80319 : return value;
229 : }
230 :
231 :
232 19 : double TD_VonMises::getNormalization(const double kappa, const double period) const {
233 : //
234 19 : std::vector<double> centers(1);
235 19 : centers[0] = 0.0;
236 19 : std::vector<double> kappas(1);
237 19 : kappas[0] = kappa;
238 19 : std::vector<double> periods(1);
239 19 : periods[0] = period;
240 19 : std::vector<double> norm(1);
241 19 : norm[0] = 1.0;
242 : //
243 : const unsigned int nbins = 1001;
244 19 : std::vector<double> points;
245 19 : std::vector<double> weights;
246 : double min = 0.0;
247 : double max = period;
248 19 : GridIntegrationWeights::getOneDimensionalIntegrationPointsAndWeights(points,weights,nbins,min,max);
249 : //
250 : double sum = 0.0;
251 19038 : for(unsigned int l=0; l<nbins; l++) {
252 19019 : std::vector<double> arg(1);
253 19019 : arg[0]= points[l];
254 19019 : sum += weights[l] * VonMisesDiagonal(arg,centers,kappas,periods,norm);
255 : }
256 38 : return 1.0/sum;
257 : }
258 :
259 :
260 : }
261 : }
|
