From 5ce03cbdd62b5c3469ed4232052be438c8aa55de Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Fri, 13 Sep 2024 16:54:32 -0600 Subject: [PATCH] Modify RegularGrid1D to allow for transformations. --- spiner/regular_grid_1d.hpp | 91 +++++++++++++++++++++++++++----------- 1 file changed, 66 insertions(+), 25 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index d69076d65..6753d8fbb 100644 --- a/spiner/regular_grid_1d.hpp +++ b/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 ::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 &w) const { - ix = index(x); - const auto floor = static_cast(ix) * dx_ + min_; - w[1] = idx_ * (x - floor); + const T u = Transform::forward(x); + ix = index_u(u); + const auto floor = static_cast(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 &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 &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::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(N_); status = H5LTmake_dataset(loc, name.c_str(), SP5::RG1D::RANGE_RANK, @@ -144,10 +185,10 @@ 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; @@ -155,8 +196,8 @@ class RegularGrid1D { #endif private: - T min_, max_; - T dx_, idx_; + T umin_, umax_; + T du_, inv_du_; size_t N_; };