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_UNIFORM
32 : /*
33 : Uniform target distribution (static).
34 :
35 : Using this keyword you can define a uniform target distribution which is a
36 : product of one-dimensional distributions \f$p_{k}(s_{k})\f$ that are uniform
37 : over a given interval \f$[a_{k},b_{k}]\f$
38 :
39 : \f[
40 : p_{k}(s_{k}) =
41 : \left \{\begin{array}{ll}
42 : \frac{1}{(b_{k}-a_{k})} & \mathrm{if} \ a_{k} \leq s_{k} \leq b_{k} \\
43 : &\\
44 : 0 & \mathrm{otherwise}
45 : \end{array}\right .
46 : \f]
47 :
48 : The overall distribution is then given as
49 : \f[
50 : p(\mathbf{s}) =
51 : \prod^{d}_{k} p_{k}(s_{k}) =
52 : \left\{\begin{array}{ll}
53 : \prod^{d}_{k} \frac{1}{(b_{k}-a_{k})}
54 : & \mathrm{if} \ a_{k} \leq s_{k} \leq b_{k} \ \mathrm{for\ all}\ k \\
55 : \\
56 : 0 & \mathrm{otherwise}
57 : \end{array}\right.
58 : \f]
59 : The distribution is thus uniform inside a rectangular for two arguments
60 : and a cube for a three arguments.
61 :
62 : The limits of the intervals \f$ a_{k}\f$ and \f$ b_{k}\f$ are given
63 : with the MINIMA and MAXIMA keywords, respectively. If one or both of
64 : these keywords are missing the code should automatically detect the limits.
65 :
66 :
67 : It is also possible to use one-dimensional distributions
68 : that go smoothly to zero at the boundaries.
69 : This is done by employing a function with
70 : Gaussian switching functions at the boundaries \f$a_{k}\f$ and \f$b_{k}\f$
71 : \f[
72 : f_{k}(s_{k}) =
73 : \begin{cases}
74 : \exp\left(-\frac{(s_{k}-a_{k})^2}{2 \sigma^2_{a,k}}\right)
75 : & \mathrm{if}\, s_{k} < a_{k} \\
76 : \\
77 : 1 & \mathrm{if}\, a_{k} \leq s_{k} \leq b_{k} \\
78 : \\
79 : \exp\left(-\frac{(s_{k}-b_{k})^2}{2 \sigma^2_{b,k}}\right)
80 : & \mathrm{if}\, s_{k} > b_{k}
81 : \end{cases}
82 : \f]
83 : where the standard deviation parameters \f$\sigma_{a,k}\f$
84 : and \f$\sigma_{b,k}\f$ determine how quickly the switching functions
85 : goes to zero.
86 : The overall distribution is then normalized
87 : \f[
88 : p(\mathbf{s}) =
89 : \prod^{d}_{k} p_{k}(s_{k}) =
90 : \prod^{d}_{k} \frac{f(s_{k})}{\int d s_{k} \, f(s_{k})}
91 : \f]
92 : To use this option you need to provide the standard deviation
93 : parameters \f$\sigma_{a,k}\f$ and \f$\sigma_{b,k}\f$ by using the
94 : SIGMA_MINIMA and SIGMA_MAXIMA keywords, respectively. Giving a value of
95 : 0.0 means that the boundary is sharp, which is the default behavior.
96 :
97 :
98 :
99 :
100 :
101 :
102 : \par Examples
103 :
104 : If one or both of the MINIMA or MAXIMA keywords are missing
105 : the code should automatically detect the limits not given.
106 : Therefore, if we consider a target distribution that is
107 : defined over an interval from 0.0 to 10.0 for the first
108 : argument and from 0.2 to 1.0 for the second argument are
109 : the following example
110 : \plumedfile
111 : td: TD_UNIFORM
112 : \endplumedfile
113 :
114 : is equivalent to this one
115 :
116 : \plumedfile
117 : TD_UNIFORM ...
118 : MINIMA=0.0,0.2
119 : MAXIMA=10.0,1.0
120 : LABEL=td
121 : ... TD_UNIFORM
122 : \endplumedfile
123 :
124 : and this one
125 :
126 : \plumedfile
127 : td: TD_UNIFORM MAXIMA=10.0,1.0
128 : \endplumedfile
129 :
130 : and also this one
131 :
132 : \plumedfile
133 : td: TD_UNIFORM MINIMA=0.0,0,2
134 : \endplumedfile
135 :
136 :
137 : We can also define a target distribution that goes smoothly to zero
138 : at the boundaries of the uniform distribution. In the following
139 : we consider an interval of 0 to 10 for the target distribution.
140 : The following input would result in a target distribution that
141 : would be uniform from 2 to 7 and then smoothly go to zero from
142 : 2 to 0 and from 7 to 10.
143 : \plumedfile
144 : TD_UNIFORM ...
145 : MINIMA=2.0
146 : MAXIMA=+7.0
147 : SIGMA_MINIMA=0.5
148 : SIGMA_MAXIMA=1.0
149 : LABEL=td
150 : ... TD_UNIFORM
151 : \endplumedfile
152 : It is also possible to employ a smooth switching function for just one
153 : of the boundaries as shown here where the target distribution
154 : would be uniform from 0 to 7 and then smoothly go to zero from 7 to 10.
155 : \plumedfile
156 : TD_UNIFORM ...
157 : MAXIMA=+7.0
158 : SIGMA_MAXIMA=1.0
159 : LABEL=td
160 : ... TD_UNIFORM
161 : \endplumedfile
162 : Furthermore, it is possible to employ a sharp boundary by
163 : using
164 : \plumedfile
165 : TD_UNIFORM ...
166 : MAXIMA=+7.0
167 : SIGMA_MAXIMA=0.0
168 : LABEL=td
169 : ... TD_UNIFORM
170 : \endplumedfile
171 : or
172 : \plumedfile
173 : td: TD_UNIFORM MAXIMA=+7.0
174 : \endplumedfile
175 :
176 :
177 : */
178 : //+ENDPLUMEDOC
179 :
180 : class TD_Uniform : public TargetDistribution {
181 : std::vector<double> minima_;
182 : std::vector<double> maxima_;
183 : std::vector<double> sigma_min_;
184 : std::vector<double> sigma_max_;
185 : double GaussianSwitchingFunc(const double, const double, const double) const;
186 : void setupAdditionalGrids(const std::vector<Value*>&, const std::vector<std::string>&, const std::vector<std::string>&, const std::vector<unsigned int>&) override;
187 : public:
188 : static void registerKeywords( Keywords&);
189 : explicit TD_Uniform(const ActionOptions& ao);
190 : double getValue(const std::vector<double>&) const override;
191 : };
192 :
193 :
194 10492 : PLUMED_REGISTER_ACTION(TD_Uniform,"TD_UNIFORM")
195 :
196 :
197 74 : void TD_Uniform::registerKeywords(Keywords& keys) {
198 74 : TargetDistribution::registerKeywords(keys);
199 148 : keys.add("optional","MINIMA","The minimum of the intervals where the target distribution is taken as uniform. You should give one value for each argument.");
200 148 : keys.add("optional","MAXIMA","The maximum of the intervals where the target distribution is taken as uniform. You should give one value for each argument.");
201 148 : keys.add("optional","SIGMA_MINIMA","The standard deviation parameters of the Gaussian switching functions for the minima of the intervals. You should give one value for each argument. Value of 0.0 means that switch is done without a smooth switching function, this is the default behavior.");
202 148 : keys.add("optional","SIGMA_MAXIMA","The standard deviation parameters of the Gaussian switching functions for the maximum of the intervals. You should give one value for each argument. Value of 0.0 means that switch is done without a smooth switching function, this is the default behavior.");
203 74 : }
204 :
205 :
206 73 : TD_Uniform::TD_Uniform(const ActionOptions& ao):
207 : PLUMED_VES_TARGETDISTRIBUTION_INIT(ao),
208 146 : minima_(0),
209 73 : maxima_(0),
210 73 : sigma_min_(0),
211 146 : sigma_max_(0)
212 : {
213 73 : parseVector("MINIMA",minima_);
214 73 : parseVector("MAXIMA",maxima_);
215 :
216 73 : parseVector("SIGMA_MINIMA",sigma_min_);
217 146 : parseVector("SIGMA_MAXIMA",sigma_max_);
218 73 : if(minima_.size()==0 && sigma_min_.size()>0) {plumed_merror(getName()+": you cannot give SIGMA_MINIMA if MINIMA is not given");}
219 73 : if(maxima_.size()==0 && sigma_max_.size()>0) {plumed_merror(getName()+": you cannot give SIGMA_MAXIMA if MAXIMA is not given");}
220 :
221 73 : if(minima_.size()>0 && maxima_.size()>0) {
222 : // both MINIMA and MAXIMA given, do all checks
223 58 : if(minima_.size()!=maxima_.size()) {plumed_merror(getName()+": MINIMA and MAXIMA do not have the same number of values.");}
224 58 : setDimension(minima_.size());
225 122 : for(unsigned int k=0; k<getDimension(); k++) {
226 64 : if(minima_[k]>maxima_[k]) {
227 0 : plumed_merror(getName()+": error in MINIMA and MAXIMA keywords, one of the MINIMA values is larger than the corresponding MAXIMA values");
228 : }
229 : }
230 : }
231 15 : else if(minima_.size()>0 && maxima_.size()==0) {
232 : // only MINIMA given, MAXIMA assigned later on.
233 1 : setDimension(minima_.size());
234 : }
235 14 : else if(maxima_.size()>0 && minima_.size()==0) {
236 : // only MAXIMA given, MINIMA assigned later on.
237 1 : setDimension(maxima_.size());
238 : }
239 13 : else if(maxima_.size()==0 && minima_.size()==0) {
240 : // neither MAXIMA nor MINIMA givenm, both assigned later on.
241 13 : setDimension(0);
242 : }
243 :
244 73 : if(sigma_min_.size()==0) {sigma_min_.assign(getDimension(),0.0);}
245 73 : if(sigma_max_.size()==0) {sigma_max_.assign(getDimension(),0.0);}
246 73 : if(sigma_min_.size()!=getDimension()) {plumed_merror(getName()+": SIGMA_MINIMA has the wrong number of values");}
247 73 : if(sigma_max_.size()!=getDimension()) {plumed_merror(getName()+": SIGMA_MAXIMA has the wrong number of values");}
248 : //
249 : setForcedNormalization();
250 73 : checkRead();
251 73 : }
252 :
253 :
254 73 : void TD_Uniform::setupAdditionalGrids(const std::vector<Value*>& arguments, const std::vector<std::string>& min, const std::vector<std::string>& max, const std::vector<unsigned int>& nbins) {
255 :
256 73 : if(minima_.size()==0) {
257 14 : minima_.assign(getDimension(),0.0);
258 33 : for(unsigned int k=0; k<getDimension(); k++) {Tools::convert(min[k],minima_[k]);}
259 : }
260 :
261 73 : if(maxima_.size()==0) {
262 14 : maxima_.assign(getDimension(),0.0);
263 33 : for(unsigned int k=0; k<getDimension(); k++) {Tools::convert(max[k],maxima_[k]);}
264 : }
265 :
266 73 : }
267 :
268 :
269 118654 : double TD_Uniform::getValue(const std::vector<double>& argument) const {
270 : //
271 : double value = 1.0;
272 338719 : for(unsigned int k=0; k<getDimension(); k++) {
273 : double tmp;
274 220065 : if(argument[k] < minima_[k]) {
275 15379 : tmp = GaussianSwitchingFunc(argument[k],minima_[k],sigma_min_[k]);
276 : }
277 204686 : else if(argument[k] > maxima_[k]) {
278 15566 : tmp = GaussianSwitchingFunc(argument[k],maxima_[k],sigma_max_[k]);
279 : }
280 : else {
281 : tmp = 1.0;
282 : }
283 220065 : value *= tmp;
284 : }
285 118654 : return value;
286 : }
287 :
288 : inline
289 : double TD_Uniform::GaussianSwitchingFunc(const double argument, const double center, const double sigma) const {
290 30945 : if(sigma>0.0) {
291 23278 : double arg=(argument-center)/sigma;
292 23278 : return exp(-0.5*arg*arg);
293 : }
294 : else {
295 : return 0.0;
296 : }
297 : }
298 :
299 :
300 :
301 :
302 :
303 :
304 : }
305 : }
|
