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