Line data    Source code

       1             : #ifndef __PLUMED_ttsketch_BasisFunc_h
       2             : #define __PLUMED_ttsketch_BasisFunc_h
       3             : #include "tools/Matrix.h"
       4             : #include <utility>
       5             : #include <vector>
       6             : 
       7             : namespace PLMD {
       8             : namespace ttsketch {
       9             : 
      10             : // Represents a 1D basis expansion for one CV dimension, used to express the
      11             : // bias potential as a tensor train over basis function indices. Two types are
      12             : // supported:
      13             : //
      14             : //   Fourier (kernel_=false): orthonormal sinusoidal functions on [a, b].
      15             : //     phi_1(x) = 1/sqrt(2L)
      16             : //     phi_{2k}(x) = 1/sqrt(L) * cos(pi*(x-shift)*k/L)   (k = 1, 2, ...)
      17             : //     phi_{2k+1}(x) = 1/sqrt(L) * sin(pi*(x-shift)*k/L) (k = 1, 2, ...)
      18             : //
      19             : //   Kernel (kernel_=true): periodized Gaussian kernels on a uniform grid,
      20             : //     plus a constant function phi_1(x) = 1.
      21             : //
      22             : // "Conv" (convolution) mode pre-smooths basis functions by convolving with a
      23             : // Gaussian of width w_. For Fourier bases this multiplies each harmonic by
      24             : // exp(-(pi*w*k)^2 / (2L^2)); for kernel bases it widens each kernel
      25             : // analytically. This yields smooth bias forces at no extra runtime cost.
      26             : class BasisFunc {
      27             : 
      28             : private:
      29             :   std::pair<double, double> dom_; // CV domain [a, b]
      30             :   int nbasis_;                    // number of basis functions n (including constant phi_1)
      31             :   double L_;                      // domain half-length: L = (b-a)/2
      32             :   double shift_;                  // domain midpoint: shift = (a+b)/2
      33             :   double inv_L_;                  // 1/L, precomputed to avoid repeated division in Fourier evaluations
      34             :   double sqrt_inv_L_;             // 1/sqrt(L), precomputed to avoid repeated pow() in Fourier evaluations
      35             :   double inv_sqrt_2L_;            // 1/sqrt(2L), the value of phi_1(x) for the Fourier basis
      36             :   double w_;                      // smoothing width for conv mode (Fourier: sigma in x-units; kernel: Gaussian hill width)
      37             :   bool kernel_;                   // true = kernel basis, false = Fourier basis
      38             :   double dx_;                     // kernel basis function width (std dev of each periodized Gaussian)
      39             :   std::vector<double> centers_;   // centers of kernel basis functions (size nbasis-1, uniformly spaced)
      40             :   Matrix<double> ginv_;           // pseudo-inverse of the Gram matrix (kernel basis only)
      41             :   double gramCondition_;          // ratio between the largest and the smallest eigenvalue of the Gram matrix
      42             :   int gramDropped_;               // eigen-directions of the Gram matrix left out of ginv_ (see the constructor)
      43             :   Matrix<double> gram_;           // Gram matrix: gram_(i,j) = int phi_i(x)*phi_j(x) dx (kernel basis only)
      44             :   std::vector<double> damp_;      // conv mode, Fourier basis: damp_[k] = exp(-(pi*w*k/L)^2/2), k = 0..nbasis/2
      45             : 
      46             : public:
      47             :   BasisFunc();
      48             :   BasisFunc(std::pair<double, double> dom, int nbasis, double w, bool kernel, double dx);
      49             :   // Evaluate Fourier basis function phi_pos(x) (1-indexed, not convolved).
      50             :   double fourier(double x, int pos) const;
      51             :   // Evaluate kernel basis function phi_pos(x) (1-indexed, not convolved).
      52             :   // Uses periodic images (k=-1,0,1) so the function wraps smoothly across domain boundaries.
      53             :   double gaussian(double x, int pos) const;
      54             :   // Derivative of fourier(x, pos) w.r.t. x.
      55             :   double fourierd(double x, int pos) const;
      56             :   // Derivative of gaussian(x, pos) w.r.t. x.
      57             :   double gaussiand(double x, int pos) const;
      58             :   // Evaluate basis function phi_pos(x). If conv=true, applies Gaussian smoothing of width w_.
      59             :   double operator()(double x, int pos, bool conv) const;
      60             :   // Derivative of operator()(x, pos, conv) w.r.t. x.
      61             :   double grad(double x, int pos, bool conv) const;
      62             :   // Values and derivatives of all the basis functions at x:
      63             :   // phi[pos-1] = operator()(x, pos, conv) and dphi[pos-1] = grad(x, pos, conv), pos = 1..nbasis.
      64             :   // Equivalent to nbasis calls of operator() and grad(), but the Fourier harmonics are
      65             :   // obtained from a single sin/cos pair by the angle-addition recurrence.
      66             :   void evalAll(double x, bool conv, double* phi, double* dphi) const;
      67             :   int nbasis() const {
      68          92 :     return this->nbasis_;
      69             :   }
      70             :   const std::pair<double, double>& dom() const {
      71             :     return this->dom_;
      72             :   }
      73             :   double w() const {
      74             :     return this->w_;
      75             :   }
      76             :   // Integral of phi_pos over the domain: int_a^b phi_pos(x) dx.
      77             :   // Used to compute the partition function Z in covMat().
      78             :   double int0(int pos) const;
      79             :   // Integral of x*phi_pos over the domain: int_a^b x*phi_pos(x) dx.
      80             :   // Used to compute the marginal mean E[x] in covMat().
      81             :   double int1(int pos) const;
      82             :   // Integral of x^2*phi_pos over the domain: int_a^b x^2*phi_pos(x) dx.
      83             :   // Used to compute E[x^2] for the marginal variance in covMat().
      84             :   double int2(int pos) const;
      85             :   bool kernel() const {
      86        3596 :     return this->kernel_;
      87             :   }
      88             :   const Matrix<double>& ginv() const {
      89             :     return this->ginv_;
      90             :   }
      91             :   const Matrix<double>& gram() const {
      92             :     return this->gram_;
      93             :   }
      94             :   // Eigenvalues of the Gram matrix smaller than gramTolerance() times the largest one are
      95             :   // left out of ginv(): their directions are combinations of the basis functions that
      96             :   // vanish (to working precision) and whose dual coefficients would be meaningless.
      97             :   static constexpr double gramTolerance() {
      98             :     return 1e-10;
      99             :   }
     100             :   double gramCondition() const {
     101           3 :     return this->gramCondition_;
     102             :   }
     103             :   int gramDropped() const {
     104           3 :     return this->gramDropped_;
     105             :   }
     106             :   double center(int pos) const {
     107         189 :     return this->kernel_ ? this->centers_[pos - 1] : 0.0;
     108             :   }
     109             :   double dx() const {
     110         192 :     return this->dx_;
     111             :   }
     112             : };
     113             : 
     114             : }
     115             : }
     116             : 
     117             : #endif