Use cubic spline interpolation to speed up volmap (#869)

* SPAM shouldn't be hidden.

* Skip peaks where the bin calculation would potentially overflow

* Avoid using zero bandwidth for KDE histogram when only 1 value available

* Start adding dataset that can be used to hold vectors and an associated scalar value.

* Add Sync()

* Update dependencies

* Add VECTOR_SCALAR set

* Add DataIO class for writing peaks.

* Enable peaks file output

* Use vector scalar data to store peaks, write to data file.

* Make the vector scalar set 0 dimension so it is forced to use the DataIO_Peaks format. May want to do this in a less hacky way in the future.

* Re-organize branches in volmap loop to avoid some unnecessary calculation. Slightly more efficient (1-2%) and avoids using a 'continue'.

* A line with no elements cannot be a Tinker atom line

* Start adding code to read the peaks file.

* Remove some debug info

* Add peaks read.

* Rewrite peaks read to make it more strict. Must conform to what ID_DataFormat expects.

* Use DataSet/DataFile framework for peaks in SPAM

* Ensure that the GetNextString for output file is the last possible thing parsed from the argument list

* Fix class const variable name

* Make it so that you can calculate peaks without having to write them to a file.

* Fix code comment

* Fix up dependencies

* Add class for approximating functions with splines

* Ensure table is cleared when fill is called.

* Use the splines. 36% speed increase!

* Disable some debug checking

* Woops, bring back return statement

* Remove some debug info. Add scaling factor argument.

* Add splinedx keyword

* Add code note about function calls in openmp parallel regions

* Fix up help.

* Warn if extra arguments remain after 'purewater' specified for SPAM

* Fix up SPAM documentation.

* Update the volmap manual entry.

* More manual updates for volmap. Fix up help.

* Fix up some dependencies

* Add /= operator

* Add IncrementElement() function

* Add SetGrid function

* Change SetGrid to accept index instead

* Enable DIV and MULT for grids

* Add option to read grids in as double precision

* Add function to return format type

* Add == operator for TextFormat

* Fix the operator

* Ensure debug level is set for files being written

* When writing opendx files, if the format is still the default for grids, use general format (matches previous behavior). However, if another format has been set use that instead.

* Add check for out of range value in Yval(). Save input X values for use in a more accurate function, Yval_accurate().

* Add another way of filling the table; can give better results with Yvals_accurate()

* DRR - Omit out of bounds check, rely on args being reasonable. Results in minor speedup (e.g. 5.4 fps to 6.4 fps).

* Put the range check behind a define.

* Allow volmap to be compiled for either single or double precision grid

* Add ability to compile and use exponential function.

* Add define for using the more accurate but slower table lookup

* Add some debug

* Add some debug info and a new function that uses the internal X table.

* Add spline function table check for Yval_xtable

* Add VOLMAP_USEXTABLE

* Ensure the xvalues used in the mesh are generated in the same way as during the Y eval step. Round-off makes the addition method prone to errors for small table spacings

* When enabled, make sure the out of range check happens before dx calculation so it actually takes effect

* Add splinescale arg, improve output

* Add back the ReadHelp function.

* Remove scaling argument from spline table setup. Instead, have volmap add 1.0 to the max expected value to avoid Runges phenomenon.

* Add fast exponential from N. N. Schraudolph, Neural Computation 11, 853–862 (1999)

* Add define for using fast exp from Schraudolph

* Add 64 bit version and header protection

* Fix the always inline attribute; add volmap define for using 64 bit fast exps

* Add library from https://github.com/simonpf/fastexp for testing

* Add ability to test the fast exp code from https://github.com/simonpf/fastexp

* Remove the "fast" exponential code. Was only for testing. Leave the hooks in, just in case we ever need to test it again. Fast exponential code will be kept in a separate repository.

* Remove obsolete splinescale argument

* Clean up output

* Remove old code, improve code documentation.

* Put spline info print into separate functions.

* Add debug info. Do not set up spline table if using system exp()

* Add a warning if any grid spacings are smaller than 0.4

* Use SplineFxnTable for erfc

* Make erfc things private

* Make function private

* Rename function to avoid confusion

* Make function private. Add some code docs

* Make more vars private

* Remove obsolete code

* Make more variables private

* Minor version bump for using splines to speed up volmap. In principle the defaults chosen should make it so that the results are exactly the same as with system exp(); however there is the potential for some differences, so bump the minor version instead of revision.

* Add 'splinedx' to volmap help

Co-authored-by: Daniel R. Roe <daniel.roe@nih.gov>

Author: drroe

doc/cpptraj.lyx

@@ -35435,7 +35435,7 @@ volmap [out <filename>] <mask> [radscale <factor>] [stepfac <fac>]
 \end_layout
 
 \begin_layout LyX-Code
-   	[sphere] [radii {vdw | element}]
+   	[sphere] [radii {vdw | element}] [splinedx <spacing>]
 \end_layout
 
 \begin_layout LyX-Code
@@ -35528,6 +35528,14 @@ radii
  radii.
 \end_layout
 
+\begin_layout Description
+splinedx
+\begin_inset space ~
+\end_inset
+
+<spacing> Spacing to use for cubic spline interpolation (default 0.01 Ang.).
+\end_layout
+
 \begin_layout Description
 calcpeaks If specified, peaks in the grid density will be calculated and
  saved to set <setname> with aspect 
@@ -35737,6 +35745,25 @@ buffer
 .
 \end_layout
 
+\begin_layout Standard
+The calculation is sped up by using cubic splines to interpolate the exponential
+ function when calculating the Gaussians.
+ The details are in Roe & Brooks, 
+\begin_inset Quotes eld
+\end_inset
+
+Improving the Speed of Volumetric Density Map Generation via Cubic Spline
+ Interpolation
+\begin_inset Quotes erd
+\end_inset
+
+, J.
+ Mol.
+ Graph.
+ Model.
+ (2021).
+\end_layout
+
 \begin_layout Subsection
 volume
 \end_layout

src/Action_Volmap.cpp

@@ -15,6 +15,11 @@
 #ifdef _OPENMP
 # include <omp.h>
 #endif
+//#if defined(VOLMAP_USEFASTEXPS) || defined(VOLMAP_USEFASTEXP64)
+//# incl ude "FastExp_Schraudolph.h"
+//#elif defined(VOLMAP_USEFASTEXPIEEE)
+//# incl ude "FastExp_fastexp.h"
+//#endif
 
 const double Action_Volmap::sqrt_8_pi_cubed_ = sqrt(8.0*Constants::PI*Constants::PI*Constants::PI);
 
@@ -39,7 +44,8 @@ Action_Volmap::Action_Volmap() :
   peakcut_(0.05),
   buffer_(3.0),
   radscale_(1.0),
-  stepfac_(4.1)
+  stepfac_(4.1),
+  splineDx_(0.01) // Recommendation from Roe & Brooks JMGM 2021
 {}
 
 void Action_Volmap::Help() const {
@@ -55,7 +61,7 @@ void Action_Volmap::Help() const {
 
 void Action_Volmap::RawHelp() const {
   mprintf("\t[out <filename>] <mask> [radscale <factor>] [stepfac <fac>]\n"
-          "\t[sphere] [radii {vdw | element}]\n"
+          "\t[sphere] [radii {vdw | element}] [splinedx <spacing>]\n"
           "\t[calcpeaks] [peakcut <cutoff>] [peakfile <xyzfile>]\n"
           "\t{ data <existing set> |\n"
           "\t  name <setname> <dx> [<dy> <dz>]\n"
@@ -82,6 +88,8 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d
     radscale_ = 0.5;
     stepfac_ = 1.0;
   }
+  //splineDx_ = actionArgs.getKeyDouble("splinedx", 1.0/5000.0);
+  splineDx_ = actionArgs.getKeyDouble("splinedx", 0.01);
   radscale_ = 1.0 / actionArgs.getKeyDouble("radscale", radscale_);
   stepfac_ = actionArgs.getKeyDouble("stepfac", stepfac_);
   std::string radarg = actionArgs.GetStringKey("radii");
@@ -214,7 +222,11 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d
     xmin_ = oxyz[0];
     ymin_ = oxyz[1];
     zmin_ = oxyz[2];
-  } 
+  }
+  // Warn for small grid spacing
+  if (dx_ < 0.4 || dy_ < 0.4 || dz_ < 0.4)
+    mprintf("Warning: Grid spacings smaller than 0.4 Ang. may be very slow;\n"
+            "Warning:  consider using larger grid spacings.\n");
   //std::string density = actionArgs.GetStringKey("density"); // FIXME obsolete?
   // Get the required mask
   std::string reqmask = actionArgs.GetMaskNext();
@@ -295,7 +307,26 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d
     mprintf("\tWhen smearing Gaussian, voxels farther than radii/2 will be skipped.\n");
   mprintf("\tDividing radii by %f\n", 1.0/radscale_);
   mprintf("\tFactor for determining number of bins to smear Gaussian is %f\n", stepfac_);
-
+# if defined(VOLMAP_USEEXP)
+  mprintf("\tUsing system exp() function for evaluating Gaussians.\n");
+//# elif defined(VOLMAP_USEFASTEXPS)
+//  mprintf("\tUsing exp() from N. Schraudolph, Neural Computation 11, 853–862 (1999).\n");
+//# elif defined(VOLMAP_USEFASTEXP64)
+//  mprintf("\tUsing 64 bit version of exp() from N. Schraudolph, Neural Computation 11, 853–862 (1999).\n");
+//# elif defined(VOLMAP_USEFASTEXPIEEE)
+//  mprintf("\tUsing exp() from Schraudolph & Malossi et al.\n");
+# else /* VOLMAP_USEEXP */
+  mprintf("\tExponential for Gaussians will be approximated using cubic splines\n"
+          "\t  with a spacing of %g Ang.\n", splineDx_);
+  mprintf("# Citation: Roe, D. R.; Brooks, B. R.; \"Improving the Speed of Volumetric\n"
+          "#           Density Map Generation via Cubic Spline Interpolation\".\n"
+          "#           Journal of Molecular Graphics and Modelling (2021).\n");
+# if defined(VOLMAP_USEACCURATE)
+  mprintf("\tSplines using more accurate but slower table lookup.\n");
+# elif defined(VOLMAP_USEXTABLE)
+  mprintf("\tSplines using less accurate lookup with tabled X values.\n");
+# endif
+# endif /* VOLMAP_USEEXP */
   if (outfile != 0)
     mprintf("\tDensity will wrtten to '%s'\n", outfile->DataFilename().full());
   mprintf("\tGrid dataset name is '%s'\n", grid_->legend());
@@ -381,8 +412,14 @@ Action::RetType Action_Volmap::Setup(ActionSetup& setup) {
   //mprintf("DEBUG: nx= %i  ny= %i  nz= %i\n", nxstep, nystep, nzstep);
   //mprintf("DEBUG: %g %g %g %g\n", maxx, maxy, maxz, maxDist);
   maxDist *= (-1.0 / (2.0 * maxRad * maxRad));
-  if (debug_ > 1)
-    mprintf("DEBUG: maxDist= %g\n", maxDist);
+  //mprintf("DEBUG: max= %g\n", maxDist);
+# ifndef VOLMAP_USEEXP
+  // Set up the interpolation table
+  if (table_.FillTable( exp, splineDx_, maxDist, 1.0 )) return Action::ERR;
+  table_.PrintMemUsage("\t");
+  if (debug_ > 0)
+    table_.PrintTableInfo("DEBUG: ");
+# endif
 
   if ((int)Atoms_.size() < densitymask_.Nselected())
     mprintf("Warning: %i atoms have 0.0 radii and will be skipped.\n",
@@ -511,9 +548,39 @@ Action::RetType Action_Volmap::DoAction(int frameNum, ActionFrame& frm) {
                 if (dist2 < rcut2) {
                   //mprintf("DEBUG: rhalf= %g  dist2= %g  exfac= %g  exp= %g\n", rhalf, dist2, exfac, exfac*dist2);
 #                 ifdef _OPENMP
+                  // NOTE: It is OK to call table_.Yval() here because in OpenMP
+                  //       local variables of called functions are private.
+#                 if defined(VOLMAP_USEEXP)
                   GRID_THREAD_[mythread].incrementBy(xval, yval, zval, norm * exp(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXPS)
+//                  GRID_THREAD[mythread].incrementBy(xval, yval, zval, norm * FASTEXPS(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXP64)
+//                  GRID_THREAD[mythread].incrementBy(xval, yval, zval, norm * fast_exps_64(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXPIEEE)
+//                  GRID_THREAD[mythread].incrementBy(xval, yval, zval, norm * fastexp::exp<double, fastexp::IEEE, 4>(exfac * dist2));
+#                 elif defined(VOLMAP_USEACCURATE)
+                  GRID_THREAD_[mythread].incrementBy(xval, yval, zval, norm * table_.Yval_accurate(exfac * dist2));
+#                 elif defined(VOLMAP_USEXTABLE)
+                  GRID_THREAD_[mythread].incrementBy(xval, yval, xval, norm * table_.Yval_xtable(exfac * dist2));
+#                 else
+                  GRID_THREAD_[mythread].incrementBy(xval, yval, zval, norm * table_.Yval(exfac * dist2));
+#                 endif
 #                 else /* OPENMP */
+#                 if defined(VOLMAP_USEEXP)
                   grid_->Increment(xval, yval, zval, norm * exp(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXPS)
+//                  grid_->Increment(xval, yval, zval, norm * FASTEXPS(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXP64)
+//                  grid_->Increment(xval, yval, zval, norm * fast_exps_64(exfac * dist2));
+//#                 elif defined(VOLMAP_USEFASTEXPIEEE)
+//                  grid_->Increment(xval, yval, zval, norm * fastexp::exp<double, fastexp::IEEE, 4>(exfac * dist2));
+#                 elif defined(VOLMAP_USEACCURATE)
+                  grid_->Increment(xval, yval, zval, norm * table_.Yval_accurate(exfac * dist2));
+#                 elif defined(VOLMAP_USEXTABLE)
+                  grid_->Increment(xval, yval, zval, norm * table_.Yval_xtable(exfac * dist2));
+#                 else
+                  grid_->Increment(xval, yval, zval, norm * table_.Yval(exfac * dist2));
+#                 endif
 #                 endif /* OPENMP */
                 }
               } // END loop over zval
@@ -531,6 +598,7 @@ Action::RetType Action_Volmap::DoAction(int frameNum, ActionFrame& frm) {
 }
 
 #ifdef MPI
+/** Sync grid to the master process. */
 int Action_Volmap::SyncAction() {
 # ifdef _OPENMP
   CombineGridThreads();

src/Action_Volmap.h

@@ -2,6 +2,7 @@
 #define INC_ACTION_VOLMAP_H
 #include "Action.h"
 #include "Grid.h"
+#include "SplineFxnTable.h"
 #ifdef VOLMAP_DOUBLE
 # define VOLMAP_DS_T DataSet_GridDbl
 # define VOLMAP_T double
@@ -10,6 +11,14 @@
 # define VOLMAP_T float
 #endif
 class VOLMAP_DS_T;
+/// Calculate atomic volumetric density maps from trajectory data.
+/** By default the grid type used is single-precision, mostly to save space.
+  * A double-precision grid can be used by compiling with the
+  * VOLMAP_DOUBLE define.
+  * Also by default the exp() function will be approximated with cubic spline
+  * interpolation. To use the system exp() function, compile with the
+  * VOLMAP_USEEXP define.
+  */
 class Action_Volmap : public Action {
   public:
     Action_Volmap();
@@ -51,6 +60,8 @@ class Action_Volmap : public Action {
     double radscale_;       ///< The scaling factor to divide all radii by
     double stepfac_;        ///< Factor for determining how many steps to smear Gaussian
     static const double sqrt_8_pi_cubed_;
+    SplineFxnTable table_;
+    double splineDx_;
 #   ifdef _OPENMP
     typedef std::vector< Grid<VOLMAP_T> > Garray;
     Garray GRID_THREAD_;

src/Ewald.cpp

@@ -13,18 +13,17 @@
 
 /// CONSTRUCTOR
 Ewald::Ewald() :
-  sumq_(0.0),
-  sumq2_(0.0),
   ew_coeff_(0.0),
   lw_coeff_(0.0),
   switch_width_(0.0),
   cutoff_(0.0),
   cut2_(0.0),
   cut2_0_(0.0),
   dsumTol_(0.0),
-  erfcTableDx_(0.0),
-  one_over_Dx_(0.0),
-  debug_(0)
+  debug_(0),
+  sumq_(0.0),
+  sumq2_(0.0),
+  Vdw_Recip_term_(0)
 {
 # ifdef DEBUG_EWALD
   // Save fractional translations for 1 cell in each direction (and primary cell).
@@ -91,44 +90,9 @@ double Ewald::erfc_func(double xIn) {
   return erfc;
 }
 
-// Ewald::FillErfcTable()
-void Ewald::FillErfcTable(double cutoffIn, double dxdr) {
-  one_over_Dx_ = 1.0 / erfcTableDx_;
-  unsigned int erfcTableSize = (unsigned int)(dxdr * one_over_Dx_ * cutoffIn * 1.5);
-  Darray erfc_X, erfc_Y;
-  erfc_X.reserve( erfcTableSize );
-  erfc_Y.reserve( erfcTableSize );
-  // Save X and Y values so we can calc the spline coefficients
-  double xval = 0.0;
-  for (unsigned int i = 0; i != erfcTableSize; i++) {
-    double yval = erfc_func( xval );
-    erfc_X.push_back( xval );
-    erfc_Y.push_back( yval );
-    xval += erfcTableDx_;
-  }
-  Spline cspline;
-  cspline.CubicSpline_Coeff(erfc_X, erfc_Y);
-  erfc_X.clear();
-  // Store values in Spline table
-  erfc_table_.reserve( erfcTableSize * 4 ); // Y B C D
-  for (unsigned int i = 0; i != erfcTableSize; i++) {
-    erfc_table_.push_back( erfc_Y[i] );
-    erfc_table_.push_back( cspline.B_coeff()[i] );
-    erfc_table_.push_back( cspline.C_coeff()[i] );
-    erfc_table_.push_back( cspline.D_coeff()[i] );
-  }
-  // Memory saved Y values plus spline B, C, and D coefficient arrays.
-  mprintf("\tMemory used by Erfc table and splines: %s\n",
-          ByteString(erfc_table_.size() * sizeof(double), BYTE_DECIMAL).c_str());
-}
-
 // Ewald::ERFC()
 double Ewald::ERFC(double xIn) const {
-  int xidx = ((int)(one_over_Dx_ * xIn));
-  double dx = xIn - ((double)xidx * erfcTableDx_);
-  xidx *= 4;
-  return erfc_table_[xidx] + 
-         dx*(erfc_table_[xidx+1] + dx*(erfc_table_[xidx+2] + dx*erfc_table_[xidx+3]));
+  return table_.Yval( xIn);
 }
 
 /** Determine Ewald coefficient from cutoff and direct sum tolerance.
@@ -200,7 +164,7 @@ void Ewald::CalculateC6params(Topology const& topIn, AtomMask const& maskIn) {
 }
 
 /** Set up exclusion lists for selected atoms. */
-void Ewald::SetupExcluded(Topology const& topIn, AtomMask const& maskIn)
+void Ewald::SetupExclusionList(Topology const& topIn, AtomMask const& maskIn)
 {
   // Use distance of 4 (up to dihedrals)
   if (Excluded_.SetupExcluded(topIn.Atoms(), maskIn, 4,
@@ -223,7 +187,7 @@ int Ewald::CheckInput(Box const& boxIn, int debugIn, double cutoffIn, double dsu
   ew_coeff_ = ew_coeffIn;
   lw_coeff_ = lw_coeffIn;
   switch_width_ = switch_widthIn;
-  erfcTableDx_ = erfcTableDxIn;
+  double erfcTableDx = erfcTableDxIn;
   // Check input
   if (cutoff_ < Constants::SMALL) {
     mprinterr("Error: Direct space cutoff (%g) is too small.\n", cutoff_);
@@ -252,9 +216,14 @@ int Ewald::CheckInput(Box const& boxIn, int debugIn, double cutoffIn, double dsu
     dsumTol_ = 1E-5;
   if (DABS(ew_coeff_) < Constants::SMALL)
     ew_coeff_ = FindEwaldCoefficient( cutoff_, dsumTol_ );
-  if (erfcTableDx_ <= 0.0) erfcTableDx_ = 1.0 / 5000;
+  if (erfcTableDx <= 0.0) erfcTableDx = 1.0 / 5000;
   // TODO make this optional
-  FillErfcTable( cutoff_, ew_coeff_ ); 
+  if (table_.FillTable( erfc_func, erfcTableDx, 0.0, cutoff_*ew_coeff_*1.5 )) {
+    mprinterr("Error: Could not set up spline table for ERFC\n");
+    return 1;
+  }
+  table_.PrintMemUsage("\t");
+  table_.PrintTableInfo("\t");
   // TODO do for C6 as well
   // TODO for C6 correction term
   if (lw_coeff_ < 0.0)