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 122 : 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 4882 : PLUMED_REGISTER_ACTION(TD_Uniform,"TD_UNIFORM")
187 :
188 :
189 62 : void TD_Uniform::registerKeywords(Keywords& keys) {
190 62 : TargetDistribution::registerKeywords(keys);
191 62 : 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 62 : 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 62 : 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 62 : 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 61 : parseVector("MINIMA",minima_);
206 61 : parseVector("MAXIMA",maxima_);
207 :
208 61 : parseVector("SIGMA_MINIMA",sigma_min_);
209 61 : parseVector("SIGMA_MAXIMA",sigma_max_);
210 61 : if(minima_.size()==0 && sigma_min_.size()>0) {plumed_merror(getName()+": you cannot give SIGMA_MINIMA if MINIMA is not given");}
211 61 : if(maxima_.size()==0 && sigma_max_.size()>0) {plumed_merror(getName()+": you cannot give SIGMA_MAXIMA if MAXIMA is not given");}
212 :
213 61 : 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 98 : for(unsigned int k=0; k<getDimension(); k++) {
218 52 : 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 15 : else if(minima_.size()>0 && maxima_.size()==0) {
224 : // only MINIMA given, MAXIMA assigned later on.
225 1 : setDimension(minima_.size());
226 : }
227 14 : else if(maxima_.size()>0 && minima_.size()==0) {
228 : // only MAXIMA given, MINIMA assigned later on.
229 1 : setDimension(maxima_.size());
230 : }
231 13 : 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 61 : if(sigma_min_.size()==0) {sigma_min_.assign(getDimension(),0.0);}
237 61 : 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 61 : 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 14 : minima_.assign(getDimension(),0.0);
250 14 : for(unsigned int k=0; k<getDimension(); k++) {Tools::convert(min[k],minima_[k]);}
251 : }
252 :
253 61 : if(maxima_.size()==0) {
254 14 : maxima_.assign(getDimension(),0.0);
255 14 : 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 115046 : double value = 1.0;
264 331503 : for(unsigned int k=0; k<getDimension(); k++) {
265 : double tmp;
266 216457 : 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 185668 : tmp = 1.0;
274 : }
275 216457 : value *= tmp;
276 : }
277 115046 : return value;
278 : }
279 :
280 : inline
281 30789 : 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 7511 : return 0.0;
288 : }
289 : }
290 :
291 :
292 :
293 :
294 :
295 :
296 : }
297 4821 : }
|
