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 :
25 : #include "core/ActionRegister.h"
26 :
27 :
28 : namespace PLMD {
29 : namespace ves {
30 :
31 : //+PLUMEDOC VES_TARGETDIST TD_GENERALIZED_NORMAL
32 : /*
33 : Target distribution given by a sum of generalized normal distributions (static).
34 :
35 : Employ a target distribution that is given by a sum where each
36 : term is a product of one-dimensional
37 : [generalized normal distributions](https://en.wikipedia.org/wiki/Generalized_normal_distribution)
38 : (version 1, also know as an exponential power distribution), defined as
39 : \f[
40 : p(\mathbf{s}) = \sum_{i} \, w_{i}
41 : \prod_{k}^{d}
42 : \frac{\beta_{k,i}}{2 \, \alpha_{k,i} \, \Gamma(1/\beta_{k,i})}
43 : \exp\left( -\left\vert \frac{s_{k}-\mu_{k,i}}{\alpha_{k,i}} \right\vert^{\beta_{k,i}} \right)
44 : \f]
45 : where \f$(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})\f$
46 : are the centers of the distributions,
47 : \f$(\alpha_{1,i},\alpha_{2,i},\ldots,\alpha_{d,i})\f$ are the scale
48 : parameters of the distributions,
49 : \f$(\beta_{1,i},\beta_{2,i},\ldots,\beta_{d,i})\f$ are the shape
50 : parameters of the distributions, and \f$\Gamma(x)\f$ is the
51 : gamma function.
52 : The weights \f$w_{i}\f$ are normalized to 1, \f$\sum_{i}w_{i}=1\f$.
53 :
54 : Employing \f$\beta=2\f$ results in a
55 : Gaussian (normal) distributions with mean
56 : \f$\mu\f$ and variance \f$\alpha^2/2\f$,
57 : \f$\beta=1\f$ gives the Laplace distribution, and
58 : the limit \f$\beta \to \infty\f$ results in a
59 : uniform distribution on the interval \f$[\mu-\alpha,\mu+\alpha]\f$.
60 :
61 : The centers \f$(\mu_{1,i},\mu_{2,i},\ldots,\mu_{d,i})\f$
62 : are given using the numbered CENTER keywords, the scale
63 : parameters \f$(\alpha_{1,i},\alpha_{2,i},\ldots,\alpha_{d,i})\f$
64 : using the numbered SCALE keywords, and the shape parameters
65 : \f$(\beta_{1,i},\beta_{2,i},\ldots,\beta_{d,i})\f$ using the
66 : numbered SHAPE keywords.
67 : The weights are given using the WEIGHTS keywords, if no weights are
68 : given are all terms weighted equally.
69 :
70 : \par Examples
71 :
72 : A generalized normal distribution in one-dimensional
73 : \plumedfile
74 : td1: TD_GENERALIZED_NORMAL CENTER1=+20.0 ALPHA1=5.0 BETA1=4.0
75 : \endplumedfile
76 :
77 : A sum of two one-dimensional generalized normal distributions
78 : \plumedfile
79 : TD_GENERALIZED_NORMAL ...
80 : CENTER1=+20.0 ALPHA1=5.0 BETA1=4.0
81 : CENTER2=-20.0 ALPHA2=5.0 BETA2=3.0
82 : LABEL=td1
83 : ... TD_GENERALIZED_NORMAL
84 : \endplumedfile
85 :
86 : A sum of two two-dimensional generalized normal distributions
87 : \plumedfile
88 : TD_GENERALIZED_NORMAL ...
89 : CENTER1=-20.0,-20.0 ALPHA1=5.0,3.0 BETA1=2.0,4.0
90 : CENTER2=-20.0,+20.0 ALPHA2=3.0,5.0 BETA2=4.0,2.0
91 : WEIGHTS=2.0,1.0
92 : LABEL=td1
93 : ... TD_GENERALIZED_NORMAL
94 : \endplumedfile
95 :
96 : */
97 : //+ENDPLUMEDOC
98 :
99 : class TD_GeneralizedNormal: public TargetDistribution {
100 : std::vector< std::vector<double> > centers_;
101 : std::vector< std::vector<double> > alphas_;
102 : std::vector< std::vector<double> > betas_;
103 : std::vector< std::vector<double> > normalization_;
104 : std::vector<double> weights_;
105 : unsigned int ncenters_;
106 : double ExponentialPowerDiagonal(const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&, const std::vector<double>&) const;
107 : public:
108 : static void registerKeywords(Keywords&);
109 : explicit TD_GeneralizedNormal(const ActionOptions& ao);
110 : double getValue(const std::vector<double>&) const override;
111 : };
112 :
113 :
114 : PLUMED_REGISTER_ACTION(TD_GeneralizedNormal,"TD_GENERALIZED_NORMAL")
115 :
116 :
117 10 : void TD_GeneralizedNormal::registerKeywords(Keywords& keys) {
118 10 : TargetDistribution::registerKeywords(keys);
119 10 : keys.add("numbered","CENTER","The center of each generalized normal distribution.");
120 10 : keys.add("numbered","ALPHA","The alpha parameters for each generalized normal distribution.");
121 10 : keys.add("numbered","BETA","The beta parameters for each generalized normal distribution.");
122 10 : keys.add("optional","WEIGHTS","The weights of the generalized normal distribution. By default all are weighted equally.");
123 10 : keys.use("WELLTEMPERED_FACTOR");
124 10 : keys.use("SHIFT_TO_ZERO");
125 10 : keys.use("NORMALIZE");
126 10 : }
127 :
128 :
129 8 : TD_GeneralizedNormal::TD_GeneralizedNormal(const ActionOptions& ao):
130 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
131 16 : centers_(0),
132 8 : alphas_(0),
133 8 : betas_(0),
134 8 : normalization_(0),
135 8 : weights_(0),
136 16 : ncenters_(0) {
137 15 : for(unsigned int i=1;; i++) {
138 : std::vector<double> tmp_center;
139 46 : if(!parseNumberedVector("CENTER",i,tmp_center) ) {
140 : break;
141 : }
142 15 : centers_.push_back(tmp_center);
143 15 : }
144 15 : for(unsigned int i=1;; i++) {
145 : std::vector<double> tmp_alpha;
146 46 : if(!parseNumberedVector("ALPHA",i,tmp_alpha) ) {
147 : break;
148 : }
149 35 : for(unsigned int k=0; k<tmp_alpha.size(); k++) {
150 20 : if(tmp_alpha[k]<=0.0) {
151 0 : plumed_merror(getName()+": the values given in ALPHA should be positive");
152 : }
153 : }
154 15 : alphas_.push_back(tmp_alpha);
155 15 : }
156 15 : for(unsigned int i=1;; i++) {
157 : std::vector<double> tmp_beta;
158 46 : if(!parseNumberedVector("BETA",i,tmp_beta) ) {
159 : break;
160 : }
161 35 : for(unsigned int k=0; k<tmp_beta.size(); k++) {
162 20 : if(tmp_beta[k]<=0.0) {
163 0 : plumed_merror(getName()+": the values given in BETA should be positive");
164 : }
165 : }
166 15 : betas_.push_back(tmp_beta);
167 15 : }
168 : //
169 8 : if(centers_.size()==0) {
170 0 : plumed_merror(getName()+": CENTER keywords seem to be missing. Note that numbered keywords start at CENTER1.");
171 : }
172 : //
173 8 : if(centers_.size()!=alphas_.size() || centers_.size()!=betas_.size() ) {
174 0 : plumed_merror(getName()+": there has to be an equal amount of CENTER, ALPHA, and BETA keywords");
175 : }
176 : //
177 8 : setDimension(centers_[0].size());
178 8 : ncenters_ = centers_.size();
179 : //
180 : // check centers and sigmas
181 23 : for(unsigned int i=0; i<ncenters_; i++) {
182 15 : if(centers_[i].size()!=getDimension()) {
183 0 : plumed_merror(getName()+": one of the CENTER keyword does not match the given dimension");
184 : }
185 15 : if(alphas_[i].size()!=getDimension()) {
186 0 : plumed_merror(getName()+": one of the ALPHA keyword does not match the given dimension");
187 : }
188 15 : if(betas_[i].size()!=getDimension()) {
189 0 : plumed_merror(getName()+": one of the BETA keyword does not match the given dimension");
190 : }
191 : }
192 : //
193 16 : parseVector("WEIGHTS",weights_);
194 8 : if(weights_.size()==0) {
195 4 : weights_.assign(centers_.size(),1.0);
196 : }
197 8 : if(centers_.size()!=weights_.size()) {
198 0 : plumed_merror(getName()+": there has to be as many weights given in WEIGHTS as numbered CENTER keywords");
199 : }
200 : //
201 : double sum_weights=0.0;
202 23 : for(unsigned int i=0; i<weights_.size(); i++) {
203 15 : sum_weights+=weights_[i];
204 : }
205 23 : for(unsigned int i=0; i<weights_.size(); i++) {
206 15 : weights_[i]/=sum_weights;
207 : }
208 : //
209 8 : normalization_.resize(ncenters_);
210 23 : for(unsigned int i=0; i<ncenters_; i++) {
211 15 : normalization_[i].resize(getDimension());
212 35 : for(unsigned int k=0; k<getDimension(); k++) {
213 20 : normalization_[i][k] = 0.5*betas_[i][k]/(alphas_[i][k]*tgamma(1.0/betas_[i][k]));
214 : }
215 : }
216 8 : checkRead();
217 8 : }
218 :
219 :
220 21608 : double TD_GeneralizedNormal::getValue(const std::vector<double>& argument) const {
221 : double value=0.0;
222 74623 : for(unsigned int i=0; i<ncenters_; i++) {
223 53015 : value+=weights_[i]*ExponentialPowerDiagonal(argument,centers_[i],alphas_[i],betas_[i],normalization_[i]);
224 : }
225 21608 : return value;
226 : }
227 :
228 :
229 53015 : double TD_GeneralizedNormal::ExponentialPowerDiagonal(const std::vector<double>& argument, const std::vector<double>& center, const std::vector<double>& alpha, const std::vector<double>& beta, const std::vector<double>& normalization) const {
230 : double value = 1.0;
231 157035 : for(unsigned int k=0; k<argument.size(); k++) {
232 104020 : double arg=(std::abs(argument[k]-center[k]))/alpha[k];
233 104020 : arg = pow(arg,beta[k]);
234 104020 : value*=normalization[k]*exp(-arg);
235 : }
236 53015 : return value;
237 : }
238 :
239 :
240 :
241 : }
242 : }
|
