Merge pull request #4875 from PDoakORNL/introducing_ParseGridInput

Introducing ParseGridInput

src/Estimators/CMakeLists.txt

@@ -18,6 +18,7 @@ set(QMCEST_SRC
     RMCLocalEnergyEstimator.cpp
     EstimatorInput.cpp
     EstimatorManagerInput.cpp
+    ParseGridInput.cpp
     ScalarEstimatorInputs.cpp
     SpinDensityInput.cpp
     MomentumDistributionInput.cpp

src/Estimators/ParseGridInput.cpp

@@ -0,0 +1,225 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// This file is distributed under the University of Illinois/NCSA Open Source License.
+// See LICENSE file in top directory for details.
+//
+// Copyright (c) 2023 QMCPACK developers.
+//
+// File developed by: Peter  W. Doak, doakpw@ornl.gov, Oak Ridge National Lab
+//
+// File created by: Peter  W. Doak, doakpw@ornl.gov, Oak Ridge National Lab
+//////////////////////////////////////////////////////////////////////////////////////
+
+#include "ParseGridInput.hpp"
+
+#include <cmath>
+
+#include <ModernStringUtils.hpp>
+#include "Message/UniformCommunicateError.h"
+
+namespace qmcplusplus
+{
+
+template<typename REAL>
+AxisGrid<REAL> parseGridInput(std::istringstream& grid_input_stream)
+{
+  using namespace modernstrutil;
+  const REAL utol = 1e-5;
+
+  std::string grid_line;
+  std::getline(grid_input_stream, grid_line);
+  // This refactored from QMCHamiltonian/SpaceGrid::initializeRectilinear
+  std::vector<std::string_view> tokens = split(grid_line, " ");
+
+  //  count the number of intervals
+  int nintervals = 0;
+  for (int i = 0; i < tokens.size(); i++)
+    if (tokens[i][0] != '(')
+      nintervals++;
+  nintervals--;
+  //  allocate temporary interval variables
+  AxisGrid<REAL> agr;
+  agr.ndom_int.resize(nintervals);
+  agr.du_int.resize(nintervals);
+  agr.ndu_int.resize(nintervals);
+  //  determine number of domains in each interval and the width of each domain
+  REAL u1 = string2Real<REAL>(tokens[0]);
+
+  agr.umin = u1;
+  if (std::abs(u1) > 1.0000001)
+  {
+    std::ostringstream msg;
+    msg << "  interval endparticles cannot be greater than " << 1 << '\n' << "  endpoint provided: " << u1;
+    throw UniformCommunicateError(msg.str());
+  }
+  bool is_int        = false;
+  bool has_paren_val = false;
+  REAL du_i;
+  int ndom_i   = 1;
+  int interval = -1;
+  for (int i = 1; i < tokens.size(); i++)
+  {
+    if (tokens[i][0] != '(')
+    {
+      REAL u2  = string2Real<REAL>(tokens[i]);
+      agr.umax = u2;
+      if (!has_paren_val)
+        du_i = u2 - u1;
+      has_paren_val = false;
+      interval++;
+      if (u2 < u1)
+      {
+        std::ostringstream msg;
+
+        msg << "  interval (" << u1 << "," << u2 << ") is negative" << std::endl;
+        throw UniformCommunicateError(msg.str());
+      }
+      if (std::abs(u2) > 1.0000001)
+      {
+        std::ostringstream msg;
+
+        msg << "  interval endparticles cannot be greater than " << 1 << std::endl;
+        msg << "  endpoint provided: " << u2 << std::endl;
+        throw UniformCommunicateError(msg.str());
+      }
+      if (is_int)
+      {
+        agr.du_int[interval]   = (u2 - u1) / ndom_i;
+        agr.ndom_int[interval] = ndom_i;
+      }
+      else
+      {
+        agr.du_int[interval]   = du_i;
+        agr.ndom_int[interval] = std::floor((u2 - u1) / du_i + .5);
+        if (std::abs(u2 - u1 - du_i * agr.ndom_int[interval]) > utol)
+        {
+          std::ostringstream msg;
+
+          msg << "  interval (" << u1 << "," << u2 << ") not divisible by du=" << du_i << std::endl;
+          throw UniformCommunicateError(msg.str());
+        }
+      }
+      u1 = u2;
+    }
+    else
+    {
+      has_paren_val = true;
+      std::string paren_val{tokens[i].substr(1, tokens[i].length() - 2)};
+      is_int = tokens[i].find(".") == std::string::npos;
+      if (is_int)
+      {
+        ndom_i = string2Int<int>(paren_val);
+        du_i   = -1.0;
+      }
+      else
+      {
+        ndom_i = 0;
+        du_i   = string2Real<REAL>(paren_val);
+      }
+    }
+  }
+  // find the smallest domain width
+  REAL du_min = 1.0;
+  for (int i = 0; i < agr.du_int.size(); i++)
+    du_min = std::min(du_min, agr.du_int[i]);
+  agr.odu = 1.0 / du_min;
+  // make sure it divides into all other domain widths
+  for (int i = 0; i < agr.du_int.size(); i++)
+  {
+    agr.ndu_int[i] = floor(agr.du_int[i] / du_min + .5);
+    if (std::abs(agr.du_int[i] - agr.ndu_int[i] * du_min) > utol)
+    {
+      std::ostringstream msg;
+      msg << "interval " << i + 1 << " of axis is not divisible by smallest subinterval " << du_min;
+      throw UniformCommunicateError(msg.str());
+    }
+  }
+  // set up the interval map such that gmap[u/du]==domain index
+  agr.gmap.resize(floor((agr.umax - agr.umin) * agr.odu + .5));
+  int n  = 0;
+  int nd = -1;
+  for (int i = 0; i < agr.ndom_int.size(); i++)
+    for (int j = 0; j < agr.ndom_int[i]; j++)
+    {
+      nd++;
+      for (int k = 0; k < agr.ndu_int[i]; k++)
+      {
+        //app_log()<<"        accessing gmap "<<iaxis<<" "<<n<< std::endl;
+        agr.gmap[n] = nd;
+        n++;
+      }
+    }
+  agr.dimensions = nd + 1;
+  //end read in the grid contents
+  //save interval width information
+  int ndom_tot = 0;
+  for (int i = 0; i < agr.ndom_int.size(); i++)
+    ndom_tot += agr.ndom_int[i];
+  agr.ndu_per_interval.resize(ndom_tot);
+  int idom = 0;
+  for (int i = 0; i < agr.ndom_int.size(); i++)
+    for (int ii = 0; ii < agr.ndom_int[i]; ii++)
+    {
+      agr.ndu_per_interval[idom] = agr.ndu_int[i];
+      idom++;
+    }
+  return agr;
+}
+
+template<typename REAL>
+AxisGrid<REAL>::AxisGrid(std::vector<int>&& rhs_ndom_int,
+                         std::vector<int>&& rhs_ndu_int,
+                         std::vector<REAL>&& rhs_du_int,
+                         REAL rhs_umin,
+                         REAL rhs_umax,
+                         REAL rhs_odu,
+                         std::vector<int>&& rhs_gmap,
+                         std::vector<int>&& rhs_ndu_per_interval,
+                         int rhs_dimensions)
+{
+  ndom_int         = rhs_ndom_int;
+  ndu_int          = rhs_ndu_int;
+  du_int           = rhs_du_int;
+  umin             = rhs_umin;
+  umax             = rhs_umax;
+  odu              = rhs_odu;
+  gmap             = rhs_gmap;
+  ndu_per_interval = rhs_ndu_per_interval;
+  dimensions       = rhs_dimensions;
+}
+
+template<typename REAL>
+AxisGrid<REAL>::AxisGrid(const AxisGrid& rhs)
+{
+  ndom_int         = rhs.ndom_int;
+  ndu_int          = rhs.ndu_int;
+  du_int           = rhs.du_int;
+  umin             = rhs.umin;
+  umax             = rhs.umax;
+  odu              = rhs.odu;
+  gmap             = rhs.gmap;
+  ndu_per_interval = rhs.ndu_per_interval;
+  dimensions       = rhs.dimensions;
+}
+
+template<typename REAL>
+AxisGrid<REAL>& AxisGrid<REAL>::operator=(AxisGrid<REAL> rhs)
+{
+  std::swap(ndom_int, rhs.ndom_int);
+  std::swap(ndu_int, rhs.ndu_int);
+  std::swap(du_int, rhs.du_int);
+  std::swap(umin, rhs.umin);
+  std::swap(umax, rhs.umax);
+  std::swap(odu, rhs.odu);
+  std::swap(gmap, rhs.gmap);
+  std::swap(ndu_per_interval, rhs.ndu_per_interval);
+  std::swap(dimensions, rhs.dimensions);
+  return *this;
+}
+
+
+// explicit instantiation
+template struct AxisGrid<float>;
+template struct AxisGrid<double>;
+template AxisGrid<double> parseGridInput<double>(std::istringstream& grid_input_stream);
+template AxisGrid<float> parseGridInput<float>(std::istringstream& grid_input_stream);
+} // namespace qmcplusplus

src/Estimators/ParseGridInput.hpp

@@ -0,0 +1,89 @@
+//////////////////////////////////////////////////////////////////////////////////////
+// This file is distributed under the University of Illinois/NCSA Open Source License.
+// See LICENSE file in top directory for details.
+//
+// Copyright (c) 2023 QMCPACK developers.
+//
+// File developed by: Peter  W. Doak, doakpw@ornl.gov, Oak Ridge National Lab
+//
+// File created by: Peter  W. Doak, doakpw@ornl.gov, Oak Ridge National Lab
+//////////////////////////////////////////////////////////////////////////////////////
+
+/** \file
+ *  The AxisGrid data structure and the ParseGridInput factor parsing in a manner usable in acustom handler
+ *  method of any InputSection subtype. Right now it is only required by
+ *  SpaceGridInput::SpaceGridAxisInputSection::setFromStreamCustom.
+ *  This takes the one dimensional grid input specified for EnergyDensityEstimator (described in the for
+ *  the legacy energy density estimator). To the suprising complex input variables required by SpaceGrid.
+ *  AxisGrid should be subjected to analysis at a later date for which of these data elements are actually
+ *  required to successfully construct a NESpaceGrid instance.
+ */
+
+#ifndef QMCPLUSPLUS_PARSE_GRID_INPUT_HPP
+#define QMCPLUSPLUS_PARSE_GRID_INPUT_HPP
+#include <vector>
+#include <string>
+#include <sstream>
+#include <iostream>
+
+namespace qmcplusplus
+{
+
+template<typename REAL>
+struct AxisGrid
+{
+  AxisGrid() = default;
+
+  /** "aggregate" style constructor
+   */
+  AxisGrid(std::vector<int>&& rhs_ndom_int,
+           std::vector<int>&& rhs_ndu_int,
+           std::vector<REAL>&& rhs_du_int,
+           REAL rhs_umin,
+           REAL rhs_umax,
+           REAL rhs_odu,
+           std::vector<int>&& rhs_gmap,
+           std::vector<int>&& rhs_ndu_per_interval,
+           int rhs_dimensions);
+
+  AxisGrid(const AxisGrid& rhs);
+
+  AxisGrid& operator=(AxisGrid rhs);
+
+  /** equality operator, all values must be equal.
+   *  In C++20 this will just be the defaulted operator==.
+   */
+  bool operator==(const AxisGrid& ag) const
+  {
+    return (ag.ndom_int == ndom_int && ag.ndu_int == ndu_int && ag.du_int == du_int && ag.umin == umin &&
+            ag.umax == umax && ag.odu == odu && ag.gmap == gmap && ag.ndu_per_interval == ndu_per_interval &&
+            ag.dimensions == dimensions);
+  }
+
+  std::vector<int> ndom_int;
+  std::vector<int> ndu_int;
+  std::vector<REAL> du_int;
+  REAL umin;
+  REAL umax;
+  REAL odu;
+  std::vector<int> gmap;
+  std::vector<int> ndu_per_interval;
+  int dimensions;
+};
+
+/** Parses the one dimensional grid specification format originated for EnergyDensity
+ *  This has be refactored from QMCHamiltonian/SpaceGrid.cpp.
+ *  My evidence for correctness is that it can pass the unit test I wrote, there was no legacy test.
+ */
+template<typename REAL>
+AxisGrid<REAL> parseGridInput(std::istringstream& grid_input_stream);
+
+// explicit instantiation declaration.
+extern template struct AxisGrid<float>;
+extern template struct AxisGrid<double>;
+extern template AxisGrid<double> parseGridInput<double>(std::istringstream& grid_input_stream);
+extern template AxisGrid<float> parseGridInput<float>(std::istringstream& grid_input_stream);
+
+} // namespace qmcplusplus
+
+#endif

src/Estimators/tests/CMakeLists.txt

@@ -42,6 +42,7 @@ set(SRCS
     test_EstimatorManagerCrowd.cpp
     test_MagnetizationDensityInput.cpp
     test_MagnetizationDensity.cpp
+    test_ParseGridInput.cpp
     test_PerParticleHamiltonianLogger.cpp
     test_ReferencePointsInput.cpp
     test_ReferencePoints.cpp