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 "KernelFunctions.h"
23 : #include <limits>
24 : #include "IFile.h"
25 : #include <cmath>
26 :
27 : namespace PLMD {
28 :
29 : //+PLUMEDOC INTERNAL kernelfunctions
30 : /*
31 : Functions that are used to construct histograms
32 :
33 : Constructing histograms is something you learned to do relatively early in life. You perform an experiment a number of times,
34 : count the number of times each result comes up and then draw a bar graph that describes how often each of the results came up.
35 : This only works when there are a finite number of possible results. If the result a number between 0 and 1 the bar chart is
36 : less easy to draw as there are as many possible results as there are numbers between zero and one - an infinite number.
37 : To resolve this problem we replace probability, \f$P\f$ with probability density, \f$\pi\f$, and write the probability of getting
38 : a number between \f$a\f$ and \f$b\f$ as:
39 :
40 : \f[
41 : P = \int_{a}^b \textrm{d}x \pi(x)
42 : \f]
43 :
44 : To calculate probability densities from a set of results we use a process called kernel density estimation.
45 : Histograms are accumulated by adding up kernel functions, \f$K\f$, with finite spatial extent, that integrate to one.
46 : These functions are centered on each of the \f$n\f$-dimensional data points, \f$\mathbf{x}_i\f$. The overall effect of this
47 : is that each result we obtain in our experiments contributes to the probability density in a finite sized region of the space.
48 :
49 : Expressing all this mathematically in kernel density estimation we write the probability density as:
50 :
51 : \f[
52 : \pi(\mathbf{x}) = \sum_i K\left[ (\mathbf{x} - \mathbf{x}_i)^T \Sigma (\mathbf{x} - \mathbf{x}_i) \right]
53 : \f]
54 :
55 : where \f$\Sigma\f$ is an \f$n \times n\f$ matrix called the bandwidth that controls the spatial extent of
56 : the kernel. Whenever we accumulate a histogram (e.g. in \ref HISTOGRAM or in \ref METAD) we use this
57 : technique.
58 :
59 : There is thus some flexibility in the particular function we use for \f$K[\mathbf{r}]\f$ in the above.
60 : The following variants are available.
61 :
62 : <table align=center frame=void width=95%% cellpadding=5%%>
63 : <tr>
64 : <td> TYPE </td> <td> FUNCTION </td>
65 : </tr> <tr>
66 : <td> gaussian </td> <td> \f$f(r) = \frac{1}{(2 \pi)^{n} \sqrt{|\Sigma^{-1}|}} \exp\left(-0.5 r^2 \right)\f$ </td>
67 : </tr> <tr>
68 : <td> truncated-gaussian </td> <td> \f$f(r) = \frac{1}{(2 \pi)^{n} \sqrt{|\Sigma^{-1}|} \left(\frac{\mathrm{erf}(-6.25/sqrt{2}) - \mathrm{erf}(-6.25/sqrt{2})}{2}\right)^n} \exp\left(-0.5 r^2 \right)\f$ </td>
69 : </tr> <tr>
70 : <td> triangular </td> <td> \f$f(r) = \frac{3}{V} ( 1 - | r | )H(1-|r|) \f$ </td>
71 : </tr> <tr>
72 : <td> uniform </td> <td> \f$f(r) = \frac{1}{V}H(1-|r|)\f$ </td>
73 : </tr>
74 : </table>
75 :
76 : In the above \f$H(y)\f$ is a function that is equal to one when \f$y>0\f$ and zero when \f$y \le 0\f$. \f$n\f$ is
77 : the dimensionality of the vector \f$\mathbf{x}\f$ and \f$V\f$ is the volume of an ellipse in an \f$n\f$ dimensional
78 : space which is given by:
79 :
80 : \f{eqnarray*}{
81 : V &=& | \Sigma^{-1} | \frac{ \pi^{\frac{n}{2}} }{\left( \frac{n}{2} \right)! } \qquad \textrm{for even} \quad n \\
82 : V &=& | \Sigma^{-1} | \frac{ 2^{\frac{n+1}{2}} \pi^{\frac{n-1}{2}} }{ n!! }
83 : \f}
84 :
85 : In \ref METAD the normalization constants are ignored so that the value of the function at \f$r=0\f$ is equal
86 : to one. In addition in \ref METAD we must be able to differentiate the bias in order to get forces. This limits
87 : the kernels we can use in this method. Notice also that Gaussian kernels should have infinite support. When used
88 : with grids, however, they are assumed to only be non-zero over a finite range. The difference between the
89 : truncated-gaussian and regular gaussian is that the truncated gaussian is scaled so that its integral over the grid
90 : is equal to one when it is normalized. The integral of a regular gaussian when it is evaluated on a grid will be
91 : slightly less that one because of the truncation of a function that should have infinite support.
92 : */
93 : //+ENDPLUMEDOC
94 :
95 0 : KernelFunctions::KernelFunctions( const std::string& input ) {
96 :
97 0 : if(!dp2cutoffNoStretch()) {
98 0 : stretchA=dp2cutoffA;
99 0 : stretchB=dp2cutoffB;
100 : }
101 :
102 0 : std::vector<std::string> data=Tools::getWords(input);
103 0 : std::string name=data[0];
104 : data.erase(data.begin());
105 :
106 0 : std::vector<double> at;
107 0 : bool foundc = Tools::parseVector(data,"CENTER",at);
108 0 : if(!foundc) {
109 0 : plumed_merror("failed to find center keyword in definition of kernel");
110 : }
111 0 : std::vector<double> sig;
112 0 : bool founds = Tools::parseVector(data,"SIGMA",sig);
113 0 : if(!founds) {
114 0 : plumed_merror("failed to find sigma keyword in definition of kernel");
115 : }
116 :
117 0 : bool multi=false;
118 0 : Tools::parseFlag(data,"MULTIVARIATE",multi);
119 0 : bool vonmises=false;
120 0 : Tools::parseFlag(data,"VON-MISSES",vonmises);
121 0 : if( center.size()==1 && multi ) {
122 0 : plumed_merror("one dimensional kernel cannot be multivariate");
123 : }
124 0 : if( center.size()==1 && vonmises ) {
125 0 : plumed_merror("one dimensional kernel cannot be von-misses");
126 : }
127 0 : if( center.size()==1 && sig.size()!=1 ) {
128 0 : plumed_merror("size mismatch between center size and sigma size");
129 : }
130 0 : if( multi && center.size()>1 && sig.size()!=0.5*center.size()*(center.size()-1) ) {
131 0 : plumed_merror("size mismatch between center size and sigma size");
132 : }
133 0 : if( !multi && center.size()>1 && sig.size()!=center.size() ) {
134 0 : plumed_merror("size mismatch between center size and sigma size");
135 : }
136 :
137 : double h;
138 0 : bool foundh = Tools::parse(data,"HEIGHT",h);
139 0 : if( !foundh) {
140 0 : h=1.0;
141 : }
142 :
143 0 : if( multi ) {
144 0 : setData( at, sig, name, "MULTIVARIATE", h );
145 0 : } else if( vonmises ) {
146 0 : setData( at, sig, name, "VON-MISSES", h );
147 : } else {
148 0 : setData( at, sig, name, "DIAGONAL", h );
149 : }
150 0 : }
151 :
152 6143 : KernelFunctions::KernelFunctions( const std::vector<double>& at,
153 : const std::vector<double>& sig,
154 : const std::string& type,
155 : const std::string& mtype,
156 6143 : const double& w ) {
157 :
158 6143 : if(!dp2cutoffNoStretch()) {
159 6143 : stretchA=dp2cutoffA;
160 6143 : stretchB=dp2cutoffB;
161 : }
162 :
163 6143 : setData( at, sig, type, mtype, w );
164 6143 : }
165 :
166 0 : KernelFunctions::KernelFunctions( const KernelFunctions* in ):
167 0 : dtype(in->dtype),
168 0 : ktype(in->ktype),
169 0 : center(in->center),
170 0 : width(in->width),
171 0 : height(in->height) {
172 0 : if(!dp2cutoffNoStretch()) {
173 0 : stretchA=dp2cutoffA;
174 0 : stretchB=dp2cutoffB;
175 : }
176 0 : }
177 :
178 6143 : void KernelFunctions::setData( const std::vector<double>& at,
179 : const std::vector<double>& sig,
180 : const std::string& type,
181 : const std::string& mtype,
182 : const double& w ) {
183 :
184 6143 : height=w;
185 6143 : center.resize( at.size() );
186 18208 : for(unsigned i=0; i<at.size(); ++i) {
187 12065 : center[i]=at[i];
188 : }
189 6143 : width.resize( sig.size() );
190 22576 : for(unsigned i=0; i<sig.size(); ++i) {
191 16433 : width[i]=sig[i];
192 : }
193 6143 : if( mtype=="MULTIVARIATE" ) {
194 4368 : dtype=matrixType::multi;
195 1775 : } else if( mtype=="VON-MISSES" ) {
196 0 : dtype=matrixType::vonmises;
197 1775 : } else if( mtype=="DIAGONAL" ) {
198 1775 : dtype=matrixType::diagonal;
199 : } else {
200 0 : plumed_merror(mtype + " is not a valid metric type");
201 : }
202 :
203 : // Setup the kernel type
204 6143 : if(type=="GAUSSIAN" || type=="gaussian" ) {
205 0 : ktype=kernelType::gaussian;
206 6143 : } else if(type=="STRETCHED-GAUSSIAN" || type=="stretched-gaussian" ) {
207 6143 : ktype=kernelType::stretchedgaussian;
208 0 : } else if(type=="TRUNCATED-GAUSSIAN" || type=="truncated-gaussian" ) {
209 0 : ktype=kernelType::truncatedgaussian;
210 0 : } else if(type=="UNIFORM" || type=="uniform") {
211 0 : ktype=kernelType::uniform;
212 0 : } else if(type=="TRIANGULAR" || type=="triangular") {
213 0 : ktype=kernelType::triangular;
214 : } else {
215 0 : plumed_merror(type+" is an invalid kernel type\n");
216 : }
217 6143 : }
218 :
219 0 : void KernelFunctions::normalize( const std::vector<Value*>& myvals ) {
220 :
221 : double det=1.;
222 : unsigned ncv=ndim();
223 0 : if(dtype==matrixType::diagonal) {
224 0 : for(unsigned i=0; i<width.size(); ++i) {
225 0 : det*=width[i]*width[i];
226 : }
227 0 : } else if(dtype==matrixType::multi) {
228 : // det is the determinant of the inverse of mymatrix
229 0 : Matrix<double> mymatrix( getMatrix() );
230 : double logd;
231 0 : logdet( mymatrix, logd );
232 0 : det=std::exp(-logd);
233 : }
234 0 : if( dtype==matrixType::diagonal || dtype==matrixType::multi ) {
235 : double volume;
236 0 : if( ktype==kernelType::gaussian || ktype==kernelType::stretchedgaussian ) {
237 0 : volume=pow( 2*pi, 0.5*ncv ) * pow( det, 0.5 );
238 0 : } else if( ktype==kernelType::truncatedgaussian ) {
239 : // This makes it so the gaussian integrates to one over the range over which it has support
240 : const double DP2CUTOFF=std::sqrt(6.25);
241 0 : volume=pow( 2*pi, 0.5*ncv ) * pow( det, 0.5 ) * pow( 0.5 * ( erf(DP2CUTOFF) - erf(-DP2CUTOFF) ), ncv);
242 0 : } else if( ktype==kernelType::uniform || ktype==kernelType::triangular ) {
243 0 : if( ncv%2==1 ) {
244 : double dfact=1;
245 0 : for(unsigned i=1; i<ncv; i+=2) {
246 0 : dfact*=static_cast<double>(i);
247 : }
248 0 : volume=( pow( pi, (ncv-1)/2 ) ) * ( pow( 2., (ncv+1)/2 ) ) / dfact;
249 : } else {
250 : double fact=1.;
251 0 : for(unsigned i=1; i<ncv/2; ++i) {
252 0 : fact*=static_cast<double>(i);
253 : }
254 0 : volume=pow( pi,ncv/2 ) / fact;
255 : }
256 0 : if(ktype==kernelType::uniform) {
257 0 : volume*=det;
258 0 : } else if(ktype==kernelType::triangular) {
259 0 : volume*=det / 3.;
260 : }
261 : } else {
262 0 : plumed_merror("not a valid kernel type");
263 : }
264 0 : height /= volume;
265 0 : return;
266 : }
267 0 : plumed_assert( dtype==matrixType::vonmises && ktype==kernelType::gaussian );
268 : // Now calculate determinant for aperiodic variables
269 : unsigned naper=0;
270 0 : for(unsigned i=0; i<ncv; ++i) {
271 0 : if( !myvals[i]->isPeriodic() ) {
272 0 : naper++;
273 : }
274 : }
275 : // Now construct sub matrix
276 : double volume=1;
277 0 : if( naper>0 ) {
278 : unsigned isub=0;
279 0 : Matrix<double> mymatrix( getMatrix() ), mysub( naper, naper );
280 0 : for(unsigned i=0; i<ncv; ++i) {
281 0 : if( myvals[i]->isPeriodic() ) {
282 : continue;
283 : }
284 : unsigned jsub=0;
285 0 : for(unsigned j=0; j<ncv; ++j) {
286 0 : if( myvals[j]->isPeriodic() ) {
287 0 : continue;
288 : }
289 0 : mysub( isub, jsub ) = mymatrix( i, j );
290 0 : jsub++;
291 : }
292 0 : isub++;
293 : }
294 0 : Matrix<double> myisub( naper, naper );
295 : double logd;
296 0 : Invert( mysub, myisub );
297 0 : logdet( myisub, logd );
298 0 : det=std::exp(logd);
299 0 : volume=pow( 2*pi, 0.5*ncv ) * pow( det, 0.5 );
300 : }
301 :
302 : // Calculate volume of periodic variables
303 : unsigned nper=0;
304 0 : for(unsigned i=0; i<ncv; ++i) {
305 0 : if( myvals[i]->isPeriodic() ) {
306 0 : nper++;
307 : }
308 : }
309 :
310 : // Now construct sub matrix
311 0 : if( nper>0 ) {
312 : unsigned isub=0;
313 0 : Matrix<double> mymatrix( getMatrix() ), mysub( nper, nper );
314 0 : for(unsigned i=0; i<ncv; ++i) {
315 0 : if( !myvals[i]->isPeriodic() ) {
316 : continue;
317 : }
318 : unsigned jsub=0;
319 0 : for(unsigned j=0; j<ncv; ++j) {
320 0 : if( !myvals[j]->isPeriodic() ) {
321 0 : continue;
322 : }
323 0 : mysub( isub, jsub ) = mymatrix( i, j );
324 0 : jsub++;
325 : }
326 0 : isub++;
327 : }
328 0 : Matrix<double> eigvec( nper, nper );
329 0 : std::vector<double> eigval( nper );
330 0 : diagMat( mysub, eigval, eigvec );
331 : unsigned iper=0;
332 : volume=1;
333 0 : for(unsigned i=0; i<ncv; ++i) {
334 0 : if( myvals[i]->isPeriodic() ) {
335 0 : volume *= myvals[i]->getMaxMinusMin()*Tools::bessel0(eigval[iper])*std::exp(-eigval[iper]);
336 0 : iper++;
337 : }
338 : }
339 : }
340 0 : height /= volume;
341 : }
342 :
343 10985 : double KernelFunctions::getCutoff(const double baseWidth) const {
344 : const double DP2CUTOFF=6.25;
345 10985 : if( ktype==kernelType::gaussian
346 : || ktype==kernelType::truncatedgaussian
347 10985 : || ktype==kernelType::stretchedgaussian ) {
348 10985 : return std::sqrt(2.0*DP2CUTOFF)*baseWidth;
349 0 : } else if(ktype==kernelType::triangular ) {
350 : return baseWidth;
351 0 : } else if(ktype==kernelType::uniform) {
352 : return baseWidth;
353 : } else {
354 0 : plumed_merror("No valid kernel type");
355 : }
356 : return 0.0;
357 : }
358 :
359 5609 : std::vector<double> KernelFunctions::getContinuousSupport( ) const {
360 : unsigned ncv=ndim();
361 5609 : std::vector<double> support( ncv );
362 5609 : if(dtype==matrixType::diagonal) {
363 3490 : for(unsigned i=0; i<ncv; ++i) {
364 2249 : support[i]=getCutoff(width[i]);
365 : }
366 4368 : } else if(dtype==matrixType::multi) {
367 4368 : Matrix<double> mymatrix( getMatrix() ), myinv( ncv,ncv );
368 4368 : Invert(mymatrix,myinv);
369 4368 : Matrix<double> myautovec(ncv,ncv);
370 4368 : std::vector<double> myautoval(ncv);
371 4368 : diagMat(myinv,myautoval,myautovec);
372 : double maxautoval=0.;
373 : unsigned ind_maxautoval=0;
374 13104 : for (unsigned i=0; i<ncv; i++) {
375 8736 : if(myautoval[i]>maxautoval) {
376 : maxautoval=myautoval[i];
377 : ind_maxautoval=i;
378 : }
379 : }
380 13104 : for(unsigned i=0; i<ncv; ++i) {
381 8736 : double extent=std::fabs(std::sqrt(maxautoval)*myautovec(i,ind_maxautoval));
382 8736 : support[i]=getCutoff( extent );
383 : }
384 : } else {
385 0 : plumed_merror("cannot find support if metric is not multi or diagonal type");
386 : }
387 5609 : return support;
388 : }
389 :
390 3787 : std::vector<unsigned> KernelFunctions::getSupport( const std::vector<double>& dx ) const {
391 3787 : plumed_assert( ndim()==dx.size() );
392 3787 : std::vector<unsigned> support( dx.size() );
393 3787 : std::vector<double> vv=getContinuousSupport( );
394 11128 : for(unsigned i=0; i<dx.size(); ++i) {
395 7341 : support[i]=static_cast<unsigned>(ceil( vv[i]/dx[i] ));
396 : }
397 3787 : return support;
398 : }
399 :
400 970504 : double KernelFunctions::evaluate( const std::vector<Value*>& pos,
401 : std::vector<double>& derivatives,
402 : const bool usederiv,
403 : const bool doInt,
404 : const double lowI_,
405 : const double uppI_) const {
406 : plumed_dbg_assert( pos.size()==ndim() && derivatives.size()==ndim() );
407 : #ifndef NDEBUG
408 : if( usederiv ) {
409 : plumed_massert( ktype!=kernelType::uniform, "step function can not be differentiated" );
410 : }
411 : #endif
412 970504 : if(doInt) {
413 : plumed_dbg_assert(center.size()==1);
414 0 : if(pos[0]->get()<lowI_) {
415 0 : pos[0]->set(lowI_);
416 : }
417 0 : if(pos[0]->get()>uppI_) {
418 0 : pos[0]->set(uppI_);
419 : }
420 : }
421 : double r2=0;
422 970504 : if(dtype==matrixType::diagonal) {
423 0 : for(unsigned i=0; i<ndim(); ++i) {
424 0 : derivatives[i]=-pos[i]->difference( center[i] ) / width[i];
425 0 : r2+=derivatives[i]*derivatives[i];
426 0 : derivatives[i] /= width[i];
427 : }
428 970504 : } else if(dtype==matrixType::multi) {
429 970504 : Matrix<double> mymatrix( getMatrix() );
430 2911512 : for(unsigned i=0; i<mymatrix.nrows(); ++i) {
431 : double dp_i, dp_j;
432 1941008 : derivatives[i]=0;
433 1941008 : dp_i=-pos[i]->difference( center[i] );
434 5823024 : for(unsigned j=0; j<mymatrix.ncols(); ++j) {
435 3882016 : if(i==j) {
436 : dp_j=dp_i;
437 : } else {
438 1941008 : dp_j=-pos[j]->difference( center[j] );
439 : }
440 :
441 3882016 : derivatives[i]+=mymatrix(i,j)*dp_j;
442 3882016 : r2+=dp_i*dp_j*mymatrix(i,j);
443 : }
444 : }
445 0 : } else if(dtype==matrixType::vonmises) {
446 0 : std::vector<double> costmp( ndim() ), sintmp( ndim() ), sinout( ndim(), 0.0 );
447 0 : for(unsigned i=0; i<ndim(); ++i) {
448 0 : if( pos[i]->isPeriodic() ) {
449 0 : sintmp[i]=std::sin( 2.*pi*(pos[i]->get() - center[i])/pos[i]->getMaxMinusMin() );
450 0 : costmp[i]=std::cos( 2.*pi*(pos[i]->get() - center[i])/pos[i]->getMaxMinusMin() );
451 : } else {
452 0 : sintmp[i]=pos[i]->get() - center[i];
453 0 : costmp[i]=1.0;
454 : }
455 : }
456 :
457 0 : Matrix<double> mymatrix( getMatrix() );
458 0 : for(unsigned i=0; i<mymatrix.nrows(); ++i) {
459 0 : derivatives[i]=0;
460 0 : if( pos[i]->isPeriodic() ) {
461 0 : r2+=2*( 1 - costmp[i] )*mymatrix(i,i);
462 : } else {
463 0 : r2+=sintmp[i]*sintmp[i]*mymatrix(i,i);
464 : }
465 0 : for(unsigned j=0; j<mymatrix.ncols(); ++j) {
466 0 : if( i!=j ) {
467 0 : sinout[i]+=mymatrix(i,j)*sintmp[j];
468 : }
469 : }
470 0 : derivatives[i] = mymatrix(i,i)*sintmp[i] + sinout[i]*costmp[i];
471 0 : if( pos[i]->isPeriodic() ) {
472 0 : derivatives[i] *= (2*pi/pos[i]->getMaxMinusMin());
473 : }
474 : }
475 0 : for(unsigned i=0; i<sinout.size(); ++i) {
476 0 : r2+=sintmp[i]*sinout[i];
477 : }
478 : }
479 : double kderiv, kval;
480 970504 : if(ktype==kernelType::gaussian || ktype==kernelType::truncatedgaussian) {
481 0 : kval=height*std::exp(-0.5*r2);
482 0 : kderiv=-kval;
483 970504 : } else if(ktype==kernelType::stretchedgaussian) {
484 970504 : auto dp=0.5*r2;
485 970504 : if(dp<dp2cutoff) {
486 604935 : auto ee=std::exp(-0.5*r2);
487 604935 : kval=height*(stretchA*ee+stretchB);
488 604935 : kderiv=-height*stretchA*ee;
489 : } else {
490 : kval=0.0;
491 : kderiv=0.0;
492 : }
493 : } else {
494 0 : double r=std::sqrt(r2);
495 0 : if(ktype==kernelType::triangular) {
496 0 : if( r<1.0 ) {
497 : if(r==0) {
498 : kderiv=0;
499 : }
500 : kderiv=-1;
501 0 : kval=height*( 1. - std::fabs(r) );
502 : } else {
503 : kval=0.;
504 : kderiv=0.;
505 : }
506 0 : } else if(ktype==kernelType::uniform) {
507 : kderiv=0.;
508 0 : if(r<1.0) {
509 0 : kval=height;
510 : } else {
511 : kval=0;
512 : }
513 : } else {
514 0 : plumed_merror("Not a valid kernel type");
515 : }
516 0 : kderiv*=height / r ;
517 : }
518 2911512 : for(unsigned i=0; i<ndim(); ++i) {
519 1941008 : derivatives[i]*=kderiv;
520 : }
521 970504 : if(doInt) {
522 0 : if((pos[0]->get() <= lowI_ || pos[0]->get() >= uppI_) && usederiv )
523 0 : for(unsigned i=0; i<ndim(); ++i) {
524 0 : derivatives[i]=0;
525 : }
526 : }
527 970504 : return kval;
528 : }
529 :
530 3787 : bool KernelFunctions::getSeparableGaussian(double & h, double & a, double & b, double & r2cut, std::vector<double> & w) const {
531 3787 : if(dtype!=matrixType::diagonal) {
532 : return false;
533 : }
534 1057 : if(ktype==kernelType::gaussian || ktype==kernelType::truncatedgaussian) {
535 0 : a=1.0;
536 0 : b=0.0;
537 0 : r2cut=std::numeric_limits<double>::max();
538 1057 : } else if(ktype==kernelType::stretchedgaussian) {
539 1057 : a=stretchA;
540 1057 : b=stretchB;
541 1057 : r2cut=2.0*dp2cutoff;
542 : } else {
543 : return false;
544 : }
545 1057 : h=height;
546 1057 : w=width;
547 1057 : return true;
548 : }
549 :
550 5051 : std::unique_ptr<KernelFunctions> KernelFunctions::read( IFile* ifile,
551 : const bool& cholesky,
552 : const std::vector<std::string>& valnames ) {
553 : double h;
554 10102 : if( !ifile->scanField("height",h) ) {
555 : return NULL;
556 : };
557 :
558 : std::string sss;
559 5051 : ifile->scanField("multivariate",sss);
560 5051 : std::string ktype="stretched-gaussian";
561 10102 : if( ifile->FieldExist("kerneltype") ) {
562 7918 : ifile->scanField("kerneltype",ktype);
563 : }
564 5051 : plumed_massert( sss=="false"
565 : || sss=="true"
566 : || sss=="von-misses",
567 : "multivariate flag must be either false, true or von-misses");
568 :
569 : // Read the position of the center
570 5051 : std::vector<double> cc( valnames.size() );
571 14932 : for(unsigned i=0; i<valnames.size(); ++i) {
572 9881 : ifile->scanField(valnames[i],cc[i]);
573 : }
574 :
575 5051 : std::vector<double> sig;
576 5051 : if( sss=="false" ) {
577 683 : sig.resize( valnames.size() );
578 1828 : for(unsigned i=0; i<valnames.size(); ++i) {
579 1145 : ifile->scanField("sigma_"+valnames[i],sig[i]);
580 1145 : if( !cholesky ) {
581 0 : sig[i]=std::sqrt(sig[i]);
582 : }
583 : }
584 : return Tools::make_unique<KernelFunctions>(cc, sig, ktype, "DIAGONAL", h);
585 : }
586 :
587 4368 : unsigned ncv=valnames.size();
588 4368 : sig.resize( (ncv*(ncv+1))/2 );
589 4368 : Matrix<double> upper(ncv,ncv);
590 4368 : Matrix<double> lower(ncv,ncv);
591 4368 : Matrix<double> mymult(ncv,ncv);
592 4368 : Matrix<double> invmatrix(ncv,ncv);
593 13104 : for(unsigned i=0; i<ncv; ++i) {
594 21840 : for(unsigned j=0; j<ncv-i; j++) {
595 39312 : ifile->scanField("sigma_" +valnames[j+i] + "_" + valnames[j], lower(j+i,j) );
596 13104 : upper(j,j+i)=lower(j+i,j);
597 13104 : mymult(j+i,j)=mymult(j,j+i)=lower(j+i,j);
598 : }
599 : }
600 4368 : if( cholesky ) {
601 4368 : mult(lower,upper,mymult);
602 : }
603 4368 : Invert( mymult, invmatrix );
604 : unsigned k=0;
605 13104 : for(unsigned i=0; i<ncv; i++) {
606 21840 : for(unsigned j=i; j<ncv; j++) {
607 13104 : sig[k]=invmatrix(i,j);
608 13104 : k++;
609 : }
610 : }
611 4368 : if( sss=="true" ) {
612 : return Tools::make_unique<KernelFunctions>(cc, sig, ktype, "MULTIVARIATE", h);
613 : }
614 : return Tools::make_unique<KernelFunctions>( cc, sig, ktype, "VON-MISSES", h );
615 : }
616 :
617 : }
|
