Port potential Databox enhancements into Spiner from Singe -- transformations

spiner/piecewise_grid_1d.hpp

@@ -35,6 +35,7 @@
 namespace Spiner {
 
 template <typename T = Real, int NGRIDSMAX = 5,
+          typename Transform = TransformLinear,
           typename =
               typename std::enable_if<std::is_arithmetic<T>::value, bool>::type>
 class PiecewiseGrid1D {
@@ -46,7 +47,7 @@ class PiecewiseGrid1D {
   // This is functionally equivalent because grids_ will
   // be initialized to default values
   PORTABLE_INLINE_FUNCTION PiecewiseGrid1D() {}
-  PiecewiseGrid1D(const std::vector<RegularGrid1D<T>> grids) {
+  PiecewiseGrid1D(const std::vector<RegularGrid1D<T, Transform>> grids) {
     NGRIDS_ = grids.size();
     PORTABLE_ALWAYS_REQUIRE(
         NGRIDS_ <= NGRIDSMAX,
@@ -65,8 +66,8 @@ class PiecewiseGrid1D {
       }
     }
   }
-  PiecewiseGrid1D(std::initializer_list<RegularGrid1D<T>> grids)
-      : PiecewiseGrid1D(std::vector<RegularGrid1D<T>>(grids)) {}
+  PiecewiseGrid1D(std::initializer_list<RegularGrid1D<T, Transform>> grids)
+      : PiecewiseGrid1D(std::vector<RegularGrid1D<T, Transform>>(grids)) {}
 
   template <typename F>
   PORTABLE_INLINE_FUNCTION int findGrid(const F &direction) const {
@@ -124,14 +125,14 @@ class PiecewiseGrid1D {
   }
 
   PORTABLE_INLINE_FUNCTION
-  bool operator==(const PiecewiseGrid1D<T, NGRIDSMAX> &other) const {
+  bool operator==(const PiecewiseGrid1D<T, NGRIDSMAX, Transform> &other) const {
     for (int ig = 0; ig < NGRIDS_; ++ig) {
       if (grids_[ig] != other.grids_[ig]) return false;
     }
     return true;
   }
   PORTABLE_INLINE_FUNCTION
-  bool operator!=(const PiecewiseGrid1D<T, NGRIDSMAX> &other) const {
+  bool operator!=(const PiecewiseGrid1D<T, NGRIDSMAX, Transform> &other) const {
     return !(*this == other);
   }
   PORTABLE_INLINE_FUNCTION T min() const { return grids_[0].min(); }
@@ -225,7 +226,7 @@ class PiecewiseGrid1D {
            SP5::H1D::GRID_FORMAT[1];
   }
 
-  RegularGrid1D<T> grids_[NGRIDSMAX];
+  RegularGrid1D<T, Transform> grids_[NGRIDSMAX];
   int pointTotals_[NGRIDSMAX];
   int NGRIDS_;
 };

spiner/regular_grid_1d.hpp

@@ -31,6 +31,7 @@
 #include "ports-of-call/portable_errors.hpp"
 #include "sp5.hpp"
 #include "spiner_types.hpp"
+#include "transformations.hpp"
 
 namespace Spiner {
 
@@ -44,7 +45,7 @@ struct weights_t {
   }
 };
 
-template <typename T = Real,
+template <typename T = Real, typename Transform = TransformLinear,
           typename std::enable_if<std::is_arithmetic<T>::value, bool>::type =
               true>
 class RegularGrid1D {
@@ -55,11 +56,11 @@ class RegularGrid1D {
 
   // Constructors
   PORTABLE_INLINE_FUNCTION RegularGrid1D()
-      : min_(rNaN), max_(rNaN), dx_(rNaN), idx_(rNaN), N_(iNaN) {}
-  PORTABLE_INLINE_FUNCTION RegularGrid1D(T min, T max, size_t N)
-      : min_(min), max_(max), dx_((max - min) / ((T)(N - 1))), idx_(1 / dx_),
-        N_(N) {
-    PORTABLE_ALWAYS_REQUIRE(min_ < max_ && N_ > 0, "Valid grid");
+      : umin_(rNaN), umax_(rNaN), du_(rNaN), inv_du_(rNaN), N_(iNaN) {}
+  PORTABLE_INLINE_FUNCTION RegularGrid1D(T xmin, T xmax, size_t N)
+      : umin_(Transform::forward(xmin)), umax_(Transform::forward(xmax)),
+        du_((umax_ - umin_) / ((T)(N - 1))), inv_du_(1 / du_), N_(N) {
+    PORTABLE_ALWAYS_REQUIRE(umin_ < umax_ && N_ > 0, "Valid grid");
   }
 
   // Forces x in the interval
@@ -71,18 +72,27 @@ class RegularGrid1D {
     return ix;
   }
 
-  // Gets real value at index
-  PORTABLE_INLINE_FUNCTION T x(const int i) const { return i * dx_ + min_; }
+  // Translate between u (transformed variable) coordinate and index
+  PORTABLE_INLINE_FUNCTION T u(const int i) const { return i * du_ + umin_; }
+  PORTABLE_INLINE_FUNCTION int index_u(const T u) const {
+    return bound(inv_du_ * (u - umin_));
+  }
+
+  // Translate between x coordinate and index
+  PORTABLE_INLINE_FUNCTION T x(const int i) const {
+    return Transform::reverse(u(i));
+  }
   PORTABLE_INLINE_FUNCTION int index(const T x) const {
-    return bound(idx_ * (x - min_));
+    return index_u(Transform::forward(x));
   }
 
   // Returns closest index and weights for interpolation
   PORTABLE_INLINE_FUNCTION void weights(const T &x, int &ix,
                                         weights_t<T> &w) const {
-    ix = index(x);
-    const auto floor = static_cast<T>(ix) * dx_ + min_;
-    w[1] = idx_ * (x - floor);
+    const T u = Transform::forward(x);
+    ix = index_u(u);
+    const auto floor = static_cast<T>(ix) * du_ + umin_;
+    w[1] = inv_du_ * (u - floor);
     w[0] = (1. - w[1]);
   }
 
@@ -98,20 +108,21 @@ class RegularGrid1D {
   // utitilies
   PORTABLE_INLINE_FUNCTION bool
   operator==(const RegularGrid1D<T> &other) const {
-    return (min_ == other.min_ && max_ == other.max_ && dx_ == other.dx_ &&
-            idx_ == other.idx_ && N_ == other.N_);
+    return (umin_ == other.umin_ && umax_ == other.umax_ && du_ == other.du_ &&
+            inv_du_ == other.inv_du_ && N_ == other.N_);
   }
   PORTABLE_INLINE_FUNCTION bool
   operator!=(const RegularGrid1D<T> &other) const {
     return !(*this == other);
   }
-  PORTABLE_INLINE_FUNCTION T min() const { return min_; }
-  PORTABLE_INLINE_FUNCTION T max() const { return max_; }
-  PORTABLE_INLINE_FUNCTION T dx() const { return dx_; }
+  PORTABLE_INLINE_FUNCTION T umin() const { return umin_; }
+  PORTABLE_INLINE_FUNCTION T umax() const { return umax_; }
+  PORTABLE_INLINE_FUNCTION T min() const { return Transform::reverse(umin_); }
+  PORTABLE_INLINE_FUNCTION T max() const { return Transform::reverse(umax_); }
   PORTABLE_INLINE_FUNCTION size_t nPoints() const { return N_; }
   PORTABLE_INLINE_FUNCTION bool isnan() const {
-    return (std::isnan(min_) || std::isnan(max_) || std::isnan(dx_) ||
-            std::isnan(idx_) || std::isnan((T)N_));
+    return (std::isnan(umin_) || std::isnan(umax_) || std::isnan(du_) ||
+            std::isnan(inv_du_) || std::isnan((T)N_));
   }
   PORTABLE_INLINE_FUNCTION bool isWellFormed() const { return !isnan(); }
 
@@ -123,7 +134,7 @@ class RegularGrid1D {
     auto H5T_T =
         std::is_same<T, double>::value ? H5T_NATIVE_DOUBLE : H5T_NATIVE_FLOAT;
     herr_t status;
-    T range[] = {min_, max_, dx_};
+    T range[] = {umin_, umax_, du_};
     hsize_t range_dims[] = {3};
     int n = static_cast<int>(N_);
     status = H5LTmake_dataset(loc, name.c_str(), SP5::RG1D::RANGE_RANK,
@@ -144,19 +155,19 @@ class RegularGrid1D {
     T range[3];
     int n;
     status = H5LTread_dataset(loc, name.c_str(), H5T_T, range);
-    min_ = range[0];
-    max_ = range[1];
-    dx_ = range[2];
-    idx_ = 1. / dx_;
+    umin_ = range[0];
+    umax_ = range[1];
+    du_ = range[2];
+    inv_du_ = 1. / du_;
     status += H5LTget_attribute_int(loc, name.c_str(), SP5::RG1D::N, &n);
     N_ = n;
     return status;
   }
 #endif
 
  private:
-  T min_, max_;
-  T dx_, idx_;
+  T umin_, umax_;
+  T du_, inv_du_;
   size_t N_;
 };
 

spiner/transformations.hpp

@@ -0,0 +1,60 @@
+#ifndef _SPINER_TRANSFORM_HPP_
+#define _SPINER_TRANSFORM_HPP_
+//======================================================================
+// © (or copyright) 2019-2021. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+//======================================================================
+
+#include "ports-of-call/portability.hpp"
+
+#include <cmath>
+#include <limits>
+
+namespace Spiner {
+
+// Note on notation:
+// -- "real" space is called x
+// -- "transformed" space is called u
+// -- u = forward(x)
+// -- x = reverse(u)
+
+// linear transformation (aka no-op): y = x
+struct TransformLinear {
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION static T forward(const T x) {
+    return x;
+  }
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION static T reverse(const T u) {
+    return u;
+  }
+};
+
+// logarithmic transformation: y = log(x + small)
+struct TransformLogarithmic {
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION static T forward(const T x) {
+    return std::log(x + std::numeric_limits<T>::denorm_min());
+  }
+  template <typename T>
+  PORTABLE_INLINE_FUNCTION static T reverse(const T u) {
+    return std::exp(u) - std::numeric_limits<T>::denorm_min();
+  }
+};
+
+// TODO: log_NQT and arcsinh_NQT, but these require adding a dependency on
+//       https://github.com/lanl/not-quite-transcendental.  I may leave this for
+//       Jonah ;-)
+
+} // namespace Spiner
+#endif // _SPINER_TRANSFORM_HPP_

test/CMakeLists.txt

@@ -57,7 +57,7 @@ target_link_libraries(spiner_test
     ${spinerTest_POPULATED_TARGETS}
 )
 # add unit tests
-add_executable(test.bin test.cpp)
+add_executable(test.bin test.cpp regular_grid_1d.cpp transformations.cpp)
 target_link_libraries(test.bin
   PRIVATE
     spiner::spiner

test/benchmark.cpp

@@ -74,8 +74,10 @@ int main(int argc, char *argv[]) {
 
     RegularGrid1D gcoarse(xmin, xmax, ncoarse);
     std::vector<RegularGrid1D> gfine;
+    std::vector<Real> dxfine;
     for (const auto &n : nfine) {
       gfine.push_back(RegularGrid1D(xmin, xmax, n));
+      dxfine.push_back((xmax - xmin) / n);
     }
 
     std::cout << "# ncoarse = " << ncoarse << std::endl;
@@ -98,7 +100,8 @@ int main(int argc, char *argv[]) {
     for (int ifine = 0; ifine < nfine.size(); ++ifine) {
       auto n = nfine[ifine];
       auto g = gfine[ifine];
-      Real d3x = g.dx() * g.dx() * g.dx();
+      auto dx = dxfine[ifine];
+      Real d3x = dx * dx * dx;
       Real L2_error;
 #ifdef PORTABILITY_STRATEGY_KOKKOS
       Kokkos::fence();

test/regular_grid_1d.cpp

@@ -0,0 +1,88 @@
+// © (or copyright) 2019-2021. Triad National Security, LLC. All rights
+// reserved.  This program was produced under U.S. Government contract
+// 89233218CNA000001 for Los Alamos National Laboratory (LANL), which is
+// operated by Triad National Security, LLC for the U.S.  Department of
+// Energy/National Nuclear Security Administration. All rights in the
+// program are reserved by Triad National Security, LLC, and the
+// U.S. Department of Energy/National Nuclear Security
+// Administration. The Government is granted for itself and others acting
+// on its behalf a nonexclusive, paid-up, irrevocable worldwide license
+// in this material to reproduce, prepare derivative works, distribute
+// copies to the public, perform publicly and display publicly, and to
+// permit others to do so.
+
+#include <ports-of-call/portability.hpp>
+#include <spiner/regular_grid_1d.hpp>
+
+#include <catch2/catch_session.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include <cmath>
+
+using Catch::Matchers::WithinRel;
+
+TEST_CASE("RegularGrid1D", "[RegularGrid1D]") {
+  SECTION("A regular grid 1d emits appropriate metadata") {
+    constexpr Real min = -1;
+    constexpr Real max = 1;
+    constexpr size_t N = 10;
+    Spiner::RegularGrid1D<Real> g(min, max, N);
+    REQUIRE(g.min() == min);
+    REQUIRE(g.max() == max);
+    REQUIRE(g.nPoints() == N);
+  }
+}
+
+TEST_CASE("RegularGrid1D with transformations",
+          "[RegularGrid1D][transformations]") {
+  constexpr double min = 1;
+  constexpr double max = 1024;
+  constexpr size_t N = 11;
+  Spiner::RegularGrid1D<double, Spiner::TransformLinear> glin(min, max, N);
+  Spiner::RegularGrid1D<double, Spiner::TransformLogarithmic> glog(min, max, N);
+  REQUIRE(glin.min() == min);
+  REQUIRE(glog.min() == min);
+  REQUIRE(glin.max() == max);
+  REQUIRE(glog.max() == max);
+  REQUIRE(glin.nPoints() == N);
+  REQUIRE(glog.nPoints() == N);
+
+  // Check all fixed points (lin)
+  for (std::size_t n = 0; n < N; ++n) {
+    const double xx = 102.3 * double(n) + 1;
+    const double uu = xx;
+    CHECK_THAT(glin.u(n), WithinRel(uu, 1.0e-12));
+    CHECK_THAT(glin.x(n), WithinRel(xx, 1.0e-12));
+  }
+
+  // Check all fixed points (log)
+  for (std::size_t n = 0; n < N; ++n) {
+    const double xx = std::pow(double(2), n);
+    const double uu = std::log(xx);
+    CHECK_THAT(glog.u(n), WithinRel(uu, 1.0e-12));
+    CHECK_THAT(glog.x(n), WithinRel(xx, 1.0e-12));
+  }
+
+  // Check weights (lin)
+  for (std::size_t n = 1; n < N; ++n) {
+    const double xlin = 0.5 * (glin.x(n - 1) + glin.x(n));
+    int ilin;
+    Spiner::weights_t<double> wlin;
+    glin.weights(xlin, ilin, wlin);
+    CHECK(ilin == n - 1);
+    CHECK_THAT(wlin.first, WithinRel(0.5));
+    CHECK_THAT(wlin.second, WithinRel(0.5));
+  }
+
+  // Check weights (log)
+  for (std::size_t n = 1; n < N; ++n) {
+    const double xlog = std::sqrt(glog.x(n - 1) * glog.x(n));
+    int ilog;
+    Spiner::weights_t<double> wlog;
+    glog.weights(xlog, ilog, wlog);
+    CHECK(ilog == n - 1);
+    CHECK_THAT(wlog.first, WithinRel(0.5));
+    CHECK_THAT(wlog.second, WithinRel(0.5));
+  }
+}