Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed 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 : plumed 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 plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "HistogramBead.h"
23 : #include <vector>
24 : #include <limits>
25 : #include "Tools.h"
26 : #include "Keywords.h"
27 :
28 : namespace PLMD {
29 :
30 : //+PLUMEDOC INTERNAL histogrambead
31 : /*
32 : A function that can be used to calculate whether quantities are between fixed upper and lower bounds.
33 :
34 : If we have multiple instances of a variable we can estimate the probability density function
35 : for that variable using a process called kernel density estimation:
36 :
37 : \f[
38 : P(s) = \sum_i K\left( \frac{s - s_i}{w} \right)
39 : \f]
40 :
41 : In this equation \f$K\f$ is a symmetric function that must integrate to one that is often
42 : called a kernel function and \f$w\f$ is a smearing parameter. From a probability density function calculated using
43 : kernel density estimation we can calculate the number/fraction of values between an upper and lower
44 : bound using:
45 :
46 : \f[
47 : w(s) = \int_a^b \sum_i K\left( \frac{s - s_i}{w} \right)
48 : \f]
49 :
50 : All the input to calculate a quantity like \f$w(s)\f$ is generally provided through a single
51 : keyword that will have the following form:
52 :
53 : KEYWORD={TYPE UPPER=\f$a\f$ LOWER=\f$b\f$ SMEAR=\f$\frac{w}{b-a}\f$}
54 :
55 : This will calculate the number of values between \f$a\f$ and \f$b\f$. To calculate
56 : the fraction of values you add the word NORM to the input specification. If the
57 : function keyword SMEAR is not present \f$w\f$ is set equal to \f$0.5(b-a)\f$. Finally,
58 : type should specify one of the kernel types that is present in plumed. These are listed
59 : in the table below:
60 :
61 : <table align=center frame=void width=95%% cellpadding=5%%>
62 : <tr>
63 : <td> TYPE </td> <td> FUNCTION </td>
64 : </tr> <tr>
65 : <td> GAUSSIAN </td> <td> \f$\frac{1}{\sqrt{2\pi}w} \exp\left( -\frac{(s-s_i)^2}{2w^2} \right)\f$ </td>
66 : </tr> <tr>
67 : <td> TRIANGULAR </td> <td> \f$ \frac{1}{2w} \left( 1. - \left| \frac{s-s_i}{w} \right| \right) \quad \frac{s-s_i}{w}<1 \f$ </td>
68 : </tr>
69 : </table>
70 :
71 : Some keywords can also be used to calculate a discrete version of the histogram. That
72 : is to say the number of values between \f$a\f$ and \f$b\f$, the number of values between
73 : \f$b\f$ and \f$c\f$ and so on. A keyword that specifies this sort of calculation would look
74 : something like
75 :
76 : KEYWORD={TYPE UPPER=\f$a\f$ LOWER=\f$b\f$ NBINS=\f$n\f$ SMEAR=\f$\frac{w}{n(b-a)}\f$}
77 :
78 : This specification would calculate the following vector of quantities:
79 :
80 : \f[
81 : w_j(s) = \int_{a + \frac{j-1}{n}(b-a)}^{a + \frac{j}{n}(b-a)} \sum_i K\left( \frac{s - s_i}{w} \right)
82 : \f]
83 :
84 : */
85 : //+ENDPLUMEDOC
86 :
87 0 : void HistogramBead::registerKeywords( Keywords& keys ) {
88 0 : keys.add("compulsory","LOWER","the lower boundary for this particular bin");
89 0 : keys.add("compulsory","UPPER","the upper boundary for this particular bin");
90 0 : keys.add("compulsory","SMEAR","0.5","the amount to smear the Gaussian for each value in the distribution");
91 0 : }
92 :
93 260467 : HistogramBead::HistogramBead():
94 260467 : init(false),
95 260467 : lowb(0.0),
96 260467 : highb(0.0),
97 260467 : width(0.0),
98 260467 : cutoff(std::numeric_limits<double>::max()),
99 260467 : type(gaussian),
100 260467 : periodicity(unset),
101 260467 : min(0.0),
102 260467 : max(0.0),
103 260467 : max_minus_min(0.0),
104 260467 : inv_max_minus_min(0.0) {
105 260467 : }
106 :
107 53 : std::string HistogramBead::description() const {
108 53 : std::ostringstream ostr;
109 53 : ostr<<"between "<<lowb<<" and "<<highb<<" width of gaussian window equals "<<width;
110 53 : return ostr.str();
111 53 : }
112 :
113 0 : void HistogramBead::generateBins( const std::string& params, std::vector<std::string>& bins ) {
114 0 : std::vector<std::string> data=Tools::getWords(params);
115 0 : plumed_massert(data.size()>=1,"There is no input for this keyword");
116 :
117 0 : std::string name=data[0];
118 :
119 : unsigned nbins;
120 0 : std::vector<double> range(2);
121 : std::string smear;
122 0 : bool found_nb=Tools::parse(data,"NBINS",nbins);
123 0 : plumed_massert(found_nb,"Number of bins in histogram not found");
124 0 : bool found_r=Tools::parse(data,"LOWER",range[0]);
125 0 : plumed_massert(found_r,"Lower bound for histogram not specified");
126 0 : found_r=Tools::parse(data,"UPPER",range[1]);
127 0 : plumed_massert(found_r,"Upper bound for histogram not specified");
128 0 : plumed_massert(range[0]<range[1],"Range specification is dubious");
129 0 : bool found_b=Tools::parse(data,"SMEAR",smear);
130 0 : if(!found_b) {
131 0 : Tools::convert(0.5,smear);
132 : }
133 :
134 : std::string lb,ub;
135 0 : double delr = ( range[1]-range[0] ) / static_cast<double>( nbins );
136 0 : for(unsigned i=0; i<nbins; ++i) {
137 0 : Tools::convert( range[0]+i*delr, lb );
138 0 : Tools::convert( range[0]+(i+1)*delr, ub );
139 0 : bins.push_back( name + " " + "LOWER=" + lb + " " + "UPPER=" + ub + " " + "SMEAR=" + smear );
140 : }
141 0 : plumed_assert(bins.size()==nbins);
142 0 : }
143 :
144 53 : void HistogramBead::set( const std::string& params, std::string& errormsg ) {
145 53 : std::vector<std::string> data=Tools::getWords(params);
146 53 : if(data.size()<1) {
147 : errormsg="No input has been specified";
148 : return;
149 : }
150 :
151 53 : std::string name=data[0];
152 : const double DP2CUTOFF=6.25;
153 53 : if(name=="GAUSSIAN") {
154 53 : type=gaussian;
155 53 : cutoff=std::sqrt(2.0*DP2CUTOFF);
156 0 : } else if(name=="TRIANGULAR") {
157 0 : type=triangular;
158 0 : cutoff=1.;
159 : } else {
160 0 : plumed_merror("cannot understand kernel type " + name );
161 : }
162 :
163 : double smear;
164 53 : bool found_r=Tools::parse(data,"LOWER",lowb);
165 53 : if( !found_r ) {
166 : errormsg="Lower bound has not been specified use LOWER";
167 : }
168 53 : found_r=Tools::parse(data,"UPPER",highb);
169 53 : if( !found_r ) {
170 : errormsg="Upper bound has not been specified use UPPER";
171 : }
172 53 : if( lowb>=highb ) {
173 : errormsg="Lower bound is higher than upper bound";
174 : }
175 :
176 53 : smear=0.5;
177 53 : Tools::parse(data,"SMEAR",smear);
178 53 : width=smear*(highb-lowb);
179 53 : init=true;
180 53 : }
181 :
182 2668474 : void HistogramBead::set( double l, double h, double w) {
183 2668474 : init=true;
184 2668474 : lowb=l;
185 2668474 : highb=h;
186 2668474 : width=w;
187 : const double DP2CUTOFF=6.25;
188 2668474 : if( type==gaussian ) {
189 733122 : cutoff=std::sqrt(2.0*DP2CUTOFF);
190 1935352 : } else if( type==triangular ) {
191 1935352 : cutoff=1.;
192 : } else {
193 0 : plumed_error();
194 : }
195 2668474 : }
196 :
197 260148 : void HistogramBead::setKernelType( const std::string& ktype ) {
198 260148 : if(ktype=="gaussian") {
199 233896 : type=gaussian;
200 26252 : } else if(ktype=="triangular") {
201 26252 : type=triangular;
202 : } else {
203 0 : plumed_merror("cannot understand kernel type " + ktype );
204 : }
205 260148 : }
206 :
207 251184 : double HistogramBead::calculate( double x, double& df ) const {
208 : plumed_dbg_assert(init && periodicity!=unset );
209 : double lowB, upperB, f;
210 251184 : if( type==gaussian ) {
211 251184 : lowB = difference( x, lowb ) / ( std::sqrt(2.0) * width );
212 251184 : upperB = difference( x, highb ) / ( std::sqrt(2.0) * width );
213 251184 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( std::sqrt(2*pi)*width );
214 251184 : f = 0.5*( erf( upperB ) - erf( lowB ) );
215 0 : } else if( type==triangular ) {
216 0 : lowB = ( difference( x, lowb ) / width );
217 0 : upperB = ( difference( x, highb ) / width );
218 0 : df=0;
219 0 : if( std::fabs(lowB)<1. ) {
220 0 : df = (1 - std::fabs(lowB)) / width;
221 : }
222 0 : if( std::fabs(upperB)<1. ) {
223 0 : df -= (1 - std::fabs(upperB)) / width;
224 : }
225 0 : if (upperB<=-1. || lowB >=1.) {
226 : f=0.;
227 : } else {
228 : double ia, ib;
229 0 : if( lowB>-1.0 ) {
230 : ia=lowB;
231 : } else {
232 : ia=-1.0;
233 : }
234 0 : if( upperB<1.0 ) {
235 : ib=upperB;
236 : } else {
237 : ib=1.0;
238 : }
239 0 : f = (ib*(2.-std::fabs(ib))-ia*(2.-std::fabs(ia)))*0.5;
240 : }
241 : } else {
242 0 : plumed_merror("function type does not exist");
243 : }
244 251184 : return f;
245 : }
246 :
247 920596 : double HistogramBead::calculateWithCutoff( double x, double& df ) const {
248 : plumed_dbg_assert(init && periodicity!=unset );
249 :
250 : double lowB, upperB, f;
251 920596 : lowB = difference( x, lowb ) / width ;
252 920596 : upperB = difference( x, highb ) / width;
253 920596 : if( upperB<=-cutoff || lowB>=cutoff ) {
254 87679 : df=0;
255 87679 : return 0;
256 : }
257 :
258 832917 : if( type==gaussian ) {
259 446395 : lowB /= std::sqrt(2.0);
260 446395 : upperB /= std::sqrt(2.0);
261 446395 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( std::sqrt(2*pi)*width );
262 446395 : f = 0.5*( erf( upperB ) - erf( lowB ) );
263 386522 : } else if( type==triangular ) {
264 386522 : df=0;
265 386522 : if( std::fabs(lowB)<1. ) {
266 98488 : df = (1 - std::fabs(lowB)) / width;
267 : }
268 386522 : if( std::fabs(upperB)<1. ) {
269 98488 : df -= (1 - std::fabs(upperB)) / width;
270 : }
271 386522 : if (upperB<=-1. || lowB >=1.) {
272 : f=0.;
273 : } else {
274 : double ia, ib;
275 386522 : if( lowB>-1.0 ) {
276 : ia=lowB;
277 : } else {
278 : ia=-1.0;
279 : }
280 386522 : if( upperB<1.0 ) {
281 : ib=upperB;
282 : } else {
283 : ib=1.0;
284 : }
285 386522 : f = (ib*(2.-std::fabs(ib))-ia*(2.-std::fabs(ia)))*0.5;
286 : }
287 : } else {
288 0 : plumed_merror("function type does not exist");
289 : }
290 : return f;
291 : }
292 :
293 4680 : double HistogramBead::lboundDerivative( const double& x ) const {
294 4680 : if( type==gaussian ) {
295 4680 : double lowB = difference( x, lowb ) / ( std::sqrt(2.0) * width );
296 4680 : return exp( -lowB*lowB ) / ( std::sqrt(2*pi)*width );
297 0 : } else if ( type==triangular ) {
298 0 : plumed_error();
299 : // lowB = fabs( difference( x, lowb ) / width );
300 : // if( lowB<1 ) return ( 1 - (lowB) ) / 2*width;
301 : // else return 0;
302 : } else {
303 0 : plumed_merror("function type does not exist");
304 : }
305 : return 0;
306 : }
307 :
308 4680 : double HistogramBead::uboundDerivative( const double& x ) const {
309 : plumed_dbg_assert(init && periodicity!=unset );
310 4680 : if( type==gaussian ) {
311 4680 : double upperB = difference( x, highb ) / ( std::sqrt(2.0) * width );
312 4680 : return exp( -upperB*upperB ) / ( std::sqrt(2*pi)*width );
313 0 : } else if ( type==triangular ) {
314 0 : plumed_error();
315 : // upperB = fabs( difference( x, highb ) / width );
316 : // if( upperB<1 ) return ( 1 - (upperB) ) / 2*width;
317 : // else return 0;
318 : } else {
319 0 : plumed_merror("function type does not exist");
320 : }
321 : return 0;
322 : }
323 :
324 : }
|
