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