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 : {
106 260467 : }
107 :
108 53 : std::string HistogramBead::description() const {
109 53 : std::ostringstream ostr;
110 53 : ostr<<"between "<<lowb<<" and "<<highb<<" width of gaussian window equals "<<width;
111 53 : return ostr.str();
112 53 : }
113 :
114 0 : void HistogramBead::generateBins( const std::string& params, std::vector<std::string>& bins ) {
115 0 : std::vector<std::string> data=Tools::getWords(params);
116 0 : plumed_massert(data.size()>=1,"There is no input for this keyword");
117 :
118 0 : std::string name=data[0];
119 :
120 0 : unsigned nbins; std::vector<double> range(2); std::string smear;
121 0 : bool found_nb=Tools::parse(data,"NBINS",nbins);
122 0 : plumed_massert(found_nb,"Number of bins in histogram not found");
123 0 : bool found_r=Tools::parse(data,"LOWER",range[0]);
124 0 : plumed_massert(found_r,"Lower bound for histogram not specified");
125 0 : found_r=Tools::parse(data,"UPPER",range[1]);
126 0 : plumed_massert(found_r,"Upper bound for histogram not specified");
127 0 : plumed_massert(range[0]<range[1],"Range specification is dubious");
128 0 : bool found_b=Tools::parse(data,"SMEAR",smear);
129 0 : if(!found_b) { Tools::convert(0.5,smear); }
130 :
131 0 : std::string lb,ub; double delr = ( range[1]-range[0] ) / static_cast<double>( nbins );
132 0 : for(unsigned i=0; i<nbins; ++i) {
133 0 : Tools::convert( range[0]+i*delr, lb );
134 0 : Tools::convert( range[0]+(i+1)*delr, ub );
135 0 : bins.push_back( name + " " + "LOWER=" + lb + " " + "UPPER=" + ub + " " + "SMEAR=" + smear );
136 : }
137 0 : plumed_assert(bins.size()==nbins);
138 0 : }
139 :
140 53 : void HistogramBead::set( const std::string& params, std::string& errormsg ) {
141 53 : std::vector<std::string> data=Tools::getWords(params);
142 53 : if(data.size()<1) {
143 : errormsg="No input has been specified";
144 : return;
145 : }
146 :
147 53 : std::string name=data[0]; const double DP2CUTOFF=6.25;
148 53 : if(name=="GAUSSIAN") { type=gaussian; cutoff=std::sqrt(2.0*DP2CUTOFF); }
149 0 : else if(name=="TRIANGULAR") { type=triangular; cutoff=1.; }
150 0 : else plumed_merror("cannot understand kernel type " + name );
151 :
152 : double smear;
153 53 : bool found_r=Tools::parse(data,"LOWER",lowb);
154 53 : if( !found_r ) errormsg="Lower bound has not been specified use LOWER";
155 53 : found_r=Tools::parse(data,"UPPER",highb);
156 53 : if( !found_r ) errormsg="Upper bound has not been specified use UPPER";
157 53 : if( lowb>=highb ) errormsg="Lower bound is higher than upper bound";
158 :
159 53 : smear=0.5; Tools::parse(data,"SMEAR",smear);
160 53 : width=smear*(highb-lowb); init=true;
161 53 : }
162 :
163 2668474 : void HistogramBead::set( double l, double h, double w) {
164 2668474 : init=true; lowb=l; highb=h; width=w;
165 : const double DP2CUTOFF=6.25;
166 2668474 : if( type==gaussian ) cutoff=std::sqrt(2.0*DP2CUTOFF);
167 1935352 : else if( type==triangular ) cutoff=1.;
168 0 : else plumed_error();
169 2668474 : }
170 :
171 260148 : void HistogramBead::setKernelType( const std::string& ktype ) {
172 260148 : if(ktype=="gaussian") type=gaussian;
173 26252 : else if(ktype=="triangular") type=triangular;
174 0 : else plumed_merror("cannot understand kernel type " + ktype );
175 260148 : }
176 :
177 251184 : double HistogramBead::calculate( double x, double& df ) const {
178 : plumed_dbg_assert(init && periodicity!=unset );
179 : double lowB, upperB, f;
180 251184 : if( type==gaussian ) {
181 251184 : lowB = difference( x, lowb ) / ( std::sqrt(2.0) * width );
182 251184 : upperB = difference( x, highb ) / ( std::sqrt(2.0) * width );
183 251184 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( std::sqrt(2*pi)*width );
184 251184 : f = 0.5*( erf( upperB ) - erf( lowB ) );
185 0 : } else if( type==triangular ) {
186 0 : lowB = ( difference( x, lowb ) / width );
187 0 : upperB = ( difference( x, highb ) / width );
188 0 : df=0;
189 0 : if( std::fabs(lowB)<1. ) df = (1 - std::fabs(lowB)) / width;
190 0 : if( std::fabs(upperB)<1. ) df -= (1 - std::fabs(upperB)) / width;
191 0 : if (upperB<=-1. || lowB >=1.) {
192 : f=0.;
193 : } else {
194 : double ia, ib;
195 0 : if( lowB>-1.0 ) { ia=lowB; } else { ia=-1.0; }
196 0 : if( upperB<1.0 ) { ib=upperB; } else { ib=1.0; }
197 0 : f = (ib*(2.-std::fabs(ib))-ia*(2.-std::fabs(ia)))*0.5;
198 : }
199 : } else {
200 0 : plumed_merror("function type does not exist");
201 : }
202 251184 : return f;
203 : }
204 :
205 920596 : double HistogramBead::calculateWithCutoff( double x, double& df ) const {
206 : plumed_dbg_assert(init && periodicity!=unset );
207 :
208 : double lowB, upperB, f;
209 920596 : lowB = difference( x, lowb ) / width ; upperB = difference( x, highb ) / width;
210 920596 : if( upperB<=-cutoff || lowB>=cutoff ) { df=0; return 0; }
211 :
212 832917 : if( type==gaussian ) {
213 446395 : lowB /= std::sqrt(2.0); upperB /= std::sqrt(2.0);
214 446395 : df = ( exp( -lowB*lowB ) - exp( -upperB*upperB ) ) / ( std::sqrt(2*pi)*width );
215 446395 : f = 0.5*( erf( upperB ) - erf( lowB ) );
216 386522 : } else if( type==triangular ) {
217 386522 : df=0;
218 386522 : if( std::fabs(lowB)<1. ) df = (1 - std::fabs(lowB)) / width;
219 386522 : if( std::fabs(upperB)<1. ) df -= (1 - std::fabs(upperB)) / width;
220 386522 : if (upperB<=-1. || lowB >=1.) {
221 : f=0.;
222 : } else {
223 : double ia, ib;
224 386522 : if( lowB>-1.0 ) { ia=lowB; } else { ia=-1.0; }
225 386522 : if( upperB<1.0 ) { ib=upperB; } else { ib=1.0; }
226 386522 : f = (ib*(2.-std::fabs(ib))-ia*(2.-std::fabs(ia)))*0.5;
227 : }
228 : } else {
229 0 : plumed_merror("function type does not exist");
230 : }
231 : return f;
232 : }
233 :
234 4680 : double HistogramBead::lboundDerivative( const double& x ) const {
235 4680 : if( type==gaussian ) {
236 4680 : double lowB = difference( x, lowb ) / ( std::sqrt(2.0) * width );
237 4680 : return exp( -lowB*lowB ) / ( std::sqrt(2*pi)*width );
238 0 : } else if ( type==triangular ) {
239 0 : plumed_error();
240 : // lowB = fabs( difference( x, lowb ) / width );
241 : // if( lowB<1 ) return ( 1 - (lowB) ) / 2*width;
242 : // else return 0;
243 : } else {
244 0 : plumed_merror("function type does not exist");
245 : }
246 : return 0;
247 : }
248 :
249 4680 : double HistogramBead::uboundDerivative( const double& x ) const {
250 : plumed_dbg_assert(init && periodicity!=unset );
251 4680 : if( type==gaussian ) {
252 4680 : double upperB = difference( x, highb ) / ( std::sqrt(2.0) * width );
253 4680 : return exp( -upperB*upperB ) / ( std::sqrt(2*pi)*width );
254 0 : } else if ( type==triangular ) {
255 0 : plumed_error();
256 : // upperB = fabs( difference( x, highb ) / width );
257 : // if( upperB<1 ) return ( 1 - (upperB) ) / 2*width;
258 : // else return 0;
259 : } else {
260 0 : plumed_merror("function type does not exist");
261 : }
262 : return 0;
263 : }
264 :
265 : }
|
