Modify RegularGrid1D to allow for transformations.

spiner/regular_grid_1d.hpp

@@ -44,7 +44,23 @@ struct weights_t {
   }
 };
 
+// TODO: Do transformations need state?
+//    -- My thinking is that transformations generally won't need state, because it's things like
+//       y = x or y = log(x).
+//    -- Based on this, I would invoke the transform as Transform::forward(x) or
+//       Transform::reverse(u).
+//    -- But this means that transformations _cannot_ have state if that use-case comes up in the
+//       future.
+//    -- The other option would be to have the transformation object be passed into the
+//       constructor.
+//       -- We could default to Transform(), but that means all transformations _must_ have a
+//          default constructor even if it makes no sense.
+//       -- We could have no default, but then users are now required to update all their
+//          constructors to add TransformLinear() as an argument, and we can't hide the addition of
+//          the template.
+
 template <typename T = Real,
+          typename Transform = TransformLinear,
           typename std::enable_if<std::is_arithmetic<T>::value, bool>::type =
               true>
 class RegularGrid1D {
@@ -55,11 +71,14 @@ 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 +90,25 @@ 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,23 +124,38 @@ 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_; }
+  // TODO: This way of constructing min() and max() implicitly asserts that the u-space
+  //       representation is the "ground truth" and the x-space representation is merely derived
+  //       from there.  We could change this and assert that the x-space representation is the
+  //       "ground truth", but we would have to make some changes.  This arises because it's not
+  //       guaranteed that every transformation is _exactly_ one-to-one reversible, so you could
+  //       introduce a small gap between xmin and reverse(umin) or between forward(xmin) and umin
+  //       (and similarly for max).
+  // TODO: Should umin() and umax() not be publicly exposed?  If so, they shouldn't exist because
+  //       they're only public getters for the private umin_ / umax_ variables.
+  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_); }
+  // TODO: Should there be a public du() function?
+  // TODO: dx() is now ill-defined -- do we need it?
+  //PORTABLE_INLINE_FUNCTION T dx() const { return dx_; }
   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(); }
 
+  // TODO: I made very simple, naive changes to saveHDF and loadHDF, because it's not clear to me
+  //       if something better should be done.
 #ifdef SPINER_USE_HDF
   inline herr_t saveHDF(hid_t loc, const std::string &name) const {
     static_assert(
@@ -123,7 +164,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 +185,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_;
 };