Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-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 "SumOfKernels.h"
23 :
24 : namespace PLMD {
25 : namespace gridtools {
26 :
27 32294 : double DiscreteKernel::calc( const DiscreteKernel& params, const DiagonalKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs ) {
28 32294 : return kp.height[0];
29 : }
30 :
31 40 : void HistogramBeadKernel::registerKeywords( Keywords& keys ) {
32 40 : keys.add("compulsory","KERNEL","GAUSSIAN","the kernel function you are using.");
33 40 : }
34 :
35 18 : void HistogramBeadKernel::read( HistogramBeadKernel& p, ActionWithArguments* action, const std::vector<Value*>& args ) {
36 : std::string mykerneltype;
37 18 : action->parse("KERNEL",mykerneltype);
38 18 : p.kerneltype=HistogramBead::getKernelType(mykerneltype);
39 18 : p.periodic_v.resize( args.size(), 0 );
40 18 : p.domainmin_v.resize( args.size(), 0 );
41 18 : p.domainmax_v.resize( args.size(), 0 );
42 18 : p.gridspacing_v.resize( args.size() );
43 18 : }
44 :
45 18 : void HistogramBeadKernel::setArgumentDomain( const unsigned& i, HistogramBeadKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 ) {
46 18 : params.gridspacing_v[i] = spacing;
47 18 : if( isp ) {
48 18 : params.periodic_v[i] = 1;
49 18 : Tools::convert( min1, params.domainmin_v[i] );
50 18 : Tools::convert( max1, params.domainmax_v[i] );
51 : }
52 18 : }
53 :
54 18 : void HistogramBeadKernel::getSupport( HistogramBeadKernel& params, const View<const double>& shape, std::vector<double>& support ) {
55 18 : HistogramBead mybead( params.kerneltype, 0, 1, 0.5 );
56 36 : for(unsigned i=0; i<support.size(); ++i) {
57 18 : if( params.periodic[i] ) {
58 18 : mybead.isPeriodic( params.domainmin[i], params.domainmax[i] );
59 : } else {
60 : mybead.isNotPeriodic();
61 : }
62 18 : mybead.set( 0, params.gridspacing[i], shape[i] );
63 18 : support[i] = mybead.getCutoff();
64 : }
65 18 : }
66 :
67 533600 : double HistogramBeadKernel::calc( const HistogramBeadKernel& params, const DiagonalKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs ) {
68 : // The derivatives with respect to the bandwidths are not set below: the derivative buffers
69 : // of the task managers are reused between tasks (and are not initialised on the GPU)
70 2134400 : for(unsigned i=0; i<paramderivs.size(); ++i) {
71 1600800 : paramderivs[i] = 0;
72 : }
73 533600 : double val = kp.height[0];
74 533600 : HistogramBead mybead( params.kerneltype, 0, 1, 0.5 );
75 1067200 : for(unsigned i=0; i<x.size(); ++i) {
76 533600 : if( params.periodic[i] ) {
77 533600 : mybead.isPeriodic( params.domainmin[i], params.domainmax[i] );
78 : } else {
79 : mybead.isNotPeriodic();
80 : }
81 533600 : paramderivs[i] = mybead.calculateWithCutoff( kp.at[i], x[i], x[i]+params.gridspacing[i], kp.sigma[i], der[i] );
82 533600 : val = val*paramderivs[i];
83 : }
84 1067200 : for(unsigned i=0; i<x.size(); ++i) {
85 533600 : if( fabs(paramderivs[i])>epsilon ) {
86 445921 : paramderivs[i] = der[i]*val / paramderivs[i];
87 445921 : der[i] = 0.0; // This derivative is set equal to zero because I am not sure what its proper value should be.
88 : }
89 : }
90 533600 : paramderivs[2*kp.at.size()] = val / kp.height[0];
91 533600 : return val;
92 : }
93 :
94 468494 : bool DiagonalKernelParams::setKernelAndCheckHeight( View<double>& params, std::size_t ndim, const std::vector<double>& argval ) {
95 : View<double> loc( params.data(), ndim );
96 468494 : if( argval.size()==ndim+1 ) {
97 124560 : for(unsigned i=0; i<ndim; ++i) {
98 63068 : loc[i] = argval[i];
99 : }
100 : } else {
101 407002 : plumed_assert( argval.size()==2*ndim+1 );
102 407002 : View<double> sig( params.data()+ndim, ndim );
103 942337 : for(unsigned i=0; i<ndim; ++i) {
104 535335 : loc[i] = argval[i];
105 535335 : sig[i] = argval[ndim+i];
106 : }
107 : }
108 468494 : View<double,1> h( params.data()+2*ndim );
109 468494 : h[0] = argval[argval.size()-1];
110 468494 : return fabs( argval[argval.size()-1] )>epsilon;
111 : }
112 :
113 91 : bool DiagonalKernelParams::bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args ) {
114 199 : for(unsigned i=ndim; i<args.size()-1; ++i) {
115 108 : if( !args[i]->isConstant() ) {
116 : return false;
117 : }
118 : }
119 : return true;
120 : }
121 :
122 91 : bool DiagonalKernelParams::bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args ) {
123 199 : for(unsigned i=ndim; i<args.size()-1; ++i) {
124 108 : if( !args[i]->allElementsEqual() ) {
125 : return false;
126 : }
127 : }
128 : return true;
129 : }
130 :
131 90 : std::size_t DiagonalKernelParams::getNumberOfShapeParameters( const std::size_t& nargs ) {
132 90 : return nargs;
133 : }
134 :
135 58 : void DiagonalKernelParams::getSigmaProjections( const View<const double>& shape, std::vector<double>& support ) {
136 144 : for(unsigned i=0; i<support.size(); ++i) {
137 86 : support[i] = shape[i];
138 : }
139 58 : }
140 :
141 38356031 : double DiagonalKernelParams::evaluateR2( const RegularKernel<DiagonalKernelParams>& p, const DiagonalKernelParams& kp, View<const double> x, View<double> paramderivs ) {
142 : double r2 = 0;
143 151979432 : for(unsigned i=0; i<x.size(); ++i) {
144 113623401 : double tmp = RegularKernel<DiagonalKernelParams>::difference( p, i, x[i], kp.at[i] );
145 113623401 : if( fabs(tmp)<epsilon ) {
146 81002 : paramderivs[i]=0;
147 : } else {
148 113542399 : paramderivs[i] = tmp/(kp.sigma[i]*kp.sigma[i]);
149 : }
150 113623401 : r2 += tmp*paramderivs[i];
151 : }
152 38356031 : return r2;
153 : }
154 :
155 51 : bool NonDiagonalKernelParams::setKernelAndCheckHeight( View<double>& params, std::size_t ndim, const std::vector<double>& argval ) {
156 : // centers (ndim), covariance matrix (ndim*ndim) and height
157 51 : plumed_assert( argval.size()==ndim*(ndim+1)+1 );
158 : View<double> loc( params.data(), ndim );
159 51 : Matrix<double> sigma( ndim, ndim ), metric( ndim, ndim );
160 153 : for(unsigned i=0; i<ndim; ++i) {
161 102 : loc[i] = argval[i];
162 306 : for(unsigned j=0; j<ndim; ++j) {
163 204 : sigma[i][j] = argval[ndim + i*ndim + j];
164 : }
165 : }
166 51 : Invert( sigma, metric );
167 51 : View2D<double> met( params.data() + ndim, ndim, ndim );
168 153 : for(unsigned i=0; i<ndim; ++i) {
169 306 : for(unsigned j=0; j<ndim; ++j) {
170 204 : met[i][j] = metric[i][j];
171 : }
172 : }
173 51 : View<double,1> h( params.data() + ndim*(1+ndim) );
174 51 : h[0] = argval[argval.size()-1];
175 102 : return fabs( argval[argval.size()-1] )>epsilon;
176 : }
177 :
178 1 : bool NonDiagonalKernelParams::bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args ) {
179 1 : return DiagonalKernelParams::bandwidthIsConstant( ndim, args );
180 : }
181 :
182 1 : bool NonDiagonalKernelParams::bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args ) {
183 1 : return DiagonalKernelParams::bandwidthsAllSame( ndim, args );
184 : }
185 :
186 1 : std::size_t NonDiagonalKernelParams::getNumberOfShapeParameters( const std::size_t& nargs ) {
187 1 : return nargs*nargs;
188 : }
189 :
190 1 : void NonDiagonalKernelParams::getSigmaProjections( const View<const double>& shape, std::vector<double>& support ) {
191 : // shape holds the metric M, i.e. the inverse of the covariance matrix C. The kernel is cut
192 : // where (x-a)^T M (x-a) reaches the cutoff, and the smallest box containing that ellipsoid
193 : // has half-widths proportional to sqrt(C_ii) (the same as sigma_i for a diagonal kernel).
194 : std::size_t n = support.size();
195 1 : Matrix<double> metric(n,n), covar(n,n);
196 3 : for(unsigned i=0; i<n; ++i) {
197 6 : for(unsigned j=0; j<n; ++j) {
198 4 : metric[i][j] = shape[i*n + j];
199 : }
200 : }
201 1 : Invert( metric, covar );
202 3 : for(unsigned i=0; i<n; i++) {
203 2 : support[i] = sqrt( covar[i][i] );
204 : }
205 1 : }
206 :
207 3762 : double NonDiagonalKernelParams::evaluateR2( const RegularKernel<NonDiagonalKernelParams>& p, const NonDiagonalKernelParams& kp, View<const double> x, View<double> paramderivs ) {
208 : double r2 = 0;
209 : double dp_j, dp_k;
210 11286 : for(unsigned j=0; j<x.size(); ++j) {
211 7524 : dp_j = RegularKernel<NonDiagonalKernelParams>::difference( p, j, x[j], kp.at[j] );
212 22572 : for(unsigned k=0; k<x.size(); ++k) {
213 15048 : if( j==k ) {
214 : dp_k = dp_j;
215 : } else {
216 7524 : dp_k = RegularKernel<NonDiagonalKernelParams>::difference( p, k, x[k], kp.at[k] );
217 : }
218 15048 : paramderivs[j] += kp.metric[j][k]*dp_k;
219 15048 : r2 += dp_j*dp_k*kp.metric[j][k];
220 : }
221 : }
222 3762 : return r2;
223 : }
224 :
225 18858417 : double UniversalVonMisses::calc( const UniversalVonMisses& params, const VonMissesKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs ) {
226 : // The derivatives with respect to the concentration (and to the height if it is zero) are not set below:
227 : // the derivative buffers of the task managers are reused between tasks (and are not initialised on the GPU)
228 113150502 : for(unsigned i=0; i<paramderivs.size(); ++i) {
229 94292085 : paramderivs[i] = 0;
230 : }
231 18858417 : double dot=x[0]*kp.at[0] + x[1]*kp.at[1] + x[2]*kp.at[2];
232 18858417 : double newval = kp.height[0]*exp( kp.concentration[0]*dot );
233 18858417 : der[0] += newval*kp.concentration[0]*kp.at[0];
234 18858417 : der[1] += newval*kp.concentration[0]*kp.at[1];
235 18858417 : der[2] += newval*kp.concentration[0]*kp.at[2];
236 18858417 : paramderivs[0] = newval*kp.concentration[0]*x[0];
237 18858417 : paramderivs[1] = newval*kp.concentration[0]*x[1];
238 18858417 : paramderivs[2] = newval*kp.concentration[0]*x[2];
239 18858417 : if( fabs(kp.height[0])>epsilon ) {
240 16677104 : paramderivs[4] = newval / kp.height[0];
241 : }
242 18858417 : return newval;
243 : }
244 :
245 11 : bool VonMissesKernelParams::bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args ) {
246 11 : return args[3]->isConstant();
247 : }
248 :
249 11 : bool VonMissesKernelParams::bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args ) {
250 11 : return args[3]->allElementsEqual();
251 : }
252 :
253 132191 : bool VonMissesKernelParams::setKernelAndCheckHeight( View<double>& params, std::size_t ndim, const std::vector<double>& argval ) {
254 : plumed_dbg_assert( argval.size()==5 );
255 : View<double,3> loc( params.data() );
256 132191 : loc[0] = argval[0];
257 132191 : loc[1] = argval[1];
258 132191 : loc[2] = argval[2];
259 132191 : View<double,1> conc( params.data()+3 );
260 132191 : conc[0] = argval[3];
261 132191 : View<double,1> h( params.data()+4 );
262 132191 : h[0] = argval[4];
263 132191 : return true;
264 : }
265 :
266 11 : std::size_t VonMissesKernelParams::getNumberOfShapeParameters( const std::size_t& nargs ) {
267 11 : return 1;
268 : }
269 :
270 : }
271 : }
|
