Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2021 The VES code team
       3             :    (see the PEOPLE-VES file at the root of this folder for a list of names)
       4             : 
       5             :    See http://www.ves-code.org for more information.
       6             : 
       7             :    This file is part of VES code module.
       8             : 
       9             :    The VES code module 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             :    The VES code module 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 the VES code module.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : 
      23             : #include "BasisFunctions.h"
      24             : 
      25             : #include "core/ActionRegister.h"
      26             : #include "core/PlumedMain.h"
      27             : 
      28             : 
      29             : namespace PLMD {
      30             : namespace ves {
      31             : 
      32             : //+PLUMEDOC VES_BASISF BF_CUBIC_B_SPLINES
      33             : /*
      34             : Cubic B spline basis functions.
      35             : 
      36             : \attention
      37             : __These basis functions do not form orthogonal bases. We recommend using wavelets (\ref BF_WAVELETS) instead that do for orthogonal bases__.
      38             : 
      39             : A basis using cubic B spline functions according to \cite habermann_multidimensional_2007. See \cite ValssonPampel_Wavelets_2022 for full details.
      40             : 
      41             : The mathematical expression of the individual splines is given by
      42             : \f{align*}{
      43             :   h\left(x\right) =
      44             :   \begin{cases}
      45             :     \left(2 - \lvert x \rvert\right)^3, & 1 \leq \lvert x \rvert \leq 2\\
      46             :     4 - 6\lvert x \rvert^2 + 3 \lvert x \rvert^3,\qquad & \lvert x \rvert \leq 1\\
      47             :     0, & \text{elsewhere}.
      48             :   \end{cases}
      49             : \f}
      50             : 
      51             : The full basis consists of equidistant splines at positions \f$\mu_i\f$ which are optimized in their height:
      52             : \f{align*}{
      53             :   f_i\left(x\right) = h\left(\frac{x-\mu_i}{\sigma}\right)
      54             : \f}
      55             : 
      56             : Note that the distance between individual splines cannot be chosen freely but is equal to the width: \f$\mu_{i+1} = \mu_{i} + \sigma\f$.
      57             : 
      58             : 
      59             : The ORDER keyword of the basis set determines the number of equally sized sub-intervalls to be used.
      60             : On the borders of each of these sub-intervalls the mean \f$\mu_i\f$ of a spline function is placed.
      61             : 
      62             : The total number of basis functions is \f$\text{ORDER}+4\f$ as the constant \f$f_{0}(x)=1\f$, as well as the two splines with means just outside the interval are also included.
      63             : 
      64             : As an example two adjacent basis functions can be seen below.
      65             : The full basis consists of shifted splines in the full specified interval.
      66             : 
      67             : \image html ves_basisf-splines.png
      68             : 
      69             : When the splines are used for a periodic CV (with the PERIODIC keyword),
      70             : the sub-intervals are chosen in the same way, but only \f$\text{ORDER}+1\f$ functions
      71             : are required to fill it (the ones at the boundary coincide and the ones outside
      72             : can be omitted).
      73             : 
      74             : To avoid 'blind' optimization of the basis functions outside the currently sampled area, it is often beneficial to use the OPTIMIZATION_THRESHOLD keyword of the VES_LINEAR_EXPANSION (set it to a small value, e.g. 1e-6)
      75             : 
      76             : \par Examples
      77             : The bias is expanded with cubic B splines in the intervall from 0.0 to 10.0 specifying an order of 20.
      78             : This results in 24 basis functions.
      79             : 
      80             : \plumedfile
      81             : bf: BF_CUBIC_B_SPLINES MINIMUM=0.0 MAXIMUM=10.0 ORDER=20
      82             : \endplumedfile
      83             : 
      84             : */
      85             : //+ENDPLUMEDOC
      86             : 
      87             : class BF_CubicBsplines : public BasisFunctions {
      88             :   double spacing_;
      89             :   double inv_spacing_;
      90             :   double inv_normfactor_;
      91             :   double spline(const double, double&) const;
      92             : public:
      93             :   static void registerKeywords( Keywords&);
      94             :   explicit BF_CubicBsplines(const ActionOptions&);
      95             :   void getAllValues(const double, double&, bool&, std::vector<double>&, std::vector<double>&) const;
      96             : };
      97             : 
      98             : 
      99             : PLUMED_REGISTER_ACTION(BF_CubicBsplines,"BF_CUBIC_B_SPLINES")
     100             : 
     101             : // See DOI 10.1007/s10614-007-9092-4 for more information;
     102             : 
     103             : 
     104           5 : void BF_CubicBsplines::registerKeywords(Keywords& keys) {
     105           5 :   BasisFunctions::registerKeywords(keys);
     106          10 :   keys.add("optional","NORMALIZATION","the normalization factor that is used to normalize the basis functions by dividing the values. By default it is 2.");
     107          10 :   keys.addFlag("PERIODIC", false, "Use periodic version of basis set.");
     108           5 :   keys.remove("NUMERICAL_INTEGRALS");
     109           5 : }
     110             : 
     111           3 : BF_CubicBsplines::BF_CubicBsplines(const ActionOptions&ao):
     112           3 :   PLUMED_VES_BASISFUNCTIONS_INIT(ao)
     113             : {
     114           3 :   log.printf("  Cubic B spline basis functions, see and cite ");
     115           6 :   log << plumed.cite("Pampel and Valsson, J. Chem. Theory Comput. 18, 4127-4141 (2022) - DOI:10.1021/acs.jctc.2c00197");
     116             : 
     117           3 :   setIntrinsicInterval(intervalMin(),intervalMax());
     118           3 :   spacing_=(intervalMax()-intervalMin())/static_cast<double>(getOrder());
     119           3 :   inv_spacing_ = 1.0/spacing_;
     120             : 
     121           3 :   bool periodic = false;
     122           3 :   parseFlag("PERIODIC",periodic);
     123           4 :   if (periodic) {addKeywordToList("PERIODIC",periodic);}
     124             : 
     125             :   // 1 constant, getOrder() on interval, 1 (left) + 2 (right) at boundaries if not periodic
     126           3 :   unsigned int num_BFs = periodic ? getOrder()+1U : getOrder()+4U;
     127           3 :   setNumberOfBasisFunctions(num_BFs);
     128             : 
     129           3 :   double normfactor_=2.0;
     130           3 :   parse("NORMALIZATION",normfactor_);
     131           3 :   if(normfactor_!=2.0) {addKeywordToList("NORMALIZATION",normfactor_);}
     132           3 :   inv_normfactor_=1.0/normfactor_;
     133             : 
     134           3 :   periodic ? setPeriodic() : setNonPeriodic();
     135             :   setIntervalBounded();
     136           3 :   setType("splines_2nd-order");
     137           3 :   setDescription("Cubic B-splines (2nd order splines)");
     138           3 :   setLabelPrefix("S");
     139           3 :   setupBF();
     140           3 :   log.printf("   normalization factor: %f\n",normfactor_);
     141           3 :   checkRead();
     142           3 : }
     143             : 
     144             : 
     145        9236 : void BF_CubicBsplines::getAllValues(const double arg, double& argT, bool& inside_range, std::vector<double>& values, std::vector<double>& derivs) const {
     146             :   // plumed_assert(values.size()==numberOfBasisFunctions());
     147             :   // plumed_assert(derivs.size()==numberOfBasisFunctions());
     148        9236 :   inside_range=true;
     149        9236 :   argT=checkIfArgumentInsideInterval(arg,inside_range);
     150             :   //
     151        9236 :   values[0]=1.0;
     152        9236 :   derivs[0]=0.0;
     153             :   //
     154      212043 :   for(unsigned int i=1; i < getNumberOfBasisFunctions(); i++) {
     155      202807 :     double argx = ((argT-intervalMin())*inv_spacing_) - (static_cast<double>(i) - 2.0);
     156      202807 :     if(arePeriodic()) { // periodic range of argx is [-intervalRange/spacing,+intervalRange/spacing]
     157       64140 :       argx /= intervalRange()*inv_spacing_;
     158       64140 :       argx = Tools::pbc(argx);
     159       64140 :       argx *= (intervalRange()*inv_spacing_);
     160             :     }
     161      202807 :     values[i] = spline(argx, derivs[i]);
     162      202807 :     derivs[i] *= inv_spacing_;
     163             :   }
     164       11156 :   if(!inside_range) {for(unsigned int i=0; i<derivs.size(); i++) {derivs[i]=0.0;}}
     165        9236 : }
     166             : 
     167             : 
     168      202807 : double BF_CubicBsplines::spline(const double arg, double& deriv) const {
     169             :   double value=0.0;
     170             :   double x=arg;
     171             :   // derivative of abs(x);
     172             :   double dx = 1.0;
     173      202807 :   if(x < 0) {
     174      101208 :     x=-x;
     175             :     dx = -1.0;
     176             :   }
     177             :   //
     178      202807 :   if(x > 2) {
     179             :     value=0.0;
     180      165556 :     deriv=0.0;
     181             :   }
     182       37251 :   else if(x >= 1) {
     183       18896 :     value = ((2.0-x)*(2.0-x)*(2.0-x));
     184       18896 :     deriv = dx*(-3.0*(2.0-x)*(2.0-x));
     185             :     // value=((2.0-x)*(2.0-x)*(2.0-x))/6.0;
     186             :     // deriv=-x*x*(2.0-x)*(2.0-x);
     187             :   }
     188             :   else {
     189       18355 :     value = 4.0-6.0*x*x+3.0*x*x*x;
     190       18355 :     deriv = dx*(-12.0*x+9.0*x*x);
     191             :     // value=x*x*x*0.5-x*x+2.0/3.0;
     192             :     // deriv=(3.0/2.0)*x*x-2.0*x;
     193             :   }
     194      202807 :   value *= inv_normfactor_;
     195      202807 :   deriv *= inv_normfactor_;
     196      202807 :   return value;
     197             : }
     198             : 
     199             : 
     200             : }
     201             : }