Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2018 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 distribution (pdf)
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 funciton that must integrate to one that is often
42 : called a kernel function and \f$w\f$ is a smearing parameter. From a pdf 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 descretized 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 72 : void HistogramBead::registerKeywords( Keywords& keys ) {
88 72 : keys.add("compulsory","LOWER","the lower boundary for this particular bin");
89 72 : keys.add("compulsory","UPPER","the upper boundary for this particular bin");
90 72 : keys.add("compulsory","SMEAR","0.5","the ammount to smear the Gaussian for each value in the distribution");
91 72 : }
92 :
93 69588 : HistogramBead::HistogramBead():
94 : init(false),
95 : lowb(0.0),
96 : highb(0.0),
97 : width(0.0),
98 69588 : cutoff(std::numeric_limits<double>::max()),
99 : type(gaussian),
100 : periodicity(unset),
101 : min(0.0),
102 : max(0.0),
103 : max_minus_min(0.0),
104 139207 : inv_max_minus_min(0.0)
105 : {
106 69619 : }
107 :
108 65 : std::string HistogramBead::description() const {
109 65 : std::ostringstream ostr;
110 65 : ostr<<"betweeen "<<lowb<<" and "<<highb<<" width of gaussian window equals "<<width;
111 65 : return ostr.str();
112 : }
113 :
114 10 : void HistogramBead::generateBins( const std::string& params, std::vector<std::string>& bins ) {
115 10 : std::vector<std::string> data=Tools::getWords(params);
116 10 : plumed_massert(data.size()>=1,"There is no input for this keyword");
117 :
118 20 : std::string name=data[0];
119 :
120 20 : unsigned nbins; std::vector<double> range(2); std::string smear;
121 10 : bool found_nb=Tools::parse(data,"NBINS",nbins);
122 10 : plumed_massert(found_nb,"Number of bins in histogram not found");
123 10 : bool found_r=Tools::parse(data,"LOWER",range[0]);
124 10 : plumed_massert(found_r,"Lower bound for histogram not specified");
125 10 : found_r=Tools::parse(data,"UPPER",range[1]);
126 10 : plumed_massert(found_r,"Upper bound for histogram not specified");
127 10 : plumed_massert(range[0]<range[1],"Range specification is dubious");
128 10 : bool found_b=Tools::parse(data,"SMEAR",smear);
129 10 : if(!found_b) { Tools::convert(0.5,smear); }
130 :
131 20 : std::string lb,ub; double delr = ( range[1]-range[0] ) / static_cast<double>( nbins );
132 40 : for(unsigned i=0; i<nbins; ++i) {
133 30 : Tools::convert( range[0]+i*delr, lb );
134 30 : Tools::convert( range[0]+(i+1)*delr, ub );
135 30 : bins.push_back( name + " " + "LOWER=" + lb + " " + "UPPER=" + ub + " " + "SMEAR=" + smear );
136 : }
137 20 : plumed_assert(bins.size()==nbins);
138 10 : }
139 :
140 65 : void HistogramBead::set( const std::string& params, std::string& errormsg ) {
141 65 : std::vector<std::string> data=Tools::getWords(params);
142 65 : if(data.size()<1) errormsg="No input has been specified";
143 :
144 130 : std::string name=data[0]; const double DP2CUTOFF=6.25;
145 65 : if(name=="GAUSSIAN") { type=gaussian; cutoff=sqrt(2.0*DP2CUTOFF); }
146 0 : else if(name=="TRIANGULAR") { type=triangular; cutoff=1.; }
147 0 : else plumed_merror("cannot understand kernel type " + name );
148 :
149 : double smear;
150 65 : bool found_r=Tools::parse(data,"LOWER",lowb);
151 65 : if( !found_r ) errormsg="Lower bound has not been specified use LOWER";
152 65 : found_r=Tools::parse(data,"UPPER",highb);
153 65 : if( !found_r ) errormsg="Upper bound has not been specified use UPPER";
154 65 : if( lowb>=highb ) errormsg="Lower bound is higher than upper bound";
155 :
156 65 : smear=0.5; Tools::parse(data,"SMEAR",smear);
157 130 : width=smear*(highb-lowb); init=true;
158 65 : }
159 :
160 147499 : void HistogramBead::set( double l, double h, double w) {
161 147499 : init=true; lowb=l; highb=h; width=w;
162 147499 : const double DP2CUTOFF=6.25;
163 147499 : if( type==gaussian ) cutoff=sqrt(2.0*DP2CUTOFF);
164 0 : else if( type==triangular ) cutoff=1.;
165 0 : else plumed_error();
166 147499 : }
167 :
168 69527 : void HistogramBead::setKernelType( const std::string& ktype ) {
169 69527 : if(ktype=="gaussian") type=gaussian;
170 0 : else if(ktype=="triangular") type=triangular;
171 0 : else plumed_merror("cannot understand kernel type " + ktype );
172 69578 : }
173 :
174 160424 : double HistogramBead::calculate( double x, double& df ) const {
175 : plumed_dbg_assert(init && periodicity!=unset );
176 : double lowB, upperB, f;
177 160424 : if( type==gaussian ) {
178 160424 : lowB = difference( x, lowb ) / ( sqrt(2.0) * width );
179 160404 : upperB = difference( x, highb ) / ( sqrt(2.0) * width );
180 160427 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( sqrt(2*pi)*width );
181 160427 : f = 0.5*( erf( upperB ) - erf( lowB ) );
182 0 : } else if( type==triangular ) {
183 0 : lowB = ( difference( x, lowb ) / width );
184 0 : upperB = ( difference( x, highb ) / width );
185 0 : df=0;
186 0 : if( fabs(lowB)<1. ) df = (1 - fabs(lowB)) / width;
187 0 : if( fabs(upperB)<1. ) df -= (1 - fabs(upperB)) / width;
188 0 : if (upperB<=-1. || lowB >=1.) {
189 0 : f=0.;
190 : } else {
191 : double ia, ib;
192 0 : if( lowB>-1.0 ) { ia=lowB; } else { ia=-1.0; }
193 0 : if( upperB<1.0 ) { ib=upperB; } else { ib=1.0; }
194 0 : f = (ib*(2.-fabs(ib))-ia*(2.-fabs(ia)))*0.5;
195 : }
196 : } else {
197 0 : plumed_merror("function type does not exist");
198 : }
199 160483 : return f;
200 : }
201 :
202 0 : double HistogramBead::calculateWithCutoff( double x, double& df ) const {
203 : plumed_dbg_assert(init && periodicity!=unset );
204 :
205 : double lowB, upperB, f;
206 0 : lowB = difference( x, lowb ) / width ; upperB = difference( x, highb ) / width;
207 0 : if( upperB<=-cutoff || lowB>=cutoff ) { df=0; return 0; }
208 :
209 0 : if( type==gaussian ) {
210 0 : lowB /= sqrt(2.0); upperB /= sqrt(2.0);
211 0 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( sqrt(2*pi)*width );
212 0 : f = 0.5*( erf( upperB ) - erf( lowB ) );
213 0 : } else if( type==triangular ) {
214 0 : df=0;
215 0 : if( fabs(lowB)<1. ) df = (1 - fabs(lowB)) / width;
216 0 : if( fabs(upperB)<1. ) df -= (1 - fabs(upperB)) / width;
217 0 : if (upperB<=-1. || lowB >=1.) {
218 0 : f=0.;
219 : } else {
220 : double ia, ib;
221 0 : if( lowB>-1.0 ) { ia=lowB; } else { ia=-1.0; }
222 0 : if( upperB<1.0 ) { ib=upperB; } else { ib=1.0; }
223 0 : f = (ib*(2.-fabs(ib))-ia*(2.-fabs(ia)))*0.5;
224 : }
225 : } else {
226 0 : plumed_merror("function type does not exist");
227 : }
228 0 : return f;
229 : }
230 :
231 4860 : double HistogramBead::lboundDerivative( const double& x ) const {
232 : double lowB;
233 4860 : if( type==gaussian ) {
234 4860 : lowB = difference( x, lowb ) / ( sqrt(2.0) * width );
235 4860 : return exp( -lowB*lowB ) / ( sqrt(2*pi)*width );
236 0 : } else if ( type==triangular ) {
237 0 : plumed_error();
238 : // lowB = fabs( difference( x, lowb ) / width );
239 : // if( lowB<1 ) return ( 1 - (lowB) ) / 2*width;
240 : // else return 0;
241 : } else {
242 0 : plumed_merror("function type does not exist");
243 : }
244 : return 0;
245 : }
246 :
247 4860 : double HistogramBead::uboundDerivative( const double& x ) const {
248 : plumed_dbg_assert(init && periodicity!=unset );
249 : double upperB;
250 4860 : if( type==gaussian ) {
251 4860 : upperB = difference( x, highb ) / ( sqrt(2.0) * width );
252 4860 : return exp( -upperB*upperB ) / ( sqrt(2*pi)*width );
253 0 : } else if ( type==triangular ) {
254 0 : plumed_error();
255 : // upperB = fabs( difference( x, highb ) / width );
256 : // if( upperB<1 ) return ( 1 - (upperB) ) / 2*width;
257 : // else return 0;
258 : } else {
259 0 : plumed_merror("function type does not exist");
260 : }
261 : return 0;
262 : }
263 :
264 2523 : }
|
