Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-2018 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 27 : 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;
109 : };
110 :
111 :
112 6461 : PLUMED_REGISTER_ACTION(TD_VonMises,"TD_VONMISES")
113 :
114 :
115 10 : void TD_VonMises::registerKeywords(Keywords& keys) {
116 10 : TargetDistribution::registerKeywords(keys);
117 40 : keys.add("numbered","CENTER","The centers of the Von Mises distributions.");
118 40 : keys.add("numbered","SIGMA","The \"standard deviations\" of the Von Mises distributions.");
119 40 : 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 40 : keys.add("hidden","PERIODS","The periods for each of the dimensions. By default they are 2*pi for each dimension.");
121 20 : keys.use("WELLTEMPERED_FACTOR");
122 20 : 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 : sigmas_(0),
130 : centers_(0),
131 : normalization_(0),
132 : weights_(0),
133 : 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 35 : for(unsigned int i=0; i<ncenters_; i++) {
157 26 : 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 57 : for(unsigned int i=0; i<sigmas_.size(); i++) {
163 13 : kappas_[i].resize(sigmas_[i].size());
164 83 : for(unsigned int k=0; k<kappas_[i].size(); k++) {
165 38 : kappas_[i][k] = 1.0/(sigmas_[i][k]*sigmas_[i][k]);
166 : }
167 : }
168 : //
169 18 : parseVector("WEIGHTS",weights_);
170 18 : 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 27 : 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 57 : for(unsigned int i=0; i<weights_.size(); i++) {sum_weights+=weights_[i];}
179 57 : for(unsigned int i=0; i<weights_.size(); i++) {weights_[i]/=sum_weights;}
180 : //
181 9 : normalization_.resize(ncenters_);
182 35 : for(unsigned int i=0; i<ncenters_; i++) {
183 26 : normalization_[i].resize(getDimension());
184 51 : for(unsigned int k=0; k<getDimension(); k++) {
185 76 : 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 153700 : for(unsigned int i=0; i<ncenters_; i++) {
195 245200 : 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 581595 : for(unsigned int k=0; k<argument.size(); k++) {
204 561276 : 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 38 : GridIntegrationWeights::getOneDimensionalIntegrationPointsAndWeights(points,weights,nbins,min,max);
228 : //
229 : double sum = 0.0;
230 38057 : for(unsigned int l=0; l<nbins; l++) {
231 38038 : 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 4839 : }
|
