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