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 : #ifndef __PLUMED_gridtools_SumOfKernels_h
23 : #define __PLUMED_gridtools_SumOfKernels_h
24 :
25 : #include "function/FunctionSetup.h"
26 : #include "tools/SwitchingFunction.h"
27 : #include "tools/HistogramBead.h"
28 : #include "tools/Matrix.h"
29 :
30 : namespace PLMD {
31 : namespace gridtools {
32 :
33 : template <class K>
34 : class RegularKernel;
35 :
36 : class DiagonalKernelParams {
37 : public:
38 : View<const double> at;
39 : View<const double> sigma;
40 : View<const double,1> height;
41 : DiagonalKernelParams( const View<const double>& params, const std::size_t ndim ):
42 : at(params.data(),ndim),
43 38921925 : sigma(params.data()+ndim,ndim),
44 38921925 : height(params.data()+2*ndim) {
45 : }
46 : static bool bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args );
47 : static bool bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args );
48 : static bool setKernelAndCheckHeight( View<double>& params, std::size_t ndim, const std::vector<double>& args );
49 : static std::size_t getNumberOfParameters( const std::size_t& nargs );
50 : static std::size_t getNumberOfShapeParameters( const std::size_t& nargs );
51 : static void getSigmaProjections( const View<const double>& shape, std::vector<double>& support );
52 : #pragma acc routine seq
53 : static double evaluateR2( const RegularKernel<DiagonalKernelParams>& p, const DiagonalKernelParams& kp, View<const double> x, View<double> paramderivs );
54 : };
55 :
56 : inline
57 : std::size_t DiagonalKernelParams::getNumberOfParameters( const std::size_t& nargs ) {
58 7568 : return 2*nargs + 1;
59 : }
60 :
61 : class NonDiagonalKernelParams {
62 : public:
63 : View<const double> at;
64 : View2D<const double> metric;
65 : View<const double,1> height;
66 3762 : NonDiagonalKernelParams( const View<const double>& params, const std::size_t ndim ):
67 : at(params.data(),ndim),
68 3762 : metric(params.data()+ndim,ndim,ndim),
69 3762 : height(params.data()+ndim*(1+ndim)) {
70 3762 : }
71 : static bool bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args );
72 : static bool bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args );
73 : static bool setKernelAndCheckHeight( View<double>& params, std::size_t ndim, const std::vector<double>& args );
74 : static std::size_t getNumberOfParameters( const std::size_t& nargs );
75 : static std::size_t getNumberOfShapeParameters( const std::size_t& nargs );
76 : static void getSigmaProjections( const View<const double>& shape, std::vector<double>& support );
77 : #pragma acc routine seq
78 : static double evaluateR2( const RegularKernel<NonDiagonalKernelParams>& p, const NonDiagonalKernelParams& kp, View<const double> x, View<double> paramderivs );
79 : };
80 :
81 : inline
82 : std::size_t NonDiagonalKernelParams::getNumberOfParameters( const std::size_t& nargs ) {
83 1321 : return nargs + nargs*nargs + 1;
84 : }
85 :
86 : class DiscreteKernel {
87 : public:
88 : static void registerKeywords( Keywords& keys ) {}
89 : static void read( DiscreteKernel& p, ActionWithArguments* action, const std::vector<Value*>& args ) {}
90 : static void setArgumentDomain( const unsigned& i, DiscreteKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 ) {}
91 : static void getSupport( DiscreteKernel& params, const View<const double>& shape, std::vector<double>& support ) {}
92 : #pragma acc routine seq
93 : static double calc( const DiscreteKernel& params, const DiagonalKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs );
94 : void update() {}
95 : void toACCDevice() const {}
96 : void removeFromACCDevice() const {}
97 : };
98 :
99 : class HistogramBeadKernel {
100 : public:
101 : HistogramBead::KernelType kerneltype;
102 : std::vector<std::size_t> periodic_v;
103 : std::size_t* periodic{nullptr};
104 : std::vector<double> domainmin_v;
105 : double* domainmin{nullptr};
106 : std::vector<double> domainmax_v;
107 : double* domainmax{nullptr};
108 : std::vector<double> gridspacing_v;
109 : double* gridspacing{nullptr};
110 : static void registerKeywords( Keywords& keys );
111 : static void read( HistogramBeadKernel& p, ActionWithArguments* action, const std::vector<Value*>& args );
112 : static void setArgumentDomain( const unsigned& i, HistogramBeadKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 );
113 : static void getSupport( HistogramBeadKernel& params, const View<const double>& shape, std::vector<double>& support );
114 : #pragma acc routine seq
115 : static double calc( const HistogramBeadKernel& params, const DiagonalKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs );
116 : void update() {
117 18 : periodic = periodic_v.data();
118 18 : domainmin = domainmin_v.data();
119 18 : domainmax = domainmax_v.data();
120 18 : gridspacing = gridspacing_v.data();
121 : }
122 : void toACCDevice() const {
123 : #pragma acc enter data copyin(this[0:1],kerneltype,periodic[0:periodic_v.size()],domainmin[0:domainmin_v.size()],domainmax[0:domainmax_v.size()],gridspacing[0:gridspacing_v.size()])
124 : }
125 : void removeFromACCDevice() const {
126 : #pragma acc exit data delete(gridspacing[0:gridspacing_v.size()],domainmax[0:domainmax_v.size()],domainmin[0:domainmin_v.size()],periodic[0:periodic_v.size()],kerneltype,this[0:1])
127 : }
128 : };
129 :
130 : template <class K>
131 : class RegularKernel {
132 : public:
133 : bool canusevol;
134 : double cutoff;
135 : #ifdef __PLUMED_HAS_OPENACC
136 : SwitchingFunctionAccelerable switchingFunction;
137 : #else
138 : SwitchingFunction switchingFunction;
139 : #endif //__PLUMED_HAS_OPENACC
140 : std::vector<int> periodic_v;
141 : int* periodic{nullptr};
142 : std::vector<double> max_minus_min_v;
143 : double* max_minus_min{nullptr};
144 : std::vector<double> inv_max_minus_min_v;
145 : double* inv_max_minus_min{nullptr};
146 : static void registerKeywords( Keywords& keys );
147 : static void read( RegularKernel& p, ActionWithArguments* action, const std::vector<Value*>& args );
148 : static void setArgumentDomain( const unsigned& i, RegularKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 );
149 : static double difference( const RegularKernel& params, unsigned i, const double& val1, const double& val2 );
150 : static void getSupport( const RegularKernel& params, const View<const double>& shape, std::vector<double>& support );
151 : static double calc( const RegularKernel& params, const K& kp, View<const double> x, View<double> der, View<double> paramderivs );
152 : void update() {
153 59 : periodic = periodic_v.data();
154 59 : max_minus_min = max_minus_min_v.data();
155 59 : inv_max_minus_min = inv_max_minus_min_v.data();
156 : }
157 : #ifdef __PLUMED_HAS_OPENACC
158 : void toACCDevice() const {
159 : #pragma acc enter data copyin(this[0:1],periodic[0:periodic_v.size()],max_minus_min[0:max_minus_min_v.size()],inv_max_minus_min[0:inv_max_minus_min_v.size()])
160 : switchingFunction.toACCDevice();
161 : }
162 : void removeFromACCDevice() const {
163 : switchingFunction.removeFromACCDevice();
164 : #pragma acc exit data delete(inv_max_minus_min[0:inv_max_minus_min_v.size()],max_minus_min[0:max_minus_min_v.size()],periodic[0:periodic_v.size()],this[0:1])
165 : }
166 : #endif //__PLUMED_HAS_OPENACC
167 : };
168 :
169 : template <class K>
170 322 : void RegularKernel<K>::registerKeywords( Keywords& keys ) {
171 : std::string strcut;
172 322 : Tools::convert( sqrt(2*dp2cutoff), strcut );
173 322 : keys.add("compulsory","KERNEL","GAUSSIAN D_MAX=" + strcut,"the kernel function you are using.");
174 322 : }
175 :
176 : template <class K>
177 66 : void RegularKernel<K>::read( RegularKernel& p, ActionWithArguments* action, const std::vector<Value*>& args ) {
178 : std::string kerneltype;
179 132 : action->parse("KERNEL",kerneltype);
180 : std::string errors;
181 815 : for(auto & c: kerneltype) {
182 749 : c = std::toupper(c);
183 : }
184 66 : p.canusevol = (kerneltype.find("GAUSSIAN")!=std::string::npos );
185 66 : p.switchingFunction.set( kerneltype + " R_0=1.0", errors );
186 66 : if( p.switchingFunction.get_dmax()<std::numeric_limits<double>::max() ) {
187 15 : p.cutoff = p.switchingFunction.get_dmax();
188 30 : action->log<<" using switching function for kernel with r0="<<p.switchingFunction.description()<<" and cutoff "<<p.cutoff<<"\n";
189 : } else {
190 51 : p.cutoff = sqrt(2*dp2cutoff);
191 102 : action->warning("using a switching function without a D_MAX parameter for the kernel is not recommended. As discussed in the manual, "
192 : "when D_MAX is unset the cutoff is set equal to sqrt(2*6.25). However, density from any given kernel "
193 : "will be accumulated at points that are more than this D_MAX from the center of the kernel. The reason why "
194 : "is explained on the manual pages for KDE in the section on selecting kernels.");
195 : }
196 66 : if( errors.length()!=0 ) {
197 0 : action->error("problem reading switching function description " + errors);
198 : }
199 66 : p.periodic_v.resize( args.size(), 0 );
200 66 : p.max_minus_min_v.resize( args.size() );
201 66 : p.inv_max_minus_min_v.resize( args.size() );
202 66 : }
203 :
204 : template <class K>
205 88 : void RegularKernel<K>::setArgumentDomain( const unsigned& i, RegularKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 ) {
206 88 : if( isp ) {
207 : double min, max;
208 36 : Tools::convert( min1, min );
209 36 : Tools::convert( max1, max );
210 36 : params.periodic_v[i] = 1;
211 36 : params.max_minus_min_v[i]=max-min;
212 36 : params.inv_max_minus_min_v[i]=1.0/params.max_minus_min_v[i];
213 : }
214 88 : }
215 :
216 : template <class K>
217 113638449 : double RegularKernel<K>::difference( const RegularKernel<K>& params, unsigned i, const double& val1, const double& val2 ) {
218 113638449 : if( params.periodic[i]==0 ) {
219 22670935 : return val1 - val2;
220 : }
221 90967514 : return params.max_minus_min[i]*Tools::pbc( params.inv_max_minus_min[i]*( val1 - val2 ) );
222 : }
223 :
224 : template <class K>
225 59 : void RegularKernel<K>::getSupport( const RegularKernel<K>& params, const View<const double>& shape, std::vector<double>& support ) {
226 59 : K::getSigmaProjections( shape, support );
227 147 : for(unsigned i=0; i<support.size(); ++i) {
228 88 : support[i] = params.cutoff*support[i];
229 : }
230 59 : }
231 :
232 : template <class K>
233 38359793 : double RegularKernel<K>::calc( const RegularKernel<K>& params, const K& kp, View<const double> x, View<double> der, View<double> paramderivs ) {
234 : // Not all the parameter derivatives are set below (e.g. those with respect to the bandwidths)
235 : // and evaluateR2 for non-diagonal kernels accumulates into them, so they must start from zero:
236 : // the derivative buffers of the task managers are reused between tasks (and are not initialised on the GPU)
237 303988960 : for(unsigned i=0; i<paramderivs.size(); ++i) {
238 265629167 : paramderivs[i] = 0;
239 : }
240 38359793 : double r2 = K::evaluateR2( params, kp, x, paramderivs );
241 38359793 : double dval, val = kp.height[0]*params.switchingFunction.calculateSqr( r2, dval );
242 38359793 : dval *= kp.height[0];
243 151990718 : for(unsigned i=0; i<der.size(); ++i) {
244 113630925 : der[i] += dval*paramderivs[i];
245 113630925 : paramderivs[i] = -dval*paramderivs[i];
246 : }
247 : // The height is always the last parameter of the kernel (DiagonalKernelParams: 2*ndim, NonDiagonalKernelParams: ndim*(1+ndim))
248 38359793 : paramderivs[paramderivs.size()-1] = val / kp.height[0];
249 38359793 : return val;
250 : }
251 :
252 : class VonMissesKernelParams {
253 : public:
254 : View<const double,3> at;
255 : View<const double,1> concentration;
256 : View<const double,1> height;
257 : VonMissesKernelParams( const View<const double>& params, const std::size_t ndim ):
258 : at(params.data()),
259 18858417 : concentration(params.data()+3),
260 18858417 : height(params.data()+4) {
261 : }
262 : static bool bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args );
263 : static bool bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args );
264 : static bool setKernelAndCheckHeight( View<double>& params, std::size_t ndim, const std::vector<double>& argval );
265 : static std::size_t getNumberOfParameters( const std::size_t& nargs );
266 : static std::size_t getNumberOfShapeParameters( const std::size_t& nargs );
267 : };
268 :
269 : inline
270 : std::size_t VonMissesKernelParams::getNumberOfParameters( const std::size_t& nargs ) {
271 : return 5;
272 : }
273 :
274 : class UniversalVonMisses {
275 : public:
276 : static void registerKeywords( Keywords& keys ) {}
277 : static void read( UniversalVonMisses& p, ActionWithArguments* action, const std::vector<Value*>& args ) {}
278 : static void setArgumentDomain( const unsigned& i, UniversalVonMisses& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 ) {}
279 : #pragma acc routine seq
280 : static double calc( const UniversalVonMisses& params, const VonMissesKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs );
281 : void update() {}
282 : void toACCDevice() const {}
283 : void removeFromACCDevice() const {}
284 : };
285 :
286 : template <class K, class P>
287 13 : class SumOfKernels {
288 : public:
289 : P params;
290 : std::vector<double> kernelParams_v;
291 : double* kernelParams{nullptr};
292 : /// This is used to setup the input gridobject's bounds with the grid data from values
293 : static void registerKeywords( Keywords& keys );
294 : static void read( SumOfKernels<K,P>& func, ActionWithArguments* action, const std::vector<Value*>& args, function::FunctionOptions& options );
295 : static void calc( View<const std::size_t> klist, const SumOfKernels<K,P>& func, View<const double> args, View<double> values, View<double> der, View<double> paramderivs );
296 : void update() {
297 2599 : kernelParams=kernelParams_v.data();
298 : }
299 : void toACCDevice() const {
300 : #pragma acc enter data copyin(this[0:1],kernelParams[0:kernelParams_v.size()])
301 : params.toACCDevice();
302 : }
303 : void removeFromACCDevice() const {
304 : params.removeFromACCDevice();
305 : #pragma acc exit data delete(kernelParams[0:kernelParams_v.size()],this[0:1])
306 : }
307 : };
308 :
309 : template <class K, class P>
310 362 : void SumOfKernels<K,P>::registerKeywords( Keywords& keys ) {
311 362 : P::registerKeywords( keys );
312 362 : }
313 :
314 : template <class K, class P>
315 84 : void SumOfKernels<K,P>::read( SumOfKernels<K,P>& func, ActionWithArguments* action, const std::vector<Value*>& args, function::FunctionOptions& options ) {
316 : // Read the universal parameters for the kernel
317 84 : P::read( func.params, action, args );
318 84 : }
319 :
320 : template <class K, class P>
321 190619 : void SumOfKernels<K,P>::calc( View<const std::size_t> klist, const SumOfKernels<K,P>& func, View<const double> args, View<double> values, View<double> der, View<double> paramderivs ) {
322 190619 : values[0] = 0;
323 675784 : for(unsigned i=0; i<der.size(); ++i) {
324 485165 : der[i] = 0;
325 : }
326 : std::size_t ndim = args.size();
327 : std::size_t nparams = K::getNumberOfParameters( ndim );
328 57974723 : for(unsigned i=0; i<klist.size(); ++i) {
329 57784104 : values[0] += P::calc( func.params, K(View<const double>(func.kernelParams + nparams*klist[i],nparams),ndim), args, der, View<double>( paramderivs.data() + i*nparams, nparams ) );
330 : }
331 190619 : }
332 :
333 : }
334 : }
335 : #endif
|
