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      178360 : class DiagonalKernelParams {
      37             : public:
      38             :   std::vector<double> at;
      39             :   std::vector<double> sigma;
      40             :   double height;
      41             :   static bool bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args );
      42             :   static bool bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args );
      43             :   static bool setKernelAndCheckHeight( DiagonalKernelParams& kp, std::size_t ndim, const std::vector<double>& args );
      44             :   static std::size_t getNumberOfParameters( const DiagonalKernelParams& kp );
      45             :   static void getSigmaProjections( const DiagonalKernelParams& kp, std::vector<double>& support );
      46             :   static double evaluateR2( const RegularKernel<DiagonalKernelParams>& p, const DiagonalKernelParams& kp, View<const double> x, View<double> paramderivs );
      47             : };
      48             : 
      49             : class NonDiagonalKernelParams {
      50             : public:
      51             :   std::vector<double> at;
      52             :   Matrix<double> sigma, metric;
      53             :   double height;
      54             :   static bool bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args );
      55             :   static bool bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args );
      56             :   static bool setKernelAndCheckHeight( NonDiagonalKernelParams& kp, std::size_t ndim, const std::vector<double>& args );
      57             :   static std::size_t getNumberOfParameters( const NonDiagonalKernelParams& kp );
      58             :   static void getSigmaProjections( const NonDiagonalKernelParams& kp, std::vector<double>& support );
      59             :   static double evaluateR2( const RegularKernel<NonDiagonalKernelParams>& p, const NonDiagonalKernelParams& kp, View<const double> x, View<double> paramderivs );
      60             : };
      61             : 
      62             : class DiscreteKernel {
      63             : public:
      64             :   static void registerKeywords( Keywords& keys ) {}
      65             :   static void read( DiscreteKernel& p, ActionWithArguments* action, const std::vector<Value*>& args ) {}
      66             :   static void setArgumentDomain( const unsigned& i, DiscreteKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 ) {}
      67             :   static void getSupport( DiscreteKernel& params, const DiagonalKernelParams& kp, double dp2cutoff, std::vector<double>& support ) {}
      68             :   static double calc( const DiscreteKernel& params, const DiagonalKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs );
      69             : };
      70             : 
      71             : class HistogramBeadKernel {
      72             : public:
      73             :   std::vector<HistogramBead> beads;
      74             :   std::vector<double> gridspacing;
      75             :   static void registerKeywords( Keywords& keys );
      76             :   static void read( HistogramBeadKernel& p, ActionWithArguments* action, const std::vector<Value*>& args );
      77             :   static void setArgumentDomain( const unsigned& i, HistogramBeadKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 );
      78             :   static void getSupport( HistogramBeadKernel& params, const DiagonalKernelParams& kp, double dp2cutoff, std::vector<double>& support );
      79             :   static double calc( const HistogramBeadKernel& params, const DiagonalKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs );
      80             : };
      81             : 
      82             : template <class K>
      83             : class RegularKernel {
      84             : public:
      85             :   bool canusevol;
      86             :   SwitchingFunction switchingFunction;
      87             :   std::vector<bool> periodic;
      88             :   std::vector<double> max_minus_min, inv_max_minus_min;
      89             :   static void registerKeywords( Keywords& keys );
      90             :   static void read( RegularKernel& p, ActionWithArguments* action, const std::vector<Value*>& args );
      91             :   static void setArgumentDomain( const unsigned& i, RegularKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 );
      92             :   static double difference( const RegularKernel& params, unsigned i, const double& val1, const double& val2 );
      93             :   static void getSupport( const RegularKernel& params, const K& kp, double dp2cutoff, std::vector<double>& support );
      94             :   static double calc( const RegularKernel& params, const K& kp, View<const double> x, View<double> der, View<double> paramderivs );
      95             : };
      96             : 
      97             : template <class K>
      98         244 : void RegularKernel<K>::registerKeywords( Keywords& keys ) {
      99         244 :   keys.add("compulsory","KERNEL","GAUSSIAN","the kernel function you are using.");
     100         244 : }
     101             : 
     102             : template <class K>
     103          48 : void RegularKernel<K>::read( RegularKernel& p, ActionWithArguments* action, const std::vector<Value*>& args ) {
     104             :   std::string kerneltype;
     105          96 :   action->parse("KERNEL",kerneltype);
     106             :   std::string errors;
     107         432 :   for(auto & c: kerneltype) {
     108         384 :     c = std::toupper(c);
     109             :   }
     110          48 :   p.canusevol = (kerneltype=="GAUSSIAN");
     111          96 :   p.switchingFunction.set( kerneltype + " R_0=1.0 NOSTRETCH", errors );
     112          48 :   if( errors.length()!=0 ) {
     113           0 :     action->error("problem reading switching function description " + errors);
     114             :   }
     115          48 :   p.periodic.resize( args.size() );
     116          48 :   p.max_minus_min.resize( args.size() );
     117          48 :   p.inv_max_minus_min.resize( args.size() );
     118          48 : }
     119             : 
     120             : template <class K>
     121          74 : void RegularKernel<K>::setArgumentDomain( const unsigned& i, RegularKernel& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 ) {
     122          74 :   params.periodic[i] = isp;
     123          74 :   if( params.periodic[i] ) {
     124             :     double min, max;
     125          35 :     Tools::convert( min1, min );
     126          35 :     Tools::convert( max1, max );
     127          35 :     params.max_minus_min[i]=max-min;
     128          35 :     params.inv_max_minus_min[i]=1.0/params.max_minus_min[i];
     129             :   }
     130          74 : }
     131             : 
     132             : template <class K>
     133   113633781 : double RegularKernel<K>::difference( const RegularKernel<K>& params, unsigned i, const double& val1, const double& val2 ) {
     134   113633781 :   if( !params.periodic[i] ) {
     135    22673791 :     return val1 - val2;
     136             :   }
     137    90959990 :   return params.max_minus_min[i]*Tools::pbc( params.inv_max_minus_min[i]*( val1 - val2 ) );
     138             : }
     139             : 
     140             : template <class K>
     141          46 : void RegularKernel<K>::getSupport( const RegularKernel<K>& params, const K& kp, double dp2cutoff, std::vector<double>& support ) {
     142          46 :   K::getSigmaProjections( kp, support );
     143         120 :   for(unsigned i=0; i<support.size(); ++i) {
     144          74 :     support[i] = sqrt(2.0*dp2cutoff)*support[i];
     145             :   }
     146          46 : }
     147             : 
     148             : template <class K>
     149    38366411 : double RegularKernel<K>::calc( const RegularKernel<K>& params, const K& kp, View<const double> x, View<double> der, View<double> paramderivs ) {
     150    38366411 :   double r2 = K::evaluateR2( params, kp, x, paramderivs );
     151    38366411 :   double dval, val = kp.height*params.switchingFunction.calculateSqr( r2, dval );
     152    38366411 :   dval *= kp.height;
     153   152000192 :   for(unsigned i=0; i<der.size(); ++i) {
     154   113633781 :     der[i] += dval*paramderivs[i];
     155   113633781 :     paramderivs[i] = -dval*paramderivs[i];
     156             :   }
     157    38366411 :   paramderivs[2*kp.at.size()] = val / kp.height;
     158    38366411 :   return val;
     159             : }
     160             : 
     161      108263 : class VonMissesKernelParams {
     162             : public:
     163             :   std::vector<double> at;
     164             :   double concentration;
     165             :   double norm;
     166             :   double height;
     167             :   static bool bandwidthIsConstant( std::size_t ndim, const std::vector<Value*>& args );
     168             :   static bool bandwidthsAllSame( std::size_t ndim, const std::vector<Value*>& args );
     169             :   static bool setKernelAndCheckHeight( VonMissesKernelParams& kp, std::size_t ndim, const std::vector<double>& argval );
     170             :   static std::size_t getNumberOfParameters( const VonMissesKernelParams& kp );
     171             : };
     172             : 
     173             : class UniversalVonMisses {
     174             : public:
     175             :   std::string kerneltype;
     176             :   SwitchingFunction switchingFunction;
     177             :   static void registerKeywords( Keywords& keys ) {}
     178             :   static void read( UniversalVonMisses& p, ActionWithArguments* action, const std::vector<Value*>& args ) {}
     179             :   static void setArgumentDomain( const unsigned& i, UniversalVonMisses& params, const double& spacing, const bool isp, const std::string& min1, const std::string& max1 ) {}
     180             :   static double calc( const UniversalVonMisses& params, const VonMissesKernelParams& kp, View<const double> x, View<double> der, View<double> paramderivs );
     181             : };
     182             : 
     183             : template <class K, class P>
     184             : class SumOfKernels {
     185             : public:
     186             :   P params;
     187             :   std::vector<K> kernelParams;
     188             : /// This is used to setup the input gridobject's bounds with the grid data from values
     189             :   static void registerKeywords( Keywords& keys );
     190             :   static void read( SumOfKernels<K,P>& func, ActionWithArguments* action, const std::vector<Value*>& args, function::FunctionOptions& options );
     191             :   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 );
     192             : };
     193             : 
     194             : template <class K, class P>
     195         282 : void SumOfKernels<K,P>::registerKeywords( Keywords& keys ) {
     196         282 :   P::registerKeywords( keys );
     197         282 : }
     198             : 
     199             : template <class K, class P>
     200          66 : void SumOfKernels<K,P>::read( SumOfKernels<K,P>& func, ActionWithArguments* action, const std::vector<Value*>& args, function::FunctionOptions& options ) {
     201             : // Read the universal parameters for the kernel
     202          66 :   P::read( func.params, action, args );
     203          66 : }
     204             : 
     205             : template <class K, class P>
     206      183915 : 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 ) {
     207      183915 :   values[0] = 0;
     208      657974 :   for(unsigned i=0; i<der.size(); ++i) {
     209      474059 :     der[i] = 0;
     210             :   }
     211      183915 :   std::size_t nparams = K::getNumberOfParameters( func.kernelParams[0] );
     212    65166809 :   for(unsigned i=0; i<klist.size(); ++i) {
     213    64982894 :     values[0] += P::calc( func.params, func.kernelParams[klist[i]], args, der, View<double>( paramderivs.data() + i*nparams, nparams ) );
     214             :   }
     215      183915 : }
     216             : 
     217             : }
     218             : }
     219             : #endif