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 : PLUMED_REGISTER_ACTION(TD_Uniform,"TD_UNIFORM")
195 :
196 :
197 75 : void TD_Uniform::registerKeywords(Keywords& keys) {
198 75 : TargetDistribution::registerKeywords(keys);
199 75 : 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 75 : 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 75 : 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 75 : 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 75 : }
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 73 : parseVector("MINIMA",minima_);
213 73 : parseVector("MAXIMA",maxima_);
214 :
215 73 : parseVector("SIGMA_MINIMA",sigma_min_);
216 146 : parseVector("SIGMA_MAXIMA",sigma_max_);
217 73 : if(minima_.size()==0 && sigma_min_.size()>0) {
218 0 : plumed_merror(getName()+": you cannot give SIGMA_MINIMA if MINIMA is not given");
219 : }
220 73 : if(maxima_.size()==0 && sigma_max_.size()>0) {
221 0 : plumed_merror(getName()+": you cannot give SIGMA_MAXIMA if MAXIMA is not given");
222 : }
223 :
224 73 : if(minima_.size()>0 && maxima_.size()>0) {
225 : // both MINIMA and MAXIMA given, do all checks
226 58 : if(minima_.size()!=maxima_.size()) {
227 0 : plumed_merror(getName()+": MINIMA and MAXIMA do not have the same number of values.");
228 : }
229 58 : setDimension(minima_.size());
230 122 : for(unsigned int k=0; k<getDimension(); k++) {
231 64 : if(minima_[k]>maxima_[k]) {
232 0 : plumed_merror(getName()+": error in MINIMA and MAXIMA keywords, one of the MINIMA values is larger than the corresponding MAXIMA values");
233 : }
234 : }
235 15 : } else if(minima_.size()>0 && maxima_.size()==0) {
236 : // only MINIMA given, MAXIMA assigned later on.
237 1 : setDimension(minima_.size());
238 14 : } else if(maxima_.size()>0 && minima_.size()==0) {
239 : // only MAXIMA given, MINIMA assigned later on.
240 1 : setDimension(maxima_.size());
241 13 : } else if(maxima_.size()==0 && minima_.size()==0) {
242 : // neither MAXIMA nor MINIMA givenm, both assigned later on.
243 13 : setDimension(0);
244 : }
245 :
246 73 : if(sigma_min_.size()==0) {
247 65 : sigma_min_.assign(getDimension(),0.0);
248 : }
249 73 : if(sigma_max_.size()==0) {
250 65 : sigma_max_.assign(getDimension(),0.0);
251 : }
252 73 : if(sigma_min_.size()!=getDimension()) {
253 0 : plumed_merror(getName()+": SIGMA_MINIMA has the wrong number of values");
254 : }
255 73 : if(sigma_max_.size()!=getDimension()) {
256 0 : plumed_merror(getName()+": SIGMA_MAXIMA has the wrong number of values");
257 : }
258 : //
259 : setForcedNormalization();
260 73 : checkRead();
261 73 : }
262 :
263 :
264 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) {
265 :
266 73 : if(minima_.size()==0) {
267 14 : minima_.assign(getDimension(),0.0);
268 33 : for(unsigned int k=0; k<getDimension(); k++) {
269 19 : Tools::convert(min[k],minima_[k]);
270 : }
271 : }
272 :
273 73 : if(maxima_.size()==0) {
274 14 : maxima_.assign(getDimension(),0.0);
275 33 : for(unsigned int k=0; k<getDimension(); k++) {
276 19 : Tools::convert(max[k],maxima_[k]);
277 : }
278 : }
279 :
280 73 : }
281 :
282 :
283 118654 : double TD_Uniform::getValue(const std::vector<double>& argument) const {
284 : //
285 : double value = 1.0;
286 338719 : for(unsigned int k=0; k<getDimension(); k++) {
287 : double tmp;
288 220065 : if(argument[k] < minima_[k]) {
289 15379 : tmp = GaussianSwitchingFunc(argument[k],minima_[k],sigma_min_[k]);
290 204686 : } else if(argument[k] > maxima_[k]) {
291 15566 : tmp = GaussianSwitchingFunc(argument[k],maxima_[k],sigma_max_[k]);
292 : } else {
293 : tmp = 1.0;
294 : }
295 220065 : value *= tmp;
296 : }
297 118654 : return value;
298 : }
299 :
300 : inline
301 : double TD_Uniform::GaussianSwitchingFunc(const double argument, const double center, const double sigma) const {
302 30945 : if(sigma>0.0) {
303 23278 : double arg=(argument-center)/sigma;
304 23278 : return exp(-0.5*arg*arg);
305 : } else {
306 : return 0.0;
307 : }
308 : }
309 :
310 :
311 :
312 :
313 :
314 :
315 : }
316 : }
|
