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