From 216d115186d061b7455a10f2f41cfd3e4f269d3e Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 26 Aug 2024 16:24:25 -0600 Subject: [PATCH 01/48] Copy databox_wrapper.hpp from Singe as a starting point for discussion. --- spiner/databox_wrapper.hpp | 72 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 spiner/databox_wrapper.hpp diff --git a/spiner/databox_wrapper.hpp b/spiner/databox_wrapper.hpp new file mode 100644 index 000000000..375d1d622 --- /dev/null +++ b/spiner/databox_wrapper.hpp @@ -0,0 +1,72 @@ +#ifndef SINGE_UTIL_DATABOXWRAPPER_HPP +#define SINGE_UTIL_DATABOXWRAPPER_HPP + +#include "singe/util/databox.hpp" +#include "singe/util/portability_macros.hpp" + +#include + +namespace singe { +namespace util { + +// ================================================================================================ + +// a wrapper to (a) adapt the interface and (b) handle extrapolation off the table +// -- note that we assume this is a 1D databox +template +class DataBoxWrapper +{ +private: + DataBox databox_; + // Avoid re-calculating temperature bounds on every call + DataType_t Tlo_; + DataType_t Thi_; + +public: + SINGE_PORTABLE_FUNCTION DataBoxWrapper(DataBox databox) + : databox_{databox} // MUST be initialized before Tlo_ and Thi_ + , Tlo_{get_temperature(0)} + , Thi_{get_temperature(size() - 1)} + {} + SINGE_PORTABLE_FUNCTION constexpr DataType_t get_logT(int const index) const + { + return databox_.range(0).x(index); + } + SINGE_PORTABLE_FUNCTION constexpr DataType_t get_temperature(int const index) const + { + return std::exp(get_logT(index)); + } + SINGE_PORTABLE_FUNCTION constexpr DataType_t get_value(int const index) const + { + return databox_(index); + } + SINGE_PORTABLE_FUNCTION constexpr DataType_t get_logV(int const index) const + { + return std::log(get_value(index)); + } + SINGE_PORTABLE_FUNCTION constexpr int size() const + { + return databox_.dim(1); + } + SINGE_PORTABLE_FUNCTION constexpr DataType_t operator()(DataType_t const temperature) const + { + if (temperature < Tlo_) { + return LowerExtrapolation::extrapolate(temperature, *this); + } else if (temperature > Thi_) { + return UpperExtrapolation::extrapolate(temperature, *this); + } else { + return databox_.interpToReal(std::log(temperature)); + } + } + SINGE_PORTABLE_FUNCTION constexpr void finalize() + { + databox_.finalize(); + } +}; + +// ================================================================================================ + +} // end namespace util +} // end namespace singe + +#endif // ifndef SINGE_UTIL_DATABOXWRAPPER_HPP From ac86b2fc9e3fe3a09a938c3dc281ddb993fc77ac Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Tue, 3 Sep 2024 08:55:25 -0600 Subject: [PATCH 02/48] Notes --- spiner/databox_wrapper.hpp | 107 +++++++++++++++++++++++++++++++------ 1 file changed, 90 insertions(+), 17 deletions(-) diff --git a/spiner/databox_wrapper.hpp b/spiner/databox_wrapper.hpp index 375d1d622..aba41cce8 100644 --- a/spiner/databox_wrapper.hpp +++ b/spiner/databox_wrapper.hpp @@ -1,16 +1,90 @@ -#ifndef SINGE_UTIL_DATABOXWRAPPER_HPP -#define SINGE_UTIL_DATABOXWRAPPER_HPP +#ifndef _SPINER_DATABOXWRAPPER_HPP_ +#define _SPINER_DATABOXWRAPPER_HPP_ -#include "singe/util/databox.hpp" -#include "singe/util/portability_macros.hpp" +#include "databox.hpp" + +#include "ports-of-call//portability.hpp" #include -namespace singe { -namespace util { +namespace Spiner { // ================================================================================================ +// The original use-case did three things: +// * Add transformations of the data. This was originally hard-coded, but we'll make it general. +// * Add extrapolations off the edge of the table. +// * Adapt the interface. We should instead follow the Databox interface. +// * DataBox() = default; +// * DataBox(T * data, Args... args) noexcept; +// * DataBox(AllocationTarget t, Args... args) noexcept; +// * DataBox(Args... args) noexcept; +// * DataBox(PortableMDArray A) noexcept; +// * DataBox(PortableMDArray& A) noexcept; +// * DataBox(const DataBox& src) noexcept; +// * DataBox(const DataBox& b, const int dim, const int indx, const int nvar) noexcept +// #ifdef SPINER_USE_HDF +// * DataBox(const std::string& filename); +// * DataBox(hid_t loc, const std::string& groupname); +// #endif // SPINER_USE_HDF +// * void setArray(PortableMDArray& A); +// * void resize(AllocationTarget t, Args... args); +// * void resize(Args... args); +// * T& operator()(Args... args); +// * T& operator()(Args... args) const; +// * DataBox slice(const int dim, const int indx, const int nvar) const; +// * DataBox slice(const int indx) const; +// * DataBox slice(const int ix2, const int ix1) const; +// * void reshape(Args... args); +// * T interpToReal(const T x) const noexcept; +// * T interpToReal(const T x2, const T x1) const noexcept; +// * T interpToReal(const T x3, const T x2, const T x1) const noexcept; +// * T interpToReal(const T x4, const T x3, const T x2, const T x1) const noexcept; +// * void interpFromDB(const DataBox& db, const T x); +// * void interpFromDB(const DataBox& db, const T x2, const T x1); +// * DataBox interpToDB(Args... args); +// * void setIndexType(int i, IndexType t); +// * void setRange(int i, Grid_t g); +// * void setRange(int i, Args&&... args); +// * void copyShape(const DataBox& db, const int ndims = 0); +// * void copyMetadata(const DataBox& src); +// #ifdef SPINER_USE_HDF +// * herr_t saveHDF() const; +// * herr_t saveHDF(const std::string& filename) const; +// * herr_t saveHDF(hid_t loc, const std::string& groupname) const; +// * herr_t loadHDF(); +// * herr_t loadHDF(const std::string& filename); +// * herr_t loadHDF(hid_t loc, const std::string& groupname); +// #endif // SPINER_USE_HDF +// * IndexType& indexType(const int i); +// * Grid_t& range(const int i); +// * DataBox& operator=(const DataBox& other); +// * void copy(const DataBox& src); +// * DataStatus dataStatus() const; +// * bool isReference(); +// * bool ownsAllocatedMemory(); +// * bool operator==(const DataBox& other) const; +// * bool operator!=(const DataBox& other) const; +// * void makeShallow(); +// * void reset(); +// * T* data() const; +// * T min() const; +// * T max() const; +// * int rank() const; +// * int size() const; +// * int sizeBytes() const; +// * int dim() const; +// * Grid_t range(int i) const; +// * IndexType indexType(const int i) const; +// * std::size_t serializedSizeInBytes() const; +// * std::size_t serialize(char* dst) const; +// * std::size_t setPointer(T* src); +// * std::size_t setPointer(char* src); +// * std::size_t deSerialize(char* src); +// * DataBox getOnDevice() const; +// * void finalize(); + + // a wrapper to (a) adapt the interface and (b) handle extrapolation off the table // -- note that we assume this is a 1D databox template @@ -23,32 +97,32 @@ class DataBoxWrapper DataType_t Thi_; public: - SINGE_PORTABLE_FUNCTION DataBoxWrapper(DataBox databox) + PORTABLE_FUNCTION DataBoxWrapper(DataBox databox) : databox_{databox} // MUST be initialized before Tlo_ and Thi_ , Tlo_{get_temperature(0)} , Thi_{get_temperature(size() - 1)} {} - SINGE_PORTABLE_FUNCTION constexpr DataType_t get_logT(int const index) const + PORTABLE_FUNCTION constexpr DataType_t get_logT(int const index) const { return databox_.range(0).x(index); } - SINGE_PORTABLE_FUNCTION constexpr DataType_t get_temperature(int const index) const + PORTABLE_FUNCTION constexpr DataType_t get_temperature(int const index) const { return std::exp(get_logT(index)); } - SINGE_PORTABLE_FUNCTION constexpr DataType_t get_value(int const index) const + PORTABLE_FUNCTION constexpr DataType_t get_value(int const index) const { return databox_(index); } - SINGE_PORTABLE_FUNCTION constexpr DataType_t get_logV(int const index) const + PORTABLE_FUNCTION constexpr DataType_t get_logV(int const index) const { return std::log(get_value(index)); } - SINGE_PORTABLE_FUNCTION constexpr int size() const + PORTABLE_FUNCTION constexpr int size() const { return databox_.dim(1); } - SINGE_PORTABLE_FUNCTION constexpr DataType_t operator()(DataType_t const temperature) const + PORTABLE_FUNCTION constexpr DataType_t operator()(DataType_t const temperature) const { if (temperature < Tlo_) { return LowerExtrapolation::extrapolate(temperature, *this); @@ -58,7 +132,7 @@ class DataBoxWrapper return databox_.interpToReal(std::log(temperature)); } } - SINGE_PORTABLE_FUNCTION constexpr void finalize() + PORTABLE_FUNCTION constexpr void finalize() { databox_.finalize(); } @@ -66,7 +140,6 @@ class DataBoxWrapper // ================================================================================================ -} // end namespace util -} // end namespace singe +} // end namespace Spiner -#endif // ifndef SINGE_UTIL_DATABOXWRAPPER_HPP +#endif // ifndef _SPINER_DATABOXWRAPPER_HPP_ From 7d0865a9b6806d5ac2539b8b0b2d60d74603d10e Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Tue, 3 Sep 2024 10:26:42 -0600 Subject: [PATCH 03/48] Notes again --- spiner/databox_wrapper.hpp | 52 +++++++++++++++++++++++++++++++++----- 1 file changed, 45 insertions(+), 7 deletions(-) diff --git a/spiner/databox_wrapper.hpp b/spiner/databox_wrapper.hpp index aba41cce8..1d5e7c345 100644 --- a/spiner/databox_wrapper.hpp +++ b/spiner/databox_wrapper.hpp @@ -84,17 +84,55 @@ namespace Spiner { // * DataBox getOnDevice() const; // * void finalize(); +// TODO: +// * In theory the extra features could be split: +// * One wrapper that adds transformations +// * One wrapper that adds extrapolation +// * If you want both, you can have Extrapolation>> +// * But would it also work if you did Transformation>>? +// * This idea seems overly complicated for little benefit. -// a wrapper to (a) adapt the interface and (b) handle extrapolation off the table -// -- note that we assume this is a 1D databox -template +// TODO: Name? +// DataBoxPlusPlus +// ExtendedDataBox +// DataBoxButCooler +// FancyPantsDataBox +// DataBoxNowWithMoreFeatures +// DataBoxBetterStrongerFaster + +// TODO: Template defaults? +// * The obvious/naive defaults would be to match an unwrapped DataBox +// * All extrapolations default to "generate an error" +// * All transformations are the identity transformation +// * But given that the entire purpose of this wrapper is to add features beyond that of the basic +// DataBox, does it make sense for it to have defaults that make this a DataBox but with extra +// layers of complexity? +// * It's also not clear what order the template arguments should be in to manage defaults. +template < + typename DataType_t, + // TODO: Right now all axes will use the same lower and upper extrapolations, which may or may + // not be the desired behavior. + typename LowerExtrapolation, typename UpperExtrapolation, + typename TransformX, typename TransformY> class DataBoxWrapper { private: - DataBox databox_; - // Avoid re-calculating temperature bounds on every call - DataType_t Tlo_; - DataType_t Thi_; + // TODO: Deal with the other templates. + using DB_t = DataBox; + DB_t databox_; + // Avoid re-calculating independent variable bounds on every call + // TODO: My original implementation in Singe did this, but after thinking about it I think we + // should throw this out. + // * When interpolating, we already have to transform the independent variable(s), so no + // extra work. + // * When extrapolating, we don't use the transformed independent variable(s), so we do + // an extra transformation. + // * Finding the _exact_ bounds, that's going to be defined by the underlying DataBox, + // which is in transformed space. If we do bounds checking in the untransformed space, + // we have the possibility of numerical error creating a small space between the + // untransformed bound and the transformed bound, which could create errors. + DataType_t [DB_t::MAXRANK] lower_bounds_; + DataType_t [DB_t::MAXRANK] upper_bounds_; public: PORTABLE_FUNCTION DataBoxWrapper(DataBox databox) From f99c24148dbe5092ebb675033488746d7f148ea6 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Fri, 13 Sep 2024 16:11:42 -0600 Subject: [PATCH 04/48] Start adding the transformations. --- spiner/transformations.hpp | 54 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) create mode 100644 spiner/transformations.hpp diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp new file mode 100644 index 000000000..b0447678c --- /dev/null +++ b/spiner/transformations.hpp @@ -0,0 +1,54 @@ +#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 + +// TODO: Write tests, an ensure that transformations are symmetric (within some tolerance, possibly +// exact for critical values depending on the transformation). + +namespace Spiner { + + // linear transformation (aka no-op): y = x + struct TransformLinear{ + template + PORTABLE_INLINE_FUNCTION T forward(const T x) { + return x; + } + PORTABLE_INLINE_FUNCTION T reverse(const T x) { + return x; + } + }; + + // logarithmic transformation: y = log(x + small) + struct TransformLogarithmic{ + template + PORTABLE_INLINE_FUNCTION T forward(const T x) { + return std::log(x + std::numeric_limits::denorm_min()); + } + PORTABLE_INLINE_FUNCTION T reverse(const T x) { + return std::exp(x) - std::numeric_limits::denorm_min(); + } + }; + + // TODO: log_NQT + + // TODO: arcsinh_NQT + +} // namespace Spiner +#endif // _SPINER_TRANSFORM_HPP_ From e0e8b0a7f6e29281b965a4d24b1cd25032a39b3d Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Fri, 13 Sep 2024 16:54:14 -0600 Subject: [PATCH 05/48] Started writing the transformations. --- spiner/transformations.hpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index b0447678c..58f536859 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -24,14 +24,20 @@ 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 PORTABLE_INLINE_FUNCTION T forward(const T x) { return x; } - PORTABLE_INLINE_FUNCTION T reverse(const T x) { - return x; + PORTABLE_INLINE_FUNCTION T reverse(const T u) { + return u; } }; @@ -41,8 +47,8 @@ namespace Spiner { PORTABLE_INLINE_FUNCTION T forward(const T x) { return std::log(x + std::numeric_limits::denorm_min()); } - PORTABLE_INLINE_FUNCTION T reverse(const T x) { - return std::exp(x) - std::numeric_limits::denorm_min(); + PORTABLE_INLINE_FUNCTION T reverse(const T u) { + return std::exp(u) - std::numeric_limits::denorm_min(); } }; From 5ce03cbdd62b5c3469ed4232052be438c8aa55de Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Fri, 13 Sep 2024 16:54:32 -0600 Subject: [PATCH 06/48] 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_; }; From be1acd1e9e3ae3168346d63bf6a75eb3ac5f01a1 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Fri, 13 Sep 2024 16:56:36 -0600 Subject: [PATCH 07/48] Forgot include --- spiner/regular_grid_1d.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index 6753d8fbb..b76a404b3 100644 --- a/spiner/regular_grid_1d.hpp +++ b/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 { From fc0f83f02f29d988aa00c8995ba6c1dff7f9baac Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 11:08:44 -0600 Subject: [PATCH 08/48] Forgot (a) some templates and (b) static. --- spiner/transformations.hpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index 58f536859..be548db62 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -33,10 +33,11 @@ namespace Spiner { // linear transformation (aka no-op): y = x struct TransformLinear{ template - PORTABLE_INLINE_FUNCTION T forward(const T x) { + PORTABLE_INLINE_FUNCTION static T forward(const T x) { return x; } - PORTABLE_INLINE_FUNCTION T reverse(const T u) { + template + PORTABLE_INLINE_FUNCTION static T reverse(const T u) { return u; } }; @@ -44,10 +45,11 @@ namespace Spiner { // logarithmic transformation: y = log(x + small) struct TransformLogarithmic{ template - PORTABLE_INLINE_FUNCTION T forward(const T x) { + PORTABLE_INLINE_FUNCTION static T forward(const T x) { return std::log(x + std::numeric_limits::denorm_min()); } - PORTABLE_INLINE_FUNCTION T reverse(const T u) { + template + PORTABLE_INLINE_FUNCTION static T reverse(const T u) { return std::exp(u) - std::numeric_limits::denorm_min(); } }; From 622c8fb4d035c1516edde1422f38abae8281f126 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 11:09:20 -0600 Subject: [PATCH 09/48] I removed dx() from RegularGrid1D, so update the tests. --- test/benchmark.cpp | 5 ++++- test/test.cpp | 1 - 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/test/benchmark.cpp b/test/benchmark.cpp index 07a7195dd..e59bcbfb8 100644 --- a/test/benchmark.cpp +++ b/test/benchmark.cpp @@ -74,8 +74,10 @@ int main(int argc, char *argv[]) { RegularGrid1D gcoarse(xmin, xmax, ncoarse); std::vector gfine; + std::vector 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(); diff --git a/test/test.cpp b/test/test.cpp index 30671ff93..453c026e9 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -98,7 +98,6 @@ TEST_CASE("RegularGrid1D", "[RegularGrid1D]") { REQUIRE(g.min() == min); REQUIRE(g.max() == max); REQUIRE(g.nPoints() == N); - REQUIRE(g.dx() == (max - min) / ((Real)(N - 1))); } } From 2eccab5d9fc37ab7f73f770d50cd78a7ad9a9f90 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 11:12:28 -0600 Subject: [PATCH 10/48] Tests for transformations. --- test/CMakeLists.txt | 2 +- test/transformations.cpp | 60 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 61 insertions(+), 1 deletion(-) create mode 100644 test/transformations.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 8e019ce3e..bfb7fa01e 100644 --- a/test/CMakeLists.txt +++ b/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 transformations.cpp) target_link_libraries(test.bin PRIVATE spiner::spiner diff --git a/test/transformations.cpp b/test/transformations.cpp new file mode 100644 index 000000000..638452394 --- /dev/null +++ b/test/transformations.cpp @@ -0,0 +1,60 @@ +// © (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 + +#include +#include +#include + +#include + +using Catch::Matchers::WithinAbs; +using Catch::Matchers::WithinRel; +using Catch::Matchers::WithinULP; + +TEST_CASE("transformation: linear", "") { + // This one is almost too simple to meaningfully test, but we can at least + // ensure that it compiles and does the trivially-obvious things. + using Transform = Spiner::TransformLinear; + const double x0 = 3.14159; + CHECK(Transform::forward(x0) == x0); + CHECK(Transform::reverse(x0) == x0); + const float x1 = 2.71828; + CHECK(Transform::forward(x1) == x1); + CHECK(Transform::reverse(x1) == x1); + // TODO: Test on GPU +} + +TEST_CASE("transformation: logarithmic", "") { + using Transform = Spiner::TransformLogarithmic; + auto matcher = [](const double d) { + //return WithinULP(d, 500); + return WithinRel(d, 6.0e-12); + }; + // Scan across most of the range + for (int p = -300; p <= 300; ++p) { + const double x = std::pow(double(10), p); + const double y = std::log(x); + // Basic checks + CHECK_THAT(Transform::forward(x), matcher(y)); + CHECK_THAT(Transform::reverse(y), matcher(x)); + // Round trip + CHECK_THAT(Transform::reverse(Transform::forward(x)), matcher(x)); + } + // Special value + CHECK(std::isfinite(Transform::forward(0))); + CHECK(Transform::reverse(Transform::forward(0)) == 0); + // TODO: Test on GPU +} + From becaae9ee51c749b80f71ac9fff4b5bbacca9df8 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 11:25:55 -0600 Subject: [PATCH 11/48] Add Transform template to DataBox (currently doesn't do anything yet). --- spiner/databox.hpp | 151 +++++++++++++++++++++++---------------------- 1 file changed, 76 insertions(+), 75 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index 793664364..0b840bde5 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -51,6 +51,7 @@ enum class DataStatus { enum class AllocationTarget { Host, Device }; template , + typename Transform = TransformLinear, typename Concept = typename std::enable_if::value, bool>::type> class DataBox { @@ -101,7 +102,7 @@ class DataBox { setAllIndexed_(); } PORTABLE_INLINE_FUNCTION - DataBox(const DataBox &src) noexcept + DataBox(const DataBox &src) noexcept : rank_(src.rank_), status_(src.status_), data_(src.data_) { setAllIndexed_(); dataView_.InitWithShallowSlice(src.dataView_, 6, 0, src.dim(6)); @@ -113,7 +114,7 @@ class DataBox { // Slice constructor PORTABLE_INLINE_FUNCTION - DataBox(const DataBox &b, const int dim, const int indx, + DataBox(const DataBox &b, const int dim, const int indx, const int nvar) noexcept : status_(DataStatus::Unmanaged), data_(b.data_) { dataView_.InitWithShallowSlice(b.dataView_, dim, indx, nvar); @@ -164,15 +165,15 @@ class DataBox { // Slice operation PORTABLE_INLINE_FUNCTION - DataBox slice(const int dim, const int indx, + DataBox slice(const int dim, const int indx, const int nvar) const { return DataBox(*this, dim, indx, nvar); } - PORTABLE_INLINE_FUNCTION DataBox + PORTABLE_INLINE_FUNCTION DataBox slice(const int indx) const { return slice(rank_, indx, 1); } - PORTABLE_INLINE_FUNCTION DataBox + PORTABLE_INLINE_FUNCTION DataBox slice(const int ix2, const int ix1) const { // DataBox a(*this, rank_, ix2, 1); // return DataBox(a, a.rank_, ix1, 1); @@ -213,13 +214,13 @@ class DataBox { // WARNING: requires memory to be pre-allocated. // TODO: add 3d and higher interpFromDB if necessary PORTABLE_INLINE_FUNCTION void - interpFromDB(const DataBox &db, const T x); + interpFromDB(const DataBox &db, const T x); PORTABLE_INLINE_FUNCTION void - interpFromDB(const DataBox &db, const T x2, const T x1); + interpFromDB(const DataBox &db, const T x2, const T x1); template - PORTABLE_INLINE_FUNCTION DataBox + PORTABLE_INLINE_FUNCTION DataBox interpToDB(Args... args) { - DataBox db; + DataBox db; db.interpFromDB(*this, std::forward(args)...); return db; } @@ -245,10 +246,10 @@ class DataBox { // Does no checks that memory is available. // Optionally copies shape of source with ndims fewer slowest-moving // dimensions - PORTABLE_INLINE_FUNCTION void copyShape(const DataBox &db, + PORTABLE_INLINE_FUNCTION void copyShape(const DataBox &db, const int ndims = 0); // Returns new databox with same memory and metadata - inline void copyMetadata(const DataBox &src); + inline void copyMetadata(const DataBox &src); #ifdef SPINER_USE_HDF inline herr_t saveHDF() const { return saveHDF(SP5::DB::FILENAME); } @@ -264,9 +265,9 @@ class DataBox { inline Grid_t &range(const int i) { return grids_[i]; } // Assignment and move, both perform shallow copy - PORTABLE_INLINE_FUNCTION DataBox & - operator=(const DataBox &other); - inline void copy(const DataBox &src); + PORTABLE_INLINE_FUNCTION DataBox & + operator=(const DataBox &other); + inline void copy(const DataBox &src); // utility info inline DataStatus dataStatus() const { return status_; } @@ -274,8 +275,8 @@ class DataBox { inline bool ownsAllocatedMemory() { return (status_ != DataStatus::Unmanaged); } - inline bool operator==(const DataBox &other) const; - inline bool operator!=(const DataBox &other) const { + inline bool operator==(const DataBox &other) const; + inline bool operator!=(const DataBox &other) const { return !(*this == other); } @@ -370,11 +371,11 @@ class DataBox { } // ------------------------------------ - DataBox + DataBox getOnDevice() const { // getOnDevice is always a deep copy if (size() == 0 || status_ == DataStatus::Empty) { // edge case for unallocated - DataBox a; + DataBox a; return a; } // create device memory (host memory if no device) @@ -382,7 +383,7 @@ class DataBox { // copy to device portableCopyToDevice(device_data, data_, sizeBytes()); // create new databox of size size - DataBox a{device_data, dim(6), dim(5), dim(4), + DataBox a{device_data, dim(6), dim(5), dim(4), dim(3), dim(2), dim(1)}; a.copyShape(*this); // set correct allocation status of the new databox @@ -439,16 +440,16 @@ class DataBox { }; // Read an array, shallow -template -inline void DataBox::setArray(PortableMDArray &A) { +template +inline void DataBox::setArray(PortableMDArray &A) { dataView_ = A; rank_ = A.GetRank(); setAllIndexed_(); } -template +template PORTABLE_INLINE_FUNCTION T -DataBox::interpToReal(const T x) const noexcept { +DataBox::interpToReal(const T x) const noexcept { assert(canInterpToReal_(1)); int ix; weights_t w; @@ -456,8 +457,8 @@ DataBox::interpToReal(const T x) const noexcept { return w[0] * dataView_(ix) + w[1] * dataView_(ix + 1); } -template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +template +PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( const T x2, const T x1) const noexcept { assert(canInterpToReal_(2)); int ix1, ix2; @@ -472,8 +473,8 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( w1[1] * dataView_(ix2 + 1, ix1 + 1))); } -template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +template +PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( const T x3, const T x2, const T x1) const noexcept { assert(canInterpToReal_(3)); int ix[3]; @@ -495,8 +496,8 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( w[0][1] * dataView_(ix[2] + 1, ix[1] + 1, ix[0] + 1)))); } -template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +template +PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( const T x3, const T x2, const T x1, const int idx) const noexcept { assert(rank_ == 4); for (int r = 1; r < rank_; ++r) { @@ -526,8 +527,8 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( // DH: this is a large function to force an inline, perhaps just make it a // suggestion to the compiler? -template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +template +PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( const T x4, const T x3, const T x2, const T x1) const noexcept { assert(canInterpToReal_(4)); T x[] = {x1, x2, x3, x4}; @@ -577,8 +578,8 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( ); } -template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +template +PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( const T x4, const T x3, const T x2, const int idx, const T x1) const noexcept { assert(rank_ == 5); @@ -640,9 +641,9 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( ); } -template +template PORTABLE_INLINE_FUNCTION void -DataBox::interpFromDB(const DataBox &db, +DataBox::interpFromDB(const DataBox &db, const T x) { assert(db.indices_[db.rank_ - 1] == IndexType::Interpolated); assert(db.grids_[db.rank_ - 1].isWellFormed()); @@ -653,7 +654,7 @@ DataBox::interpFromDB(const DataBox &db, copyShape(db, 1); db.grids_[db.rank_ - 1].weights(x, ix, w); - DataBox lower(db.slice(ix)), upper(db.slice(ix + 1)); + DataBox lower(db.slice(ix)), upper(db.slice(ix + 1)); // lower = db.slice(ix); // upper = db.slice(ix+1); for (int i = 0; i < size(); i++) { @@ -661,9 +662,9 @@ DataBox::interpFromDB(const DataBox &db, } } -template +template PORTABLE_INLINE_FUNCTION void -DataBox::interpFromDB(const DataBox &db, +DataBox::interpFromDB(const DataBox &db, const T x2, const T x1) { assert(db.rank_ >= 2); assert(db.indices_[db.rank_ - 1] == IndexType::Interpolated); @@ -678,7 +679,7 @@ DataBox::interpFromDB(const DataBox &db, db.grids_[db.rank_ - 2].weights(x1, ix1, w1); db.grids_[db.rank_ - 1].weights(x2, ix2, w2); - DataBox corners[2][2]{ + DataBox corners[2][2]{ {db.slice(ix2, ix1), db.slice(ix2 + 1, ix1)}, {db.slice(ix2, ix1 + 1), db.slice(ix2 + 1, ix1 + 1)}}; // copyShape(db,2); @@ -707,9 +708,9 @@ DataBox::interpFromDB(const DataBox &db, // Reshapes from other databox, but does not allocate memory. // Does no checks that memory is available. // Optionally copies shape of source with ndims fewer slowest-moving dimensions -template +template PORTABLE_INLINE_FUNCTION void -DataBox::copyShape(const DataBox &db, +DataBox::copyShape(const DataBox &db, const int ndims) { rank_ = db.rank_ - ndims; int dims[MAXRANK]; @@ -727,9 +728,9 @@ DataBox::copyShape(const DataBox &db, } // reallocates and then copies shape from other databox // everything but the actual copy in a deep copy -template -inline void DataBox::copyMetadata( - const DataBox &src) { +template +inline void DataBox::copyMetadata( + const DataBox &src) { AllocationTarget t = (src.status_ == DataStatus::AllocatedDevice ? AllocationTarget::Device : AllocationTarget::Host); @@ -743,9 +744,9 @@ inline void DataBox::copyMetadata( } #ifdef SPINER_USE_HDF -template +template inline herr_t -DataBox::saveHDF(const std::string &filename) const { +DataBox::saveHDF(const std::string &filename) const { herr_t status; hid_t file; @@ -755,9 +756,9 @@ DataBox::saveHDF(const std::string &filename) const { return status; } -template +template inline herr_t -DataBox::saveHDF(hid_t loc, +DataBox::saveHDF(hid_t loc, const std::string &groupname) const { hid_t group, grids; herr_t status = 0; @@ -817,18 +818,18 @@ DataBox::saveHDF(hid_t loc, return status; } -template +template inline herr_t -DataBox::loadHDF(const std::string &filename) { +DataBox::loadHDF(const std::string &filename) { herr_t status; hid_t file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT); status = loadHDF(file, SP5::DB::GRPNAME); return status; } -template +template inline herr_t -DataBox::loadHDF(hid_t loc, const std::string &groupname) { +DataBox::loadHDF(hid_t loc, const std::string &groupname) { hid_t group, grids; herr_t status = 0; std::vector index_types; @@ -883,9 +884,9 @@ DataBox::loadHDF(hid_t loc, const std::string &groupname) { #endif // SPINER_USE_HDF // Performs shallow copy by default -template -PORTABLE_INLINE_FUNCTION DataBox & -DataBox::operator=(const DataBox &src) { +template +PORTABLE_INLINE_FUNCTION DataBox & +DataBox::operator=(const DataBox &src) { if (this != &src) { rank_ = src.rank_; status_ = src.status_; @@ -900,17 +901,17 @@ DataBox::operator=(const DataBox &src) { } // Performs a deep copy -template +template inline void -DataBox::copy(const DataBox &src) { +DataBox::copy(const DataBox &src) { copyMetadata(src); for (int i = 0; i < src.size(); i++) dataView_(i) = src(i); } -template -inline bool DataBox::operator==( - const DataBox &other) const { +template +inline bool DataBox::operator==( + const DataBox &other) const { if (rank_ != other.rank_) return false; for (int i = 0; i < rank_; i++) { if (indices_[i] != other.indices_[i]) return false; @@ -923,8 +924,8 @@ inline bool DataBox::operator==( } // TODO: should this be std::reduce? -template -PORTABLE_INLINE_FUNCTION T DataBox::min() const { +template +PORTABLE_INLINE_FUNCTION T DataBox::min() const { T min = std::numeric_limits::infinity(); for (int i = 0; i < size(); i++) { min = std::min(min, dataView_(i)); @@ -932,8 +933,8 @@ PORTABLE_INLINE_FUNCTION T DataBox::min() const { return min; } -template -PORTABLE_INLINE_FUNCTION T DataBox::max() const { +template +PORTABLE_INLINE_FUNCTION T DataBox::max() const { T max = -std::numeric_limits::infinity(); for (int i = 0; i < size(); i++) { max = std::max(max, dataView_(i)); @@ -941,24 +942,24 @@ PORTABLE_INLINE_FUNCTION T DataBox::max() const { return max; } -template +template PORTABLE_INLINE_FUNCTION Grid_t -DataBox::range(int i) const { +DataBox::range(int i) const { assert(0 <= i && i < rank_); assert(indices_[i] == IndexType::Interpolated); return grids_[i]; } -template -PORTABLE_INLINE_FUNCTION void DataBox::setAllIndexed_() { +template +PORTABLE_INLINE_FUNCTION void DataBox::setAllIndexed_() { for (int i = 0; i < rank_; i++) { indices_[i] = IndexType::Indexed; } } -template +template PORTABLE_INLINE_FUNCTION bool -DataBox::canInterpToReal_(const int interpOrder) const { +DataBox::canInterpToReal_(const int interpOrder) const { if (rank_ != interpOrder) return false; for (int i = 0; i < rank_; i++) { if (indices_[i] != IndexType::Interpolated) return false; @@ -967,13 +968,13 @@ DataBox::canInterpToReal_(const int interpOrder) const { return true; } -template -inline DataBox -getOnDeviceDataBox(const DataBox &a_host) { +template +inline DataBox +getOnDeviceDataBox(const DataBox &a_host) { return a_host.getOnDevice(); } -template -inline void free(DataBox &db) { +template +inline void free(DataBox &db) { db.finalize(); } From 18cd066a5502e9e7ff8d53dc891a1e72371125a4 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 11:49:39 -0600 Subject: [PATCH 12/48] Update TODO items --- spiner/regular_grid_1d.hpp | 2 ++ spiner/transformations.hpp | 3 --- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index b76a404b3..32e0bae9d 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -45,6 +45,8 @@ struct weights_t { } }; +// TODO: Need to expand testing + // 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). diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index be548db62..dfac62466 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -19,9 +19,6 @@ #include -// TODO: Write tests, an ensure that transformations are symmetric (within some tolerance, possibly -// exact for critical values depending on the transformation). - namespace Spiner { // Note on notation: From af535496a883deac67bec0f3d3aed1bf51bdaf42 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 11:50:19 -0600 Subject: [PATCH 13/48] Add dependent variable transformation. --- spiner/databox.hpp | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index 0b840bde5..fb29f4de9 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -454,7 +454,8 @@ DataBox::interpToReal(const T x) const noexcept { int ix; weights_t w; grids_[0].weights(x, ix, w); - return w[0] * dataView_(ix) + w[1] * dataView_(ix + 1); + const auto v = w[0] * dataView_(ix) + w[1] * dataView_(ix + 1); + return Transform::reverse(v); } template @@ -467,10 +468,11 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo grids_[1].weights(x2, ix2, w2); // TODO: prefectch corners for speed? // TODO: re-order access pattern? - return (w2[0] * + const auto v = (w2[0] * (w1[0] * dataView_(ix2, ix1) + w1[1] * dataView_(ix2, ix1 + 1)) + w2[1] * (w1[0] * dataView_(ix2 + 1, ix1) + w1[1] * dataView_(ix2 + 1, ix1 + 1))); + return Transform::reverse(v); } template @@ -484,7 +486,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo grids_[2].weights(x3, ix[2], w[2]); // TODO: prefect corners for speed? // TODO: re-order access pattern? - return ( + const auto v = ( w[2][0] * (w[1][0] * (w[0][0] * dataView_(ix[2], ix[1], ix[0]) + w[0][1] * dataView_(ix[2], ix[1], ix[0] + 1)) + w[1][1] * (w[0][0] * dataView_(ix[2], ix[1] + 1, ix[0]) + @@ -494,6 +496,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo w[0][1] * dataView_(ix[2] + 1, ix[1], ix[0] + 1)) + w[1][1] * (w[0][0] * dataView_(ix[2] + 1, ix[1] + 1, ix[0]) + w[0][1] * dataView_(ix[2] + 1, ix[1] + 1, ix[0] + 1)))); + return Transform::reverse(v); } template @@ -511,7 +514,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo grids_[3].weights(x3, ix[2], w[2]); // TODO: prefect corners for speed? // TODO: re-order access pattern? - return ( + const auto v = ( w[2][0] * (w[1][0] * (w[0][0] * dataView_(ix[2], ix[1], ix[0], idx) + w[0][1] * dataView_(ix[2], ix[1], ix[0] + 1, idx)) + @@ -523,6 +526,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo w[1][1] * (w[0][0] * dataView_(ix[2] + 1, ix[1] + 1, ix[0], idx) + w[0][1] * dataView_(ix[2] + 1, ix[1] + 1, ix[0] + 1, idx)))); + return Transform::reverse(v); } // DH: this is a large function to force an inline, perhaps just make it a @@ -540,7 +544,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo // TODO(JMM): This is getty pretty gross. Should we automate? // Hand-written is probably faster, though. // Breaking line-limit to make this easier to read - return ( + const auto v = ( w[3][0] * (w[2][0] * (w[1][0] * @@ -576,6 +580,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo ix[1] + 1, ix[0] + 1)))) ); + return Transform::reverse(v); } template @@ -599,7 +604,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo // TODO(JMM): This is getty pretty gross. Should we automate? // Hand-written is probably faster, though. // Breaking line-limit to make this easier to read - return ( + const auto v = ( w[3][0] * (w[2][0] * (w[1][0] * @@ -639,6 +644,7 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo idx, ix[0] + 1)))) ); + return Transform::reverse(v); } template @@ -928,7 +934,7 @@ template PORTABLE_INLINE_FUNCTION T DataBox::min() const { T min = std::numeric_limits::infinity(); for (int i = 0; i < size(); i++) { - min = std::min(min, dataView_(i)); + min = std::min(min, Transform::reverse(dataView_(i))); } return min; } @@ -937,7 +943,7 @@ template PORTABLE_INLINE_FUNCTION T DataBox::max() const { T max = -std::numeric_limits::infinity(); for (int i = 0; i < size(); i++) { - max = std::max(max, dataView_(i)); + max = std::max(max, Transform::reverse(dataView_(i))); } return max; } From 4f8fbe2106a510d801980c0934536aa416499d94 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 11:50:59 -0600 Subject: [PATCH 14/48] Update TODO items --- spiner/databox.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index fb29f4de9..c857b0063 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -39,6 +39,8 @@ // TODO: get named indices working // TODO: If asserts are too slow, remove them. +// TODO: Need to expand testing to cover transformations (both independent and dependent). + namespace Spiner { enum class IndexType { Interpolated = 0, Named = 1, Indexed = 2 }; From 6e292a00b995fbbdda4bc47a9df3e9327d68bfdf Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 11:53:02 -0600 Subject: [PATCH 15/48] Update TODO items --- spiner/databox.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index c857b0063..b84ade6dd 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -41,6 +41,8 @@ // TODO: Need to expand testing to cover transformations (both independent and dependent). +// TODO: Should extrapolations be this MR or a separate MR? + namespace Spiner { enum class IndexType { Interpolated = 0, Named = 1, Indexed = 2 }; From fa2dd5b8925324dfa63cbbd305eb3a087b225892 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 14:35:26 -0600 Subject: [PATCH 16/48] Add GPU testing for transformations. --- test/transformations.cpp | 49 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 47 insertions(+), 2 deletions(-) diff --git a/test/transformations.cpp b/test/transformations.cpp index 638452394..3e553cca3 100644 --- a/test/transformations.cpp +++ b/test/transformations.cpp @@ -23,21 +23,50 @@ using Catch::Matchers::WithinAbs; using Catch::Matchers::WithinRel; using Catch::Matchers::WithinULP; +namespace { + PORTABLE_FUNCTION double relerr(const double z, const double zz) { + if (z == 0) { + return std::abs(zz - z); + } else { + return std::abs(double(1) - zz/z); + } + } +}; + TEST_CASE("transformation: linear", "") { // This one is almost too simple to meaningfully test, but we can at least // ensure that it compiles and does the trivially-obvious things. using Transform = Spiner::TransformLinear; + + // Test on CPU const double x0 = 3.14159; CHECK(Transform::forward(x0) == x0); CHECK(Transform::reverse(x0) == x0); const float x1 = 2.71828; CHECK(Transform::forward(x1) == x1); CHECK(Transform::reverse(x1) == x1); - // TODO: Test on GPU + + // Test on GPU (or CPU if no GPU available) + const int num_threads = 10; + double accum = 0; + portableReduce( + "parallel_reduce", 0, num_threads, + PORTABLE_LAMBDA(const int n, double &reduce) + { + const double x = static_cast(n); + const double y = x; + const double yy = Transform::forward(x); + const double xx = Transform::reverse(y); + reduce += 0.5 * (relerr(x, xx) + relerr(y, yy)); + }, + accum); + CHECK(accum / num_threads == 0); } TEST_CASE("transformation: logarithmic", "") { using Transform = Spiner::TransformLogarithmic; + + // Test on CPU auto matcher = [](const double d) { //return WithinULP(d, 500); return WithinRel(d, 6.0e-12); @@ -55,6 +84,22 @@ TEST_CASE("transformation: logarithmic", "") { // Special value CHECK(std::isfinite(Transform::forward(0))); CHECK(Transform::reverse(Transform::forward(0)) == 0); - // TODO: Test on GPU + + // Test on GPU (or CPU if no GPU available) + const int num_threads = 101; + double accum = 0; + portableReduce( + "parallel_reduce", 0, num_threads, + PORTABLE_LAMBDA(const int n, double &reduce) + { + const int p = n - num_threads / 2; + const double x = std::pow(double(10.0), p); + const double y = std::log(x); + const double yy = Transform::forward(x); + const double xx = Transform::reverse(y); + reduce += 0.5 * (relerr(x, xx) + relerr(y, yy)); + }, + accum); + CHECK(accum / num_threads <= 6.0e-12); } From c54a6a5a85a48510bb1dee0a1577324e9c9cb790 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 14:46:39 -0600 Subject: [PATCH 17/48] Delete a bunch of stuff we no longer need. --- spiner/databox_wrapper.hpp | 178 ------------------------------------- 1 file changed, 178 deletions(-) diff --git a/spiner/databox_wrapper.hpp b/spiner/databox_wrapper.hpp index 1d5e7c345..af0eb1319 100644 --- a/spiner/databox_wrapper.hpp +++ b/spiner/databox_wrapper.hpp @@ -1,183 +1,5 @@ -#ifndef _SPINER_DATABOXWRAPPER_HPP_ -#define _SPINER_DATABOXWRAPPER_HPP_ - -#include "databox.hpp" - -#include "ports-of-call//portability.hpp" - -#include - -namespace Spiner { - -// ================================================================================================ - // The original use-case did three things: // * Add transformations of the data. This was originally hard-coded, but we'll make it general. // * Add extrapolations off the edge of the table. // * Adapt the interface. We should instead follow the Databox interface. -// * DataBox() = default; -// * DataBox(T * data, Args... args) noexcept; -// * DataBox(AllocationTarget t, Args... args) noexcept; -// * DataBox(Args... args) noexcept; -// * DataBox(PortableMDArray A) noexcept; -// * DataBox(PortableMDArray& A) noexcept; -// * DataBox(const DataBox& src) noexcept; -// * DataBox(const DataBox& b, const int dim, const int indx, const int nvar) noexcept -// #ifdef SPINER_USE_HDF -// * DataBox(const std::string& filename); -// * DataBox(hid_t loc, const std::string& groupname); -// #endif // SPINER_USE_HDF -// * void setArray(PortableMDArray& A); -// * void resize(AllocationTarget t, Args... args); -// * void resize(Args... args); -// * T& operator()(Args... args); -// * T& operator()(Args... args) const; -// * DataBox slice(const int dim, const int indx, const int nvar) const; -// * DataBox slice(const int indx) const; -// * DataBox slice(const int ix2, const int ix1) const; -// * void reshape(Args... args); -// * T interpToReal(const T x) const noexcept; -// * T interpToReal(const T x2, const T x1) const noexcept; -// * T interpToReal(const T x3, const T x2, const T x1) const noexcept; -// * T interpToReal(const T x4, const T x3, const T x2, const T x1) const noexcept; -// * void interpFromDB(const DataBox& db, const T x); -// * void interpFromDB(const DataBox& db, const T x2, const T x1); -// * DataBox interpToDB(Args... args); -// * void setIndexType(int i, IndexType t); -// * void setRange(int i, Grid_t g); -// * void setRange(int i, Args&&... args); -// * void copyShape(const DataBox& db, const int ndims = 0); -// * void copyMetadata(const DataBox& src); -// #ifdef SPINER_USE_HDF -// * herr_t saveHDF() const; -// * herr_t saveHDF(const std::string& filename) const; -// * herr_t saveHDF(hid_t loc, const std::string& groupname) const; -// * herr_t loadHDF(); -// * herr_t loadHDF(const std::string& filename); -// * herr_t loadHDF(hid_t loc, const std::string& groupname); -// #endif // SPINER_USE_HDF -// * IndexType& indexType(const int i); -// * Grid_t& range(const int i); -// * DataBox& operator=(const DataBox& other); -// * void copy(const DataBox& src); -// * DataStatus dataStatus() const; -// * bool isReference(); -// * bool ownsAllocatedMemory(); -// * bool operator==(const DataBox& other) const; -// * bool operator!=(const DataBox& other) const; -// * void makeShallow(); -// * void reset(); -// * T* data() const; -// * T min() const; -// * T max() const; -// * int rank() const; -// * int size() const; -// * int sizeBytes() const; -// * int dim() const; -// * Grid_t range(int i) const; -// * IndexType indexType(const int i) const; -// * std::size_t serializedSizeInBytes() const; -// * std::size_t serialize(char* dst) const; -// * std::size_t setPointer(T* src); -// * std::size_t setPointer(char* src); -// * std::size_t deSerialize(char* src); -// * DataBox getOnDevice() const; -// * void finalize(); - -// TODO: -// * In theory the extra features could be split: -// * One wrapper that adds transformations -// * One wrapper that adds extrapolation -// * If you want both, you can have Extrapolation>> -// * But would it also work if you did Transformation>>? -// * This idea seems overly complicated for little benefit. - -// TODO: Name? -// DataBoxPlusPlus -// ExtendedDataBox -// DataBoxButCooler -// FancyPantsDataBox -// DataBoxNowWithMoreFeatures -// DataBoxBetterStrongerFaster - -// TODO: Template defaults? -// * The obvious/naive defaults would be to match an unwrapped DataBox -// * All extrapolations default to "generate an error" -// * All transformations are the identity transformation -// * But given that the entire purpose of this wrapper is to add features beyond that of the basic -// DataBox, does it make sense for it to have defaults that make this a DataBox but with extra -// layers of complexity? -// * It's also not clear what order the template arguments should be in to manage defaults. -template < - typename DataType_t, - // TODO: Right now all axes will use the same lower and upper extrapolations, which may or may - // not be the desired behavior. - typename LowerExtrapolation, typename UpperExtrapolation, - typename TransformX, typename TransformY> -class DataBoxWrapper -{ -private: - // TODO: Deal with the other templates. - using DB_t = DataBox; - DB_t databox_; - // Avoid re-calculating independent variable bounds on every call - // TODO: My original implementation in Singe did this, but after thinking about it I think we - // should throw this out. - // * When interpolating, we already have to transform the independent variable(s), so no - // extra work. - // * When extrapolating, we don't use the transformed independent variable(s), so we do - // an extra transformation. - // * Finding the _exact_ bounds, that's going to be defined by the underlying DataBox, - // which is in transformed space. If we do bounds checking in the untransformed space, - // we have the possibility of numerical error creating a small space between the - // untransformed bound and the transformed bound, which could create errors. - DataType_t [DB_t::MAXRANK] lower_bounds_; - DataType_t [DB_t::MAXRANK] upper_bounds_; - -public: - PORTABLE_FUNCTION DataBoxWrapper(DataBox databox) - : databox_{databox} // MUST be initialized before Tlo_ and Thi_ - , Tlo_{get_temperature(0)} - , Thi_{get_temperature(size() - 1)} - {} - PORTABLE_FUNCTION constexpr DataType_t get_logT(int const index) const - { - return databox_.range(0).x(index); - } - PORTABLE_FUNCTION constexpr DataType_t get_temperature(int const index) const - { - return std::exp(get_logT(index)); - } - PORTABLE_FUNCTION constexpr DataType_t get_value(int const index) const - { - return databox_(index); - } - PORTABLE_FUNCTION constexpr DataType_t get_logV(int const index) const - { - return std::log(get_value(index)); - } - PORTABLE_FUNCTION constexpr int size() const - { - return databox_.dim(1); - } - PORTABLE_FUNCTION constexpr DataType_t operator()(DataType_t const temperature) const - { - if (temperature < Tlo_) { - return LowerExtrapolation::extrapolate(temperature, *this); - } else if (temperature > Thi_) { - return UpperExtrapolation::extrapolate(temperature, *this); - } else { - return databox_.interpToReal(std::log(temperature)); - } - } - PORTABLE_FUNCTION constexpr void finalize() - { - databox_.finalize(); - } -}; - -// ================================================================================================ - -} // end namespace Spiner -#endif // ifndef _SPINER_DATABOXWRAPPER_HPP_ From 15aa48ae4a00109ee571e0110f8883511a40a0de Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 14:53:18 -0600 Subject: [PATCH 18/48] Move existing RegularGrid1D test to a new file, so I can add new testing of the new feature without further adding to the overcrowdedness of test.cpp. --- test/CMakeLists.txt | 2 +- test/regular_grid_1d.cpp | 31 +++++++++++++++++++++++++++++++ test/test.cpp | 12 ------------ 3 files changed, 32 insertions(+), 13 deletions(-) create mode 100644 test/regular_grid_1d.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index bfb7fa01e..d843689d3 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -57,7 +57,7 @@ target_link_libraries(spiner_test ${spinerTest_POPULATED_TARGETS} ) # add unit tests -add_executable(test.bin test.cpp transformations.cpp) +add_executable(test.bin test.cpp regular_grid_1d.cpp transformations.cpp) target_link_libraries(test.bin PRIVATE spiner::spiner diff --git a/test/regular_grid_1d.cpp b/test/regular_grid_1d.cpp new file mode 100644 index 000000000..51dd78ac1 --- /dev/null +++ b/test/regular_grid_1d.cpp @@ -0,0 +1,31 @@ +// © (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 +#include + +#include +#include + +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 g(min, max, N); + REQUIRE(g.min() == min); + REQUIRE(g.max() == max); + REQUIRE(g.nPoints() == N); + } +} + diff --git a/test/test.cpp b/test/test.cpp index 453c026e9..51a4bb049 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -89,18 +89,6 @@ SCENARIO("PortableMDArrays can be allocated from a pointer", } } -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; - RegularGrid1D g(min, max, N); - REQUIRE(g.min() == min); - REQUIRE(g.max() == max); - REQUIRE(g.nPoints() == N); - } -} - TEST_CASE("PiecewiseGrid1D", "[PiecewiseGrid1D]") { GIVEN("Some regular grid 1Ds") { RegularGrid1D g1(0, 0.25, 3); From 222bdae55006064b0099d1b552e69e2a2c01a6da Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 15:38:17 -0600 Subject: [PATCH 19/48] Add some testing with transformations. --- spiner/regular_grid_1d.hpp | 2 -- test/regular_grid_1d.cpp | 58 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 2 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index 32e0bae9d..b76a404b3 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -45,8 +45,6 @@ struct weights_t { } }; -// TODO: Need to expand testing - // 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). diff --git a/test/regular_grid_1d.cpp b/test/regular_grid_1d.cpp index 51dd78ac1..a02dd2eb5 100644 --- a/test/regular_grid_1d.cpp +++ b/test/regular_grid_1d.cpp @@ -16,6 +16,11 @@ #include #include +#include + +#include + +using Catch::Matchers::WithinRel; TEST_CASE("RegularGrid1D", "[RegularGrid1D]") { SECTION("A regular grid 1d emits appropriate metadata") { @@ -29,3 +34,56 @@ TEST_CASE("RegularGrid1D", "[RegularGrid1D]") { } } +TEST_CASE("RegularGrid1D with transformations", "[RegularGrid1D]") { + constexpr double min = 1; + constexpr double max = 1024; + constexpr size_t N = 11; + Spiner::RegularGrid1D glin(min, max, N); + Spiner::RegularGrid1D 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 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 wlog; + glog.weights(xlog, ilog, wlog); + CHECK(ilog == n - 1); + CHECK_THAT(wlog.first, WithinRel(0.5)); + CHECK_THAT(wlog.second, WithinRel(0.5)); + } + +} + From 0bc4e328df01c257a54ff6289e251acb97dcf06a Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 16:17:36 -0600 Subject: [PATCH 20/48] Add a TODO --- spiner/piecewise_grid_1d.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/spiner/piecewise_grid_1d.hpp b/spiner/piecewise_grid_1d.hpp index 49501c964..830fd7d03 100644 --- a/spiner/piecewise_grid_1d.hpp +++ b/spiner/piecewise_grid_1d.hpp @@ -32,6 +32,8 @@ #include "regular_grid_1d.hpp" +// TODO: This is based on RegularGrid1D, so thread through options for transformations + namespace Spiner { template Date: Thu, 19 Sep 2024 16:17:58 -0600 Subject: [PATCH 21/48] test tags --- test/regular_grid_1d.cpp | 2 +- test/transformations.cpp | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/test/regular_grid_1d.cpp b/test/regular_grid_1d.cpp index a02dd2eb5..6a8142866 100644 --- a/test/regular_grid_1d.cpp +++ b/test/regular_grid_1d.cpp @@ -34,7 +34,7 @@ TEST_CASE("RegularGrid1D", "[RegularGrid1D]") { } } -TEST_CASE("RegularGrid1D with transformations", "[RegularGrid1D]") { +TEST_CASE("RegularGrid1D with transformations", "[RegularGrid1D][transformations]") { constexpr double min = 1; constexpr double max = 1024; constexpr size_t N = 11; diff --git a/test/transformations.cpp b/test/transformations.cpp index 3e553cca3..b7a05133c 100644 --- a/test/transformations.cpp +++ b/test/transformations.cpp @@ -33,7 +33,7 @@ namespace { } }; -TEST_CASE("transformation: linear", "") { +TEST_CASE("transformation: linear", "[transformations]") { // This one is almost too simple to meaningfully test, but we can at least // ensure that it compiles and does the trivially-obvious things. using Transform = Spiner::TransformLinear; @@ -63,7 +63,7 @@ TEST_CASE("transformation: linear", "") { CHECK(accum / num_threads == 0); } -TEST_CASE("transformation: logarithmic", "") { +TEST_CASE("transformation: logarithmic", "[transformations]") { using Transform = Spiner::TransformLogarithmic; // Test on CPU From d45bcc68bf379fed7787bcf6f1c1819dc4dbe671 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 16:18:50 -0600 Subject: [PATCH 22/48] Add tests for DataBox with transformations. The current results seem to indicate that the independent variable transformations are correct but the depending variable transformation is incorrect. --- spiner/databox.hpp | 2 - test/test.cpp | 92 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 92 insertions(+), 2 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index b84ade6dd..a8075f8ef 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -39,8 +39,6 @@ // TODO: get named indices working // TODO: If asserts are too slow, remove them. -// TODO: Need to expand testing to cover transformations (both independent and dependent). - // TODO: Should extrapolations be this MR or a separate MR? namespace Spiner { diff --git a/test/test.cpp b/test/test.cpp index 51a4bb049..a6d85bc45 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -52,6 +52,16 @@ PORTABLE_INLINE_FUNCTION Real linearFunction(Real b, Real a, Real z, Real y, Real x) { return x + y + z + a + b; } +PORTABLE_INLINE_FUNCTION double log_lin_func(const double x, const double y) { + // log-log: f is linear function in log(x) and log(y) + // --> f = a + b * log(x) + c * log(y) + return 0.7 * std::log(x) + 0.3 * std::log(y) - 0.5; +} +PORTABLE_INLINE_FUNCTION double log_log_func(const double x, const double y) { + // log-log: log(f) is linear function in log(x) and log(y) + // --> log(f) = a + b * log(x) + c * log(y) + return std::exp(0.7 * std::log(x) + 0.3 * std::log(y) - 0.5); +} SCENARIO("PortableMDArrays can be allocated from a pointer", "[PortableMDArray]") { @@ -496,6 +506,88 @@ TEST_CASE("DataBox interpolation", "[DataBox]") { free(db); // free databox } +TEST_CASE("DataBox Interpolation with log-lin transformations", "[DataBox]") { + using Transform = Spiner::TransformLogarithmic; + using RG1D = Spiner::RegularGrid1D; + using DB = Spiner::DataBox; + + constexpr int RANK = 2; + const int Npt = 32; + const int Nsample = Npt * 16; + + const double xmin = 1; + const double xmax = 10; + + DB db(Spiner::AllocationTarget::Device, Npt, Npt); + for (int i = 0; i < RANK; i++) { + db.setRange(i, xmin, xmax, Npt); + } + + portableFor( + "fill databox", 0, Npt, 0, Npt, + PORTABLE_LAMBDA(const int ix, const int iy) { + RG1D grid(xmin, xmax, Npt); + const double x = grid.x(ix); + const double y = grid.x(iy); + db(ix, iy) = log_lin_func(x, y); + }); + + double error = 0; + portableReduce( + "interpolate databox", 0, Nsample, 0, Nsample, + PORTABLE_LAMBDA(const int ix, const int iy, Real &accumulate) { + RG1D grid(xmin, xmax, Nsample); + const double x = grid.x(ix); + const double y = grid.x(iy); + const double f_true = log_lin_func(x, y); + const double difference = db.interpToReal(x, y) - f_true; + accumulate += (difference * difference); + }, + error); + REQUIRE(error <= EPSTEST); +} + +TEST_CASE("DataBox Interpolation with log-log transformations", "[DataBox]") { + using Transform = Spiner::TransformLogarithmic; + using RG1D = Spiner::RegularGrid1D; + using DB = Spiner::DataBox; + + constexpr int RANK = 2; + const int Npt = 32; + const int Nsample = Npt * 16; + + const double xmin = 1; + const double xmax = 10; + + DB db(Spiner::AllocationTarget::Device, Npt, Npt); + for (int i = 0; i < RANK; i++) { + db.setRange(i, xmin, xmax, Npt); + } + + portableFor( + "fill databox", 0, Npt, 0, Npt, + PORTABLE_LAMBDA(const int ix, const int iy) { + RG1D grid(xmin, xmax, Npt); + const double x = grid.x(ix); + const double y = grid.x(iy); + db(ix, iy) = log_log_func(x, y); + }); + + double error = 0; + portableReduce( + "interpolate databox", 0, Nsample, 0, Nsample, + PORTABLE_LAMBDA(const int ix, const int iy, Real &accumulate) { + RG1D grid(xmin, xmax, Nsample); + const double x = grid.x(ix); + const double y = grid.x(iy); + const double f_true = log_log_func(x, y); + const double difference = db.interpToReal(x, y) - f_true; + accumulate += (difference * difference); + }, + error); + REQUIRE(error <= EPSTEST); +} + TEST_CASE("DataBox Interpolation with piecewise grids", "[DataBox][PiecewiseGrid1D]") { GIVEN("A piecewise grid") { From 37f38406a359db667adb1fac9195dbf32f829623 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Thu, 19 Sep 2024 17:08:47 -0600 Subject: [PATCH 23/48] Starting to fix operator(). Mostly lots of stream-of-consciousness notes, writing down ideas to resolve this as I thought of them. --- spiner/databox.hpp | 91 +++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 89 insertions(+), 2 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index a8075f8ef..36727ff06 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -153,14 +153,101 @@ class DataBox { resize(AllocationTarget::Host, std::forward(args)...); } + // TODO: The index operators are a problem. If you assign to them, you _think_ you're assigning + // the y value, but in actuality you're assigning the v value, which means the transformations + // will be all messed up. I can imagine two possible approaches. + // -- The "Safer" Solution ______________________________________________________________________ + // In order to avoid the problems that come with more-complex machinery and all the weirdness + // that can arise, the "safer" solution is to simply get rid of operator() and replace it with + // get() and set(). The downside to this is that it will break existing interfaces. + // -- We don't _have_ to call these get() and set(). For example, if we want to do something + // analogous to x() from RegularGrid1D, we could call them both y() or dependent_variable() or + // something like that. Then we would use dispatch to disambiguate: y(i,j) is a getter but + // y(value, i, j) is a setter. Of course, if T is the same type as the indices, then it + // becomes ambiguous whether you mean y(i,j,k) or y(value,i,j). And with variadic parameter + // packs, it may always be ambiguous even if T is not the same type as the indices. + // -- The "Sneaky" Solution _____________________________________________________________________ + // If we want to preserve the existing interface (where you can assign to db(i,j) and also + // read from db(i,j)), we can create a hidden helper class: db_value_ref. It is templated on + // T and Transform the same as DataBox, and it holds a reference (or pointer if you prefer) to + // a single value within the DataBox's data array. It implements two functions; + // -- operator=(const T new_value) will handle statements of the form `db(i,j) = new_value` + // -- implicit cast to T will handle statements that read the value of db(i,j) + // Problems I can already see with this approach: + // -- Our data types have more than just `operator=`; for example, `double` has `+=`, `*=`, + // `++`, etc, so to truly hide this we would need to implement all such operators. But + // since we don't _know_ the data type, we would have to use some careful SFINAE to + // implement them if they work for the type T but not compile them if they don't work for + // the type T. + // -- Because of the nature of generic programming, there could well be a lot of `auto value = + // db(i,j)`, and then using that to do some things. We can't easily control the context + // where `value` is then used, so we can't guarantee that the compiler will recognize that + // something should be an implicit cast to type T. That means some statements could break. + // It also means some statements could work, but you end up with things like + // `std::vector` when you want `std::vector`, which could cause + // significant problems somewhere else down the line. + // Because of the dangers in the "sneaky" solution, I'm more inclined towards the "safer" + // solution because we can change that more quickly and with little chance of failure (whereas + // the "sneaky" solution would take longer and may not even be workable without some serious C++ + // wizardry, or possibly even possible at all). But it does mean changing the interface for our + // users, which is unfortunate. + // We could, in theory, keep operator() as they are. Instead of being an accessor to the y + // values, they're now accessors to the v values (recall: v = Transform::forward(y)). This is + // not inherently wrong so long as we're okay with the idea of exposing the u and v (transformed) + // values, which we may not be. But where this idea runs into serious problems is that we would + // be keeping an interface but changing its meaning, which means that user code changes but + // doesn't necessarily break. And that's a big problem. + // Another option to preserve functionality, but which is a hack that has problems of its + // own, would be to keep operator() as they are but use SFINAE to only enable them if the + // transform template is something we know to be a no-op transformation (that is, we would have + // to hard-code it explicitly to TransformLinear, and a user-defined equivalent wouldn't work). + // The advantage to this is that existing code will still work: existing code won't specify a + // transformation, so it'll default to TransformLinear, so operator() will be enabled and will + // effectively mean what it always meant. But if the user starts using the new feature and adds + // any transformation other than the default, then operator() will simply disappear, which will + // prevent the user from accidentally mis-using operator() due to not understanding how it + // changed. It hides the change and keeps existing code using and is safe. But it does mean + // that the interface to DataBox changes depending on the template parameters. If a user writes + // code that has no transformation and uses operator(), then decides to add a transformation + // later on, they may in for a rude surprise when they have to go through their code and + // carefully change every operator() to get() and/or set(). But we could then highlight this + // problem and deprecate operator(), giving users time to shift their code, and then eventually + // delete operator() completely in some future release, forcing users to adopt the (hopefully + // by-then long-existing and commonly-used) get() / set() paradigm instead of the operator() + // paradigm. + // -- This may run into issues with a variadic template. But since the argument list will define + // the variadic template, then hopefully the "leftover" SFINAE will be allowed. Or maybe the + // SFINAE goes first. Or maybe the SFINAE goes on the return type. I'll have to experiment. + // -- We could get around this by explicitly typing out the different numbers of indices, because + // there's a finite set (that is, MAXRANK is known; I think it's currently 6). + // If you want to play _even more_ SFINAE shenanigans, you could keep `operator()(indices...) + // const` to always be a getter and then write `operator()(value, indices...)` to be a setter. + // Then you SFINAE on the Transform to turn on or off `operator()(indices...)` as a _setter_. + // But this gets complicated and confusing, and I think different function names is a smarter + // choice. Plus, as noted above, if the getter and setter have the same name, then there is the + // possibility of ambiguity between `operator()(i,j,k)` (a getter) and `operator()(value,i,j)` (a + // setter). + + // Data array accessors + template + PORTABLE_INLINE_FUNCTION T get(Args... args) { + return Transform::reverse(dataView_(std::forward(args)...)); + } + template + PORTABLE_INLINE_FUNCTION void set(const T new_value, Args... args) { + dataView_(std::forward(args)...) = Transform::forward(new_value); + } + // Index operators // examle calls: // T x = db(n4, n3, n2, n1); - template + template ::value>::type* = nullptr> PORTABLE_INLINE_FUNCTION T &operator()(Args... args) { return dataView_(std::forward(args)...); } - template + template ::value>::type* = nullptr> PORTABLE_INLINE_FUNCTION T &operator()(Args... args) const { return dataView_(std::forward(args)...); } From 958062e388fb181c8beb860eafcede6b1a1d1753 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 10:23:51 -0600 Subject: [PATCH 24/48] Fix some things about accessing the dependent variable values. --- spiner/databox.hpp | 19 +++++++++++++------ test/test.cpp | 4 ++-- 2 files changed, 15 insertions(+), 8 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index 36727ff06..21ae0fd1f 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -230,24 +230,31 @@ class DataBox { // Data array accessors template - PORTABLE_INLINE_FUNCTION T get(Args... args) { + PORTABLE_INLINE_FUNCTION T get_data_value(Args... args) const { return Transform::reverse(dataView_(std::forward(args)...)); } + // TODO: Should this method be const? Having this be const is in line with the second operator() + // below (the one with const), which allows the user to modify the dependent variable + // values of a const DataBox. template - PORTABLE_INLINE_FUNCTION void set(const T new_value, Args... args) { + PORTABLE_INLINE_FUNCTION void set_data_value(const T new_value, Args... args) const { dataView_(std::forward(args)...) = Transform::forward(new_value); } // Index operators // examle calls: // T x = db(n4, n3, n2, n1); - template ::value>::type* = nullptr> + template < + typename... Args, + typename transform_t=Transform, // only for SFINAE + typename std::enable_if::value>::type* = nullptr> PORTABLE_INLINE_FUNCTION T &operator()(Args... args) { return dataView_(std::forward(args)...); } - template ::value>::type* = nullptr> + template < + typename... Args, + typename transform_t=Transform, // only for SFINAE + typename std::enable_if::value>::type* = nullptr> PORTABLE_INLINE_FUNCTION T &operator()(Args... args) const { return dataView_(std::forward(args)...); } diff --git a/test/test.cpp b/test/test.cpp index a6d85bc45..9fe5c4e34 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -550,7 +550,7 @@ TEST_CASE("DataBox Interpolation with log-lin transformations", "[DataBox]") { TEST_CASE("DataBox Interpolation with log-log transformations", "[DataBox]") { using Transform = Spiner::TransformLogarithmic; using RG1D = Spiner::RegularGrid1D; - using DB = Spiner::DataBox; + using DB = Spiner::DataBox; constexpr int RANK = 2; const int Npt = 32; @@ -570,7 +570,7 @@ TEST_CASE("DataBox Interpolation with log-log transformations", "[DataBox]") { RG1D grid(xmin, xmax, Npt); const double x = grid.x(ix); const double y = grid.x(iy); - db(ix, iy) = log_log_func(x, y); + db.set_data_value(log_log_func(x,y), ix, iy); }); double error = 0; From 57ddf50122379bf56dec35eb5b1ec79b5f20b9a2 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 10:37:58 -0600 Subject: [PATCH 25/48] Revise discussion of operator() and accessors. --- spiner/databox.hpp | 112 +++++++++++++++------------------------------ 1 file changed, 38 insertions(+), 74 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index 21ae0fd1f..b199459fe 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -153,80 +153,44 @@ class DataBox { resize(AllocationTarget::Host, std::forward(args)...); } - // TODO: The index operators are a problem. If you assign to them, you _think_ you're assigning - // the y value, but in actuality you're assigning the v value, which means the transformations - // will be all messed up. I can imagine two possible approaches. - // -- The "Safer" Solution ______________________________________________________________________ - // In order to avoid the problems that come with more-complex machinery and all the weirdness - // that can arise, the "safer" solution is to simply get rid of operator() and replace it with - // get() and set(). The downside to this is that it will break existing interfaces. - // -- We don't _have_ to call these get() and set(). For example, if we want to do something - // analogous to x() from RegularGrid1D, we could call them both y() or dependent_variable() or - // something like that. Then we would use dispatch to disambiguate: y(i,j) is a getter but - // y(value, i, j) is a setter. Of course, if T is the same type as the indices, then it - // becomes ambiguous whether you mean y(i,j,k) or y(value,i,j). And with variadic parameter - // packs, it may always be ambiguous even if T is not the same type as the indices. - // -- The "Sneaky" Solution _____________________________________________________________________ - // If we want to preserve the existing interface (where you can assign to db(i,j) and also - // read from db(i,j)), we can create a hidden helper class: db_value_ref. It is templated on - // T and Transform the same as DataBox, and it holds a reference (or pointer if you prefer) to - // a single value within the DataBox's data array. It implements two functions; - // -- operator=(const T new_value) will handle statements of the form `db(i,j) = new_value` - // -- implicit cast to T will handle statements that read the value of db(i,j) - // Problems I can already see with this approach: - // -- Our data types have more than just `operator=`; for example, `double` has `+=`, `*=`, - // `++`, etc, so to truly hide this we would need to implement all such operators. But - // since we don't _know_ the data type, we would have to use some careful SFINAE to - // implement them if they work for the type T but not compile them if they don't work for - // the type T. - // -- Because of the nature of generic programming, there could well be a lot of `auto value = - // db(i,j)`, and then using that to do some things. We can't easily control the context - // where `value` is then used, so we can't guarantee that the compiler will recognize that - // something should be an implicit cast to type T. That means some statements could break. - // It also means some statements could work, but you end up with things like - // `std::vector` when you want `std::vector`, which could cause - // significant problems somewhere else down the line. - // Because of the dangers in the "sneaky" solution, I'm more inclined towards the "safer" - // solution because we can change that more quickly and with little chance of failure (whereas - // the "sneaky" solution would take longer and may not even be workable without some serious C++ - // wizardry, or possibly even possible at all). But it does mean changing the interface for our - // users, which is unfortunate. - // We could, in theory, keep operator() as they are. Instead of being an accessor to the y - // values, they're now accessors to the v values (recall: v = Transform::forward(y)). This is - // not inherently wrong so long as we're okay with the idea of exposing the u and v (transformed) - // values, which we may not be. But where this idea runs into serious problems is that we would - // be keeping an interface but changing its meaning, which means that user code changes but - // doesn't necessarily break. And that's a big problem. - // Another option to preserve functionality, but which is a hack that has problems of its - // own, would be to keep operator() as they are but use SFINAE to only enable them if the - // transform template is something we know to be a no-op transformation (that is, we would have - // to hard-code it explicitly to TransformLinear, and a user-defined equivalent wouldn't work). - // The advantage to this is that existing code will still work: existing code won't specify a - // transformation, so it'll default to TransformLinear, so operator() will be enabled and will - // effectively mean what it always meant. But if the user starts using the new feature and adds - // any transformation other than the default, then operator() will simply disappear, which will - // prevent the user from accidentally mis-using operator() due to not understanding how it - // changed. It hides the change and keeps existing code using and is safe. But it does mean - // that the interface to DataBox changes depending on the template parameters. If a user writes - // code that has no transformation and uses operator(), then decides to add a transformation - // later on, they may in for a rude surprise when they have to go through their code and - // carefully change every operator() to get() and/or set(). But we could then highlight this - // problem and deprecate operator(), giving users time to shift their code, and then eventually - // delete operator() completely in some future release, forcing users to adopt the (hopefully - // by-then long-existing and commonly-used) get() / set() paradigm instead of the operator() - // paradigm. - // -- This may run into issues with a variadic template. But since the argument list will define - // the variadic template, then hopefully the "leftover" SFINAE will be allowed. Or maybe the - // SFINAE goes first. Or maybe the SFINAE goes on the return type. I'll have to experiment. - // -- We could get around this by explicitly typing out the different numbers of indices, because - // there's a finite set (that is, MAXRANK is known; I think it's currently 6). - // If you want to play _even more_ SFINAE shenanigans, you could keep `operator()(indices...) - // const` to always be a getter and then write `operator()(value, indices...)` to be a setter. - // Then you SFINAE on the Transform to turn on or off `operator()(indices...)` as a _setter_. - // But this gets complicated and confusing, and I think different function names is a smarter - // choice. Plus, as noted above, if the getter and setter have the same name, then there is the - // possibility of ambiguity between `operator()(i,j,k)` (a getter) and `operator()(value,i,j)` (a - // setter). + // TODO: The existence and implementation of operator() is a problem. + // -- If you use operator() to assign to the dependent variable, you _think_ you're assigning + // to the y (untransformed) value, because that was the previous behavior. But in reality, + // you're assigning to the v (transformed) value, which means that you're not storing the + // right value unless you know this quirk. But that means that we're no longer hiding the + // transformations, but exposing things that users won't necessarily know about. It's also + // surprising (I wrote the code and made this mistake already myself). + // -- We could try to be clever by writing a value_reference_t class that's returned by + // operator() and has customized operators to hide the transformations. This would be a + // non-trivial amount of work, and I can imagine a number of ways that this would break. + // -- Instead, I've opted for a "safer" option: write new methods get_data_value and + // set_data_value (I'm not attached to these names), which are classic accessors to the + // dependent variable values, which allow the DataBox to handle the transformations + // correctly. + // -- The next question is how to handle existing code, which already uses operator(). + // -- If we add the accessors but leave operator(), we still have the surprising behavior. + // Users will naturally tend to use the feature they're used to (operator()) and will + // run into trouble without understanding why. + // -- If we remove operator() to avoid causing problems, we break existing code. As much + // as possible, I'd like to hide changes so that existing code continues to work. It's + // not always possible, but it's preferable when possible. + // -- I've opted for a compromise solution: + // -- Leave operator() in place, but add some SFINAE so that it only works if you use + // the default transformation (TransformLinear) where y = v (that is, the transform + // is a no-op). This will preserve existing code. + // -- Existing code also has the option to use the new accessors. + // -- If you add a transformation, operator() is removed and you _must_ use the new + // accessors. In this situation, you're writing new code using new features, so the + // "keep existing code working as-is" argument doesn't apply. + // -- The trickiest situation is a user adapting existing code by adding a + // transformation. We'll have to clearly communicate the behavior here so that they + // understand what's going on (documentation, release notes, etc.). + // -- In the future, operator() _could_ be removed from a later release of Spiner, + // although this is not required. Personal opinion: it should (eventually) be + // removed. In the case of an interpolator, operator() would most naturally be + // assumed to be the method that does the interpolation (I ran into this problem when + // I first used Spiner), so I think operator() being a way to access the data points + // isn't the most intuitive way to name things. // Data array accessors template From 4734797303974373dc7e79a88a3854064e21970b Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 10:44:44 -0600 Subject: [PATCH 26/48] Remove long TODO discussion, because it's better as a comment on the PR. --- spiner/databox.hpp | 39 --------------------------------------- 1 file changed, 39 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index b199459fe..b846bc292 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -153,45 +153,6 @@ class DataBox { resize(AllocationTarget::Host, std::forward(args)...); } - // TODO: The existence and implementation of operator() is a problem. - // -- If you use operator() to assign to the dependent variable, you _think_ you're assigning - // to the y (untransformed) value, because that was the previous behavior. But in reality, - // you're assigning to the v (transformed) value, which means that you're not storing the - // right value unless you know this quirk. But that means that we're no longer hiding the - // transformations, but exposing things that users won't necessarily know about. It's also - // surprising (I wrote the code and made this mistake already myself). - // -- We could try to be clever by writing a value_reference_t class that's returned by - // operator() and has customized operators to hide the transformations. This would be a - // non-trivial amount of work, and I can imagine a number of ways that this would break. - // -- Instead, I've opted for a "safer" option: write new methods get_data_value and - // set_data_value (I'm not attached to these names), which are classic accessors to the - // dependent variable values, which allow the DataBox to handle the transformations - // correctly. - // -- The next question is how to handle existing code, which already uses operator(). - // -- If we add the accessors but leave operator(), we still have the surprising behavior. - // Users will naturally tend to use the feature they're used to (operator()) and will - // run into trouble without understanding why. - // -- If we remove operator() to avoid causing problems, we break existing code. As much - // as possible, I'd like to hide changes so that existing code continues to work. It's - // not always possible, but it's preferable when possible. - // -- I've opted for a compromise solution: - // -- Leave operator() in place, but add some SFINAE so that it only works if you use - // the default transformation (TransformLinear) where y = v (that is, the transform - // is a no-op). This will preserve existing code. - // -- Existing code also has the option to use the new accessors. - // -- If you add a transformation, operator() is removed and you _must_ use the new - // accessors. In this situation, you're writing new code using new features, so the - // "keep existing code working as-is" argument doesn't apply. - // -- The trickiest situation is a user adapting existing code by adding a - // transformation. We'll have to clearly communicate the behavior here so that they - // understand what's going on (documentation, release notes, etc.). - // -- In the future, operator() _could_ be removed from a later release of Spiner, - // although this is not required. Personal opinion: it should (eventually) be - // removed. In the case of an interpolator, operator() would most naturally be - // assumed to be the method that does the interpolation (I ran into this problem when - // I first used Spiner), so I think operator() being a way to access the data points - // isn't the most intuitive way to name things. - // Data array accessors template PORTABLE_INLINE_FUNCTION T get_data_value(Args... args) const { From ed811b558be7638f2cc54d5e1c463a8ae4e56385 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 10:59:52 -0600 Subject: [PATCH 27/48] Add transformations to PiecewiseGrid1D. --- spiner/piecewise_grid_1d.hpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/spiner/piecewise_grid_1d.hpp b/spiner/piecewise_grid_1d.hpp index 830fd7d03..206a87ab1 100644 --- a/spiner/piecewise_grid_1d.hpp +++ b/spiner/piecewise_grid_1d.hpp @@ -32,11 +32,10 @@ #include "regular_grid_1d.hpp" -// TODO: This is based on RegularGrid1D, so thread through options for transformations - namespace Spiner { template ::value, bool>::type> class PiecewiseGrid1D { @@ -48,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> grids) { + PiecewiseGrid1D(const std::vector> grids) { NGRIDS_ = grids.size(); PORTABLE_ALWAYS_REQUIRE( NGRIDS_ <= NGRIDSMAX, @@ -67,8 +66,8 @@ class PiecewiseGrid1D { } } } - PiecewiseGrid1D(std::initializer_list> grids) - : PiecewiseGrid1D(std::vector>(grids)) {} + PiecewiseGrid1D(std::initializer_list> grids) + : PiecewiseGrid1D(std::vector>(grids)) {} template PORTABLE_INLINE_FUNCTION int findGrid(const F &direction) const { @@ -126,14 +125,14 @@ class PiecewiseGrid1D { } PORTABLE_INLINE_FUNCTION - bool operator==(const PiecewiseGrid1D &other) const { + bool operator==(const PiecewiseGrid1D &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 &other) const { + bool operator!=(const PiecewiseGrid1D &other) const { return !(*this == other); } PORTABLE_INLINE_FUNCTION T min() const { return grids_[0].min(); } @@ -227,7 +226,7 @@ class PiecewiseGrid1D { SP5::H1D::GRID_FORMAT[1]; } - RegularGrid1D grids_[NGRIDSMAX]; + RegularGrid1D grids_[NGRIDSMAX]; int pointTotals_[NGRIDSMAX]; int NGRIDS_; }; From 78594ea746ada6fec2b15c89ec2d8cf7300dc7c9 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 11:05:18 -0600 Subject: [PATCH 28/48] We'll push extrapolations to another MR. --- spiner/databox.hpp | 2 -- spiner/databox_wrapper.hpp | 5 ----- 2 files changed, 7 deletions(-) delete mode 100644 spiner/databox_wrapper.hpp diff --git a/spiner/databox.hpp b/spiner/databox.hpp index b846bc292..de0872acc 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -39,8 +39,6 @@ // TODO: get named indices working // TODO: If asserts are too slow, remove them. -// TODO: Should extrapolations be this MR or a separate MR? - namespace Spiner { enum class IndexType { Interpolated = 0, Named = 1, Indexed = 2 }; diff --git a/spiner/databox_wrapper.hpp b/spiner/databox_wrapper.hpp deleted file mode 100644 index af0eb1319..000000000 --- a/spiner/databox_wrapper.hpp +++ /dev/null @@ -1,5 +0,0 @@ -// The original use-case did three things: -// * Add transformations of the data. This was originally hard-coded, but we'll make it general. -// * Add extrapolations off the edge of the table. -// * Adapt the interface. We should instead follow the Databox interface. - From cac4be1c1dbc05da241edc2d3c14323560cf2c10 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 11:27:48 -0600 Subject: [PATCH 29/48] Moving some TODO comments to MR discussion. --- spiner/regular_grid_1d.hpp | 29 ----------------------------- 1 file changed, 29 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index b76a404b3..cb5531f3f 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -45,21 +45,6 @@ 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 = @@ -132,22 +117,10 @@ class RegularGrid1D { operator!=(const RegularGrid1D &other) const { return !(*this == other); } - // 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(umin_) || std::isnan(umax_) || std::isnan(du_) || @@ -155,8 +128,6 @@ class RegularGrid1D { } 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( From 8ab5205bb48db0774f1b01ccfc903657d19edf0d Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 11:29:49 -0600 Subject: [PATCH 30/48] Moving some TODO comments to MR discussion. --- spiner/databox.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index de0872acc..5e352c195 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -156,9 +156,6 @@ class DataBox { PORTABLE_INLINE_FUNCTION T get_data_value(Args... args) const { return Transform::reverse(dataView_(std::forward(args)...)); } - // TODO: Should this method be const? Having this be const is in line with the second operator() - // below (the one with const), which allows the user to modify the dependent variable - // values of a const DataBox. template PORTABLE_INLINE_FUNCTION void set_data_value(const T new_value, Args... args) const { dataView_(std::forward(args)...) = Transform::forward(new_value); From a178359c0866931e7a391cddb1d40f1d79696d57 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 11:34:37 -0600 Subject: [PATCH 31/48] Edit TODO. --- spiner/transformations.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index dfac62466..aca278eb9 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -51,9 +51,8 @@ namespace Spiner { } }; - // TODO: log_NQT - - // TODO: arcsinh_NQT + // 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_ From c7afa649af3dfa300292c9bafda8adf5bd1ab55d Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 11:59:42 -0600 Subject: [PATCH 32/48] format --- spiner/databox.hpp | 114 +++++++++++++++++-------------- spiner/regular_grid_1d.hpp | 14 ++-- spiner/transformations.hpp | 65 +++++++++--------- test/regular_grid_1d.cpp | 95 +++++++++++++------------- test/test.cpp | 136 ++++++++++++++++++------------------- test/transformations.cpp | 135 ++++++++++++++++++------------------ 6 files changed, 284 insertions(+), 275 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index 5e352c195..eae25e410 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -114,8 +114,8 @@ class DataBox { // Slice constructor PORTABLE_INLINE_FUNCTION - DataBox(const DataBox &b, const int dim, const int indx, - const int nvar) noexcept + DataBox(const DataBox &b, const int dim, + const int indx, const int nvar) noexcept : status_(DataStatus::Unmanaged), data_(b.data_) { dataView_.InitWithShallowSlice(b.dataView_, dim, indx, nvar); rank_ = dataView_.GetRank(); @@ -154,27 +154,28 @@ class DataBox { // Data array accessors template PORTABLE_INLINE_FUNCTION T get_data_value(Args... args) const { - return Transform::reverse(dataView_(std::forward(args)...)); + return Transform::reverse(dataView_(std::forward(args)...)); } - template - PORTABLE_INLINE_FUNCTION void set_data_value(const T new_value, Args... args) const { + template + PORTABLE_INLINE_FUNCTION void set_data_value(const T new_value, + Args... args) const { dataView_(std::forward(args)...) = Transform::forward(new_value); } // Index operators // examle calls: // T x = db(n4, n3, n2, n1); - template < - typename... Args, - typename transform_t=Transform, // only for SFINAE - typename std::enable_if::value>::type* = nullptr> + template ::value>::type * = nullptr> PORTABLE_INLINE_FUNCTION T &operator()(Args... args) { return dataView_(std::forward(args)...); } - template < - typename... Args, - typename transform_t=Transform, // only for SFINAE - typename std::enable_if::value>::type* = nullptr> + template ::value>::type * = nullptr> PORTABLE_INLINE_FUNCTION T &operator()(Args... args) const { return dataView_(std::forward(args)...); } @@ -182,7 +183,7 @@ class DataBox { // Slice operation PORTABLE_INLINE_FUNCTION DataBox slice(const int dim, const int indx, - const int nvar) const { + const int nvar) const { return DataBox(*this, dim, indx, nvar); } PORTABLE_INLINE_FUNCTION DataBox @@ -232,7 +233,8 @@ class DataBox { PORTABLE_INLINE_FUNCTION void interpFromDB(const DataBox &db, const T x); PORTABLE_INLINE_FUNCTION void - interpFromDB(const DataBox &db, const T x2, const T x1); + interpFromDB(const DataBox &db, const T x2, + const T x1); template PORTABLE_INLINE_FUNCTION DataBox interpToDB(Args... args) { @@ -262,8 +264,9 @@ class DataBox { // Does no checks that memory is available. // Optionally copies shape of source with ndims fewer slowest-moving // dimensions - PORTABLE_INLINE_FUNCTION void copyShape(const DataBox &db, - const int ndims = 0); + PORTABLE_INLINE_FUNCTION void + copyShape(const DataBox &db, + const int ndims = 0); // Returns new databox with same memory and metadata inline void copyMetadata(const DataBox &src); @@ -291,8 +294,10 @@ class DataBox { inline bool ownsAllocatedMemory() { return (status_ != DataStatus::Unmanaged); } - inline bool operator==(const DataBox &other) const; - inline bool operator!=(const DataBox &other) const { + inline bool + operator==(const DataBox &other) const; + inline bool + operator!=(const DataBox &other) const { return !(*this == other); } @@ -399,8 +404,8 @@ class DataBox { // copy to device portableCopyToDevice(device_data, data_, sizeBytes()); // create new databox of size size - DataBox a{device_data, dim(6), dim(5), dim(4), - dim(3), dim(2), dim(1)}; + DataBox a{ + device_data, dim(6), dim(5), dim(4), dim(3), dim(2), dim(1)}; a.copyShape(*this); // set correct allocation status of the new databox // note this is ALWAYS device, even if host==device. @@ -457,7 +462,8 @@ class DataBox { // Read an array, shallow template -inline void DataBox::setArray(PortableMDArray &A) { +inline void +DataBox::setArray(PortableMDArray &A) { dataView_ = A; rank_ = A.GetRank(); setAllIndexed_(); @@ -475,7 +481,8 @@ DataBox::interpToReal(const T x) const noexcept { } template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +PORTABLE_FORCEINLINE_FUNCTION T +DataBox::interpToReal( const T x2, const T x1) const noexcept { assert(canInterpToReal_(2)); int ix1, ix2; @@ -484,15 +491,16 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo grids_[1].weights(x2, ix2, w2); // TODO: prefectch corners for speed? // TODO: re-order access pattern? - const auto v = (w2[0] * - (w1[0] * dataView_(ix2, ix1) + w1[1] * dataView_(ix2, ix1 + 1)) + - w2[1] * (w1[0] * dataView_(ix2 + 1, ix1) + - w1[1] * dataView_(ix2 + 1, ix1 + 1))); + const auto v = + (w2[0] * (w1[0] * dataView_(ix2, ix1) + w1[1] * dataView_(ix2, ix1 + 1)) + + w2[1] * (w1[0] * dataView_(ix2 + 1, ix1) + + w1[1] * dataView_(ix2 + 1, ix1 + 1))); return Transform::reverse(v); } template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +PORTABLE_FORCEINLINE_FUNCTION T +DataBox::interpToReal( const T x3, const T x2, const T x1) const noexcept { assert(canInterpToReal_(3)); int ix[3]; @@ -516,7 +524,8 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo } template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +PORTABLE_FORCEINLINE_FUNCTION T +DataBox::interpToReal( const T x3, const T x2, const T x1, const int idx) const noexcept { assert(rank_ == 4); for (int r = 1; r < rank_; ++r) { @@ -548,7 +557,8 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo // DH: this is a large function to force an inline, perhaps just make it a // suggestion to the compiler? template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +PORTABLE_FORCEINLINE_FUNCTION T +DataBox::interpToReal( const T x4, const T x3, const T x2, const T x1) const noexcept { assert(canInterpToReal_(4)); T x[] = {x1, x2, x3, x4}; @@ -600,7 +610,8 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo } template -PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpToReal( +PORTABLE_FORCEINLINE_FUNCTION T +DataBox::interpToReal( const T x4, const T x3, const T x2, const int idx, const T x1) const noexcept { assert(rank_ == 5); @@ -665,8 +676,8 @@ PORTABLE_FORCEINLINE_FUNCTION T DataBox::interpTo template PORTABLE_INLINE_FUNCTION void -DataBox::interpFromDB(const DataBox &db, - const T x) { +DataBox::interpFromDB( + const DataBox &db, const T x) { assert(db.indices_[db.rank_ - 1] == IndexType::Interpolated); assert(db.grids_[db.rank_ - 1].isWellFormed()); assert(size() == (db.size() / db.dim(db.rank_))); @@ -676,7 +687,8 @@ DataBox::interpFromDB(const DataBox lower(db.slice(ix)), upper(db.slice(ix + 1)); + DataBox lower(db.slice(ix)), + upper(db.slice(ix + 1)); // lower = db.slice(ix); // upper = db.slice(ix+1); for (int i = 0; i < size(); i++) { @@ -686,8 +698,8 @@ DataBox::interpFromDB(const DataBox PORTABLE_INLINE_FUNCTION void -DataBox::interpFromDB(const DataBox &db, - const T x2, const T x1) { +DataBox::interpFromDB( + const DataBox &db, const T x2, const T x1) { assert(db.rank_ >= 2); assert(db.indices_[db.rank_ - 1] == IndexType::Interpolated); assert(db.grids_[db.rank_ - 1].isWellFormed()); @@ -731,9 +743,8 @@ DataBox::interpFromDB(const DataBox -PORTABLE_INLINE_FUNCTION void -DataBox::copyShape(const DataBox &db, - const int ndims) { +PORTABLE_INLINE_FUNCTION void DataBox::copyShape( + const DataBox &db, const int ndims) { rank_ = db.rank_ - ndims; int dims[MAXRANK]; for (int i = 0; i < MAXRANK; i++) @@ -767,8 +778,8 @@ inline void DataBox::copyMetadata( #ifdef SPINER_USE_HDF template -inline herr_t -DataBox::saveHDF(const std::string &filename) const { +inline herr_t DataBox::saveHDF( + const std::string &filename) const { herr_t status; hid_t file; @@ -779,9 +790,8 @@ DataBox::saveHDF(const std::string &filename) con } template -inline herr_t -DataBox::saveHDF(hid_t loc, - const std::string &groupname) const { +inline herr_t DataBox::saveHDF( + hid_t loc, const std::string &groupname) const { hid_t group, grids; herr_t status = 0; static_assert(std::is_same::value || std::is_same::value, @@ -851,7 +861,8 @@ DataBox::loadHDF(const std::string &filename) { template inline herr_t -DataBox::loadHDF(hid_t loc, const std::string &groupname) { +DataBox::loadHDF(hid_t loc, + const std::string &groupname) { hid_t group, grids; herr_t status = 0; std::vector index_types; @@ -908,7 +919,8 @@ DataBox::loadHDF(hid_t loc, const std::string &gr // Performs shallow copy by default template PORTABLE_INLINE_FUNCTION DataBox & -DataBox::operator=(const DataBox &src) { +DataBox::operator=( + const DataBox &src) { if (this != &src) { rank_ = src.rank_; status_ = src.status_; @@ -924,8 +936,8 @@ DataBox::operator=(const DataBox -inline void -DataBox::copy(const DataBox &src) { +inline void DataBox::copy( + const DataBox &src) { copyMetadata(src); for (int i = 0; i < src.size(); i++) dataView_(i) = src(i); @@ -973,7 +985,8 @@ DataBox::range(int i) const { } template -PORTABLE_INLINE_FUNCTION void DataBox::setAllIndexed_() { +PORTABLE_INLINE_FUNCTION void +DataBox::setAllIndexed_() { for (int i = 0; i < rank_; i++) { indices_[i] = IndexType::Indexed; } @@ -981,7 +994,8 @@ PORTABLE_INLINE_FUNCTION void DataBox::setAllInde template PORTABLE_INLINE_FUNCTION bool -DataBox::canInterpToReal_(const int interpOrder) const { +DataBox::canInterpToReal_( + const int interpOrder) const { if (rank_ != interpOrder) return false; for (int i = 0; i < rank_; i++) { if (indices_[i] != IndexType::Interpolated) return false; diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index cb5531f3f..032a59601 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -45,8 +45,7 @@ struct weights_t { } }; -template ::value, bool>::type = true> class RegularGrid1D { @@ -59,11 +58,8 @@ class RegularGrid1D { PORTABLE_INLINE_FUNCTION RegularGrid1D() : 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) { + : 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"); } @@ -83,7 +79,9 @@ class RegularGrid1D { } // Translate between x coordinate and index - PORTABLE_INLINE_FUNCTION T x(const int i) const { return Transform::reverse(u(i)); } + PORTABLE_INLINE_FUNCTION T x(const int i) const { + return Transform::reverse(u(i)); + } PORTABLE_INLINE_FUNCTION int index(const T x) const { return index_u(Transform::forward(x)); } diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index aca278eb9..03139d41b 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -21,38 +21,39 @@ 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 - PORTABLE_INLINE_FUNCTION static T forward(const T x) { - return x; - } - template - PORTABLE_INLINE_FUNCTION static T reverse(const T u) { - return u; - } - }; - - // logarithmic transformation: y = log(x + small) - struct TransformLogarithmic{ - template - PORTABLE_INLINE_FUNCTION static T forward(const T x) { - return std::log(x + std::numeric_limits::denorm_min()); - } - template - PORTABLE_INLINE_FUNCTION static T reverse(const T u) { - return std::exp(u) - std::numeric_limits::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 ;-) +// 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 + PORTABLE_INLINE_FUNCTION static T forward(const T x) { + return x; + } + template + PORTABLE_INLINE_FUNCTION static T reverse(const T u) { + return u; + } +}; + +// logarithmic transformation: y = log(x + small) +struct TransformLogarithmic { + template + PORTABLE_INLINE_FUNCTION static T forward(const T x) { + return std::log(x + std::numeric_limits::denorm_min()); + } + template + PORTABLE_INLINE_FUNCTION static T reverse(const T u) { + return std::exp(u) - std::numeric_limits::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_ diff --git a/test/regular_grid_1d.cpp b/test/regular_grid_1d.cpp index 6a8142866..5beaf5853 100644 --- a/test/regular_grid_1d.cpp +++ b/test/regular_grid_1d.cpp @@ -34,56 +34,55 @@ TEST_CASE("RegularGrid1D", "[RegularGrid1D]") { } } -TEST_CASE("RegularGrid1D with transformations", "[RegularGrid1D][transformations]") { - constexpr double min = 1; - constexpr double max = 1024; - constexpr size_t N = 11; - Spiner::RegularGrid1D glin(min, max, N); - Spiner::RegularGrid1D 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); +TEST_CASE("RegularGrid1D with transformations", + "[RegularGrid1D][transformations]") { + constexpr double min = 1; + constexpr double max = 1024; + constexpr size_t N = 11; + Spiner::RegularGrid1D glin(min, max, N); + Spiner::RegularGrid1D 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 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 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 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 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 (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 wlog; - glog.weights(xlog, ilog, wlog); - CHECK(ilog == n - 1); - CHECK_THAT(wlog.first, WithinRel(0.5)); - CHECK_THAT(wlog.second, WithinRel(0.5)); - } + // 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 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 wlog; + glog.weights(xlog, ilog, wlog); + CHECK(ilog == n - 1); + CHECK_THAT(wlog.first, WithinRel(0.5)); + CHECK_THAT(wlog.second, WithinRel(0.5)); + } } - diff --git a/test/test.cpp b/test/test.cpp index 9fe5c4e34..a4ddc6eb0 100644 --- a/test/test.cpp +++ b/test/test.cpp @@ -507,85 +507,85 @@ TEST_CASE("DataBox interpolation", "[DataBox]") { } TEST_CASE("DataBox Interpolation with log-lin transformations", "[DataBox]") { - using Transform = Spiner::TransformLogarithmic; - using RG1D = Spiner::RegularGrid1D; - using DB = Spiner::DataBox; + using Transform = Spiner::TransformLogarithmic; + using RG1D = Spiner::RegularGrid1D; + using DB = Spiner::DataBox; - constexpr int RANK = 2; - const int Npt = 32; - const int Nsample = Npt * 16; + constexpr int RANK = 2; + const int Npt = 32; + const int Nsample = Npt * 16; - const double xmin = 1; - const double xmax = 10; + const double xmin = 1; + const double xmax = 10; - DB db(Spiner::AllocationTarget::Device, Npt, Npt); - for (int i = 0; i < RANK; i++) { - db.setRange(i, xmin, xmax, Npt); - } - - portableFor( - "fill databox", 0, Npt, 0, Npt, - PORTABLE_LAMBDA(const int ix, const int iy) { - RG1D grid(xmin, xmax, Npt); - const double x = grid.x(ix); - const double y = grid.x(iy); - db(ix, iy) = log_lin_func(x, y); - }); + DB db(Spiner::AllocationTarget::Device, Npt, Npt); + for (int i = 0; i < RANK; i++) { + db.setRange(i, xmin, xmax, Npt); + } - double error = 0; - portableReduce( - "interpolate databox", 0, Nsample, 0, Nsample, - PORTABLE_LAMBDA(const int ix, const int iy, Real &accumulate) { - RG1D grid(xmin, xmax, Nsample); - const double x = grid.x(ix); - const double y = grid.x(iy); - const double f_true = log_lin_func(x, y); - const double difference = db.interpToReal(x, y) - f_true; - accumulate += (difference * difference); - }, - error); - REQUIRE(error <= EPSTEST); + portableFor( + "fill databox", 0, Npt, 0, Npt, + PORTABLE_LAMBDA(const int ix, const int iy) { + RG1D grid(xmin, xmax, Npt); + const double x = grid.x(ix); + const double y = grid.x(iy); + db(ix, iy) = log_lin_func(x, y); + }); + + double error = 0; + portableReduce( + "interpolate databox", 0, Nsample, 0, Nsample, + PORTABLE_LAMBDA(const int ix, const int iy, Real &accumulate) { + RG1D grid(xmin, xmax, Nsample); + const double x = grid.x(ix); + const double y = grid.x(iy); + const double f_true = log_lin_func(x, y); + const double difference = db.interpToReal(x, y) - f_true; + accumulate += (difference * difference); + }, + error); + REQUIRE(error <= EPSTEST); } TEST_CASE("DataBox Interpolation with log-log transformations", "[DataBox]") { - using Transform = Spiner::TransformLogarithmic; - using RG1D = Spiner::RegularGrid1D; - using DB = Spiner::DataBox; + using Transform = Spiner::TransformLogarithmic; + using RG1D = Spiner::RegularGrid1D; + using DB = Spiner::DataBox; - constexpr int RANK = 2; - const int Npt = 32; - const int Nsample = Npt * 16; + constexpr int RANK = 2; + const int Npt = 32; + const int Nsample = Npt * 16; - const double xmin = 1; - const double xmax = 10; + const double xmin = 1; + const double xmax = 10; - DB db(Spiner::AllocationTarget::Device, Npt, Npt); - for (int i = 0; i < RANK; i++) { - db.setRange(i, xmin, xmax, Npt); - } - - portableFor( - "fill databox", 0, Npt, 0, Npt, - PORTABLE_LAMBDA(const int ix, const int iy) { - RG1D grid(xmin, xmax, Npt); - const double x = grid.x(ix); - const double y = grid.x(iy); - db.set_data_value(log_log_func(x,y), ix, iy); - }); + DB db(Spiner::AllocationTarget::Device, Npt, Npt); + for (int i = 0; i < RANK; i++) { + db.setRange(i, xmin, xmax, Npt); + } - double error = 0; - portableReduce( - "interpolate databox", 0, Nsample, 0, Nsample, - PORTABLE_LAMBDA(const int ix, const int iy, Real &accumulate) { - RG1D grid(xmin, xmax, Nsample); - const double x = grid.x(ix); - const double y = grid.x(iy); - const double f_true = log_log_func(x, y); - const double difference = db.interpToReal(x, y) - f_true; - accumulate += (difference * difference); - }, - error); - REQUIRE(error <= EPSTEST); + portableFor( + "fill databox", 0, Npt, 0, Npt, + PORTABLE_LAMBDA(const int ix, const int iy) { + RG1D grid(xmin, xmax, Npt); + const double x = grid.x(ix); + const double y = grid.x(iy); + db.set_data_value(log_log_func(x, y), ix, iy); + }); + + double error = 0; + portableReduce( + "interpolate databox", 0, Nsample, 0, Nsample, + PORTABLE_LAMBDA(const int ix, const int iy, Real &accumulate) { + RG1D grid(xmin, xmax, Nsample); + const double x = grid.x(ix); + const double y = grid.x(iy); + const double f_true = log_log_func(x, y); + const double difference = db.interpToReal(x, y) - f_true; + accumulate += (difference * difference); + }, + error); + REQUIRE(error <= EPSTEST); } TEST_CASE("DataBox Interpolation with piecewise grids", diff --git a/test/transformations.cpp b/test/transformations.cpp index b7a05133c..83a389a8d 100644 --- a/test/transformations.cpp +++ b/test/transformations.cpp @@ -24,82 +24,79 @@ using Catch::Matchers::WithinRel; using Catch::Matchers::WithinULP; namespace { - PORTABLE_FUNCTION double relerr(const double z, const double zz) { - if (z == 0) { - return std::abs(zz - z); - } else { - return std::abs(double(1) - zz/z); - } - } -}; +PORTABLE_FUNCTION double relerr(const double z, const double zz) { + if (z == 0) { + return std::abs(zz - z); + } else { + return std::abs(double(1) - zz / z); + } +} +}; // namespace TEST_CASE("transformation: linear", "[transformations]") { - // This one is almost too simple to meaningfully test, but we can at least - // ensure that it compiles and does the trivially-obvious things. - using Transform = Spiner::TransformLinear; + // This one is almost too simple to meaningfully test, but we can at least + // ensure that it compiles and does the trivially-obvious things. + using Transform = Spiner::TransformLinear; - // Test on CPU - const double x0 = 3.14159; - CHECK(Transform::forward(x0) == x0); - CHECK(Transform::reverse(x0) == x0); - const float x1 = 2.71828; - CHECK(Transform::forward(x1) == x1); - CHECK(Transform::reverse(x1) == x1); + // Test on CPU + const double x0 = 3.14159; + CHECK(Transform::forward(x0) == x0); + CHECK(Transform::reverse(x0) == x0); + const float x1 = 2.71828; + CHECK(Transform::forward(x1) == x1); + CHECK(Transform::reverse(x1) == x1); - // Test on GPU (or CPU if no GPU available) - const int num_threads = 10; - double accum = 0; - portableReduce( - "parallel_reduce", 0, num_threads, - PORTABLE_LAMBDA(const int n, double &reduce) - { - const double x = static_cast(n); - const double y = x; - const double yy = Transform::forward(x); - const double xx = Transform::reverse(y); - reduce += 0.5 * (relerr(x, xx) + relerr(y, yy)); - }, - accum); - CHECK(accum / num_threads == 0); + // Test on GPU (or CPU if no GPU available) + const int num_threads = 10; + double accum = 0; + portableReduce( + "parallel_reduce", 0, num_threads, + PORTABLE_LAMBDA(const int n, double &reduce) { + const double x = static_cast(n); + const double y = x; + const double yy = Transform::forward(x); + const double xx = Transform::reverse(y); + reduce += 0.5 * (relerr(x, xx) + relerr(y, yy)); + }, + accum); + CHECK(accum / num_threads == 0); } TEST_CASE("transformation: logarithmic", "[transformations]") { - using Transform = Spiner::TransformLogarithmic; + using Transform = Spiner::TransformLogarithmic; - // Test on CPU - auto matcher = [](const double d) { - //return WithinULP(d, 500); - return WithinRel(d, 6.0e-12); - }; - // Scan across most of the range - for (int p = -300; p <= 300; ++p) { - const double x = std::pow(double(10), p); - const double y = std::log(x); - // Basic checks - CHECK_THAT(Transform::forward(x), matcher(y)); - CHECK_THAT(Transform::reverse(y), matcher(x)); - // Round trip - CHECK_THAT(Transform::reverse(Transform::forward(x)), matcher(x)); - } - // Special value - CHECK(std::isfinite(Transform::forward(0))); - CHECK(Transform::reverse(Transform::forward(0)) == 0); + // Test on CPU + auto matcher = [](const double d) { + // return WithinULP(d, 500); + return WithinRel(d, 6.0e-12); + }; + // Scan across most of the range + for (int p = -300; p <= 300; ++p) { + const double x = std::pow(double(10), p); + const double y = std::log(x); + // Basic checks + CHECK_THAT(Transform::forward(x), matcher(y)); + CHECK_THAT(Transform::reverse(y), matcher(x)); + // Round trip + CHECK_THAT(Transform::reverse(Transform::forward(x)), matcher(x)); + } + // Special value + CHECK(std::isfinite(Transform::forward(0))); + CHECK(Transform::reverse(Transform::forward(0)) == 0); - // Test on GPU (or CPU if no GPU available) - const int num_threads = 101; - double accum = 0; - portableReduce( - "parallel_reduce", 0, num_threads, - PORTABLE_LAMBDA(const int n, double &reduce) - { - const int p = n - num_threads / 2; - const double x = std::pow(double(10.0), p); - const double y = std::log(x); - const double yy = Transform::forward(x); - const double xx = Transform::reverse(y); - reduce += 0.5 * (relerr(x, xx) + relerr(y, yy)); - }, - accum); - CHECK(accum / num_threads <= 6.0e-12); + // Test on GPU (or CPU if no GPU available) + const int num_threads = 101; + double accum = 0; + portableReduce( + "parallel_reduce", 0, num_threads, + PORTABLE_LAMBDA(const int n, double &reduce) { + const int p = n - num_threads / 2; + const double x = std::pow(double(10.0), p); + const double y = std::log(x); + const double yy = Transform::forward(x); + const double xx = Transform::reverse(y); + reduce += 0.5 * (relerr(x, xx) + relerr(y, yy)); + }, + accum); + CHECK(accum / num_threads <= 6.0e-12); } - From 28bcf26df64f9bf6853d3bd89c58d5096c12755b Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 30 Sep 2024 12:01:55 -0600 Subject: [PATCH 33/48] missing header --- spiner/transformations.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index 03139d41b..5cd77daf1 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -17,6 +17,7 @@ #include "ports-of-call/portability.hpp" +#include #include namespace Spiner { From 51d11c88ee68fbb185ba62b36dca9d01b5f8641c Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" <43178897+BrendanKKrueger@users.noreply.github.com> Date: Wed, 2 Oct 2024 16:43:56 -0600 Subject: [PATCH 34/48] Update spiner/transformations.hpp Co-authored-by: Jonah Miller --- spiner/transformations.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index 5cd77daf1..bf462f5aa 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -44,7 +44,7 @@ struct TransformLinear { struct TransformLogarithmic { template PORTABLE_INLINE_FUNCTION static T forward(const T x) { - return std::log(x + std::numeric_limits::denorm_min()); + return std::log(std::abs(x) + std::numeric_limits::denorm_min()); } template PORTABLE_INLINE_FUNCTION static T reverse(const T u) { From df8f052dbbaa014bbb8cdc65c838126ebd3ac5ae Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 2 Oct 2024 16:47:01 -0600 Subject: [PATCH 35/48] force inline --- spiner/transformations.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index bf462f5aa..7946aec2f 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -31,11 +31,11 @@ namespace Spiner { // linear transformation (aka no-op): y = x struct TransformLinear { template - PORTABLE_INLINE_FUNCTION static T forward(const T x) { + PORTABLE_FORCEINLINE_FUNCTION static T forward(const T x) { return x; } template - PORTABLE_INLINE_FUNCTION static T reverse(const T u) { + PORTABLE_FORCEINLINE_FUNCTION static T reverse(const T u) { return u; } }; @@ -43,11 +43,11 @@ struct TransformLinear { // logarithmic transformation: y = log(x + small) struct TransformLogarithmic { template - PORTABLE_INLINE_FUNCTION static T forward(const T x) { - return std::log(std::abs(x) + std::numeric_limits::denorm_min()); + PORTABLE_FORCEINLINE_FUNCTION static T forward(const T x) { + return std::log(std::abs(x) + std::numeric_limits::denorm_min()); } template - PORTABLE_INLINE_FUNCTION static T reverse(const T u) { + PORTABLE_FORCEINLINE_FUNCTION static T reverse(const T u) { return std::exp(u) - std::numeric_limits::denorm_min(); } }; From 2f3915028d29c80fba6b110c394478b327905050 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 2 Oct 2024 17:01:03 -0600 Subject: [PATCH 36/48] constexpr --- spiner/transformations.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index 7946aec2f..5a98cc818 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -31,11 +31,11 @@ namespace Spiner { // linear transformation (aka no-op): y = x struct TransformLinear { template - PORTABLE_FORCEINLINE_FUNCTION static T forward(const T x) { + PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) { return x; } template - PORTABLE_FORCEINLINE_FUNCTION static T reverse(const T u) { + PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) { return u; } }; @@ -43,11 +43,11 @@ struct TransformLinear { // logarithmic transformation: y = log(x + small) struct TransformLogarithmic { template - PORTABLE_FORCEINLINE_FUNCTION static T forward(const T x) { + PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) { return std::log(std::abs(x) + std::numeric_limits::denorm_min()); } template - PORTABLE_FORCEINLINE_FUNCTION static T reverse(const T u) { + PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) { return std::exp(u) - std::numeric_limits::denorm_min(); } }; From 8d76c554444fbda78d85efeed0fbfa500dd06a6d Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Tue, 8 Oct 2024 10:50:36 -0600 Subject: [PATCH 37/48] Documentation --- doc/sphinx/src/databox.rst | 45 +++++++++++++++++++++++++++++++++++++- 1 file changed, 44 insertions(+), 1 deletion(-) diff --git a/doc/sphinx/src/databox.rst b/doc/sphinx/src/databox.rst index 624e187cf..439b615a5 100644 --- a/doc/sphinx/src/databox.rst +++ b/doc/sphinx/src/databox.rst @@ -14,7 +14,13 @@ To use databox, simply include the relevant header: #include -``DatBox`` is templated on underyling data type, which defaults to the +``DataBox`` Templates +^^^^^^^^^^^^^^^^^^^^^ + +Underlying Data Type +-------------------- + +``DataBox`` is templated on the underyling data type, which defaults to the ``Real`` type provided by ``ports-of-call``. (This is usually a ``double``.) @@ -30,6 +36,9 @@ single type, you may wish to declare a type alias such as: using DataBox = Spiner::DataBox +Interpolation Gridding +---------------------- + Spiner is also templated on how the interpolation gridding works. This template parameter is called ``Grid_t``. The available options at this time are: @@ -47,6 +56,40 @@ as, for example: More detail on the interpolation gridding is available below and in the interpolation section. +Transformations +--------------- + +Spiner performs linear interpolation of the data, but can apply transformations +to the data in order to enable interpolation methods like log-log, log-lin, +etc. Spiner will perform linear interpolation of the transformed variables. + +When constructing the gridding, ``Spiner::RegularGrid1D`` and +``Spiner::PiecewiseGrid1D`` are templated on the independent variable +transformation. Note that all independent variables will use the same +transformation. ``DataBox`` is templated on the dependent variable +transformation. In all cases, the default is a linear "transformation" (no +transformation to the data), so if you do not specify the transformation then +you get linear interpolation. The available transformations are in +``spiner/transformations.hpp`` and users can write custom transformations +following the examples there. + +For example, if the user wants to transform the independent variables with a +logarithm and the dependent variable with a custom arctangent transformation, +the declaration of the ``DataBox`` might look like + +.. code-block:: cpp + + using DataBox = Spiner::DataBox, CustomArctanTransform>; + +This would result in interpolation according to the equation + +.. math:: + + \arctan(y) = b + \sum m_i \log(x_i) + +where :math:`x_i` are the independent variables. + .. note:: In C++17 and later, you can also get the default type specialization by simply omitting the template arguments. From 09c3aedc7a6f9514280863f12b78ab6ed01a8eef Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Tue, 8 Oct 2024 10:54:32 -0600 Subject: [PATCH 38/48] documentation again --- doc/sphinx/src/databox.rst | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/src/databox.rst b/doc/sphinx/src/databox.rst index 439b615a5..fd4114161 100644 --- a/doc/sphinx/src/databox.rst +++ b/doc/sphinx/src/databox.rst @@ -86,9 +86,15 @@ This would result in interpolation according to the equation .. math:: - \arctan(y) = b + \sum m_i \log(x_i) + v &= \arctan(y) -where :math:`x_i` are the independent variables. + u_i &= \log(x_i) + + v &= b + \sum m_i u_i + +where :math:`x_i` are the independent variables. By convention, the Spiner +documentation uses :math:`u_i` as the transformed independent variables, and +:math:`v` as the transformed dependent variable. .. note:: In C++17 and later, you can also get the default type specialization From 45fe193c4732d0180526197994d45336dbb183e3 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 9 Oct 2024 11:46:00 -0600 Subject: [PATCH 39/48] Change logarithmic transform's epsilon value. --- spiner/transformations.hpp | 33 +++++++++++++++++++++++++++++++-- test/transformations.cpp | 19 +++++++++++++------ 2 files changed, 44 insertions(+), 8 deletions(-) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index 5a98cc818..be4a8241c 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -44,11 +44,40 @@ struct TransformLinear { struct TransformLogarithmic { template PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) { - return std::log(std::abs(x) + std::numeric_limits::denorm_min()); + constexpr T eps = eps_f(); + return std::log(std::abs(x) + eps); } template PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) { - return std::exp(u) - std::numeric_limits::denorm_min(); + constexpr T eps = eps_r(); + return std::exp(u) - eps; + } +private: + // When possible, we use asymetric epsilon values to ensure that + // reverse(forward(0)) is exact. In general, a performant calculation is + // more important than getting this value exactly correct, so we require that + // our epsilon values be constexpr. + template + PORTABLE_FORCEINLINE_FUNCTION static constexpr T eps_f() { + return std::numeric_limits::min(); + } + template + PORTABLE_FORCEINLINE_FUNCTION static constexpr T eps_r() { + // If std::exp and std::log were constexpr, we could explicitly calculate + // the right epsilon value as a constexpr value: + // return std::exp(forward(static_cast(0))); + // Unfortunately, we have to wait for C++26 for constexpr math. We + // hard-code certain known values, but default to a symmetric epsilon if + // that's the best that we can do. + if constexpr (std::is_same_v, float> && + std::numeric_limits::is_iec559) { + return 1.175490707446280263444352135525329e-38f; + } else if constexpr (std::is_same_v, double> && + std::numeric_limits::is_iec559) { + return 2.225073858507262647230317031903882e-308; + } else { + return eps_f(); + } } }; diff --git a/test/transformations.cpp b/test/transformations.cpp index 83a389a8d..c95259d18 100644 --- a/test/transformations.cpp +++ b/test/transformations.cpp @@ -68,21 +68,28 @@ TEST_CASE("transformation: logarithmic", "[transformations]") { // Test on CPU auto matcher = [](const double d) { // return WithinULP(d, 500); - return WithinRel(d, 6.0e-12); + return WithinRel(d, 1.0e-12); }; // Scan across most of the range - for (int p = -300; p <= 300; ++p) { + for (int p = -307; p <= 307; ++p) { const double x = std::pow(double(10), p); const double y = std::log(x); // Basic checks - CHECK_THAT(Transform::forward(x), matcher(y)); - CHECK_THAT(Transform::reverse(y), matcher(x)); + // -- The "exact" calculation (y = log(x)) won't match the transformation + // (y = log(x + epsilon) for very small values of x, so skip checks. + if (p >= -295) { + CHECK_THAT(Transform::forward(x), matcher(y)); + CHECK_THAT(Transform::reverse(y), matcher(x)); + } // Round trip CHECK_THAT(Transform::reverse(Transform::forward(x)), matcher(x)); } // Special value - CHECK(std::isfinite(Transform::forward(0))); - CHECK(Transform::reverse(Transform::forward(0)) == 0); + // -- Transform::forward(0) will infer the type to be an integer, and you + // will get a WILDLY incorrect answer for the round trip. + CHECK(std::isfinite(Transform::forward(0.0))); + CHECK(Transform::reverse(Transform::forward(static_cast(0))) == 0); + CHECK(Transform::reverse(Transform::forward(static_cast(0))) == 0); // Test on GPU (or CPU if no GPU available) const int num_threads = 101; From 8578f895c40a44b6c72ac60da96a19aebca56649 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 9 Oct 2024 11:52:02 -0600 Subject: [PATCH 40/48] format --- spiner/transformations.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/spiner/transformations.hpp b/spiner/transformations.hpp index be4a8241c..f491c3ba0 100644 --- a/spiner/transformations.hpp +++ b/spiner/transformations.hpp @@ -52,7 +52,8 @@ struct TransformLogarithmic { constexpr T eps = eps_r(); return std::exp(u) - eps; } -private: + + private: // When possible, we use asymetric epsilon values to ensure that // reverse(forward(0)) is exact. In general, a performant calculation is // more important than getting this value exactly correct, so we require that From 5942d20bf58630abe5415b8d580a95cad9f7017f Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Mon, 14 Oct 2024 13:00:54 -0600 Subject: [PATCH 41/48] Switch from SFINAE to static_assert so that we can give a better error message. --- spiner/databox.hpp | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index eae25e410..ce2c340d4 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -165,18 +165,15 @@ class DataBox { // Index operators // examle calls: // T x = db(n4, n3, n2, n1); - template ::value>::type * = nullptr> - PORTABLE_INLINE_FUNCTION T &operator()(Args... args) { - return dataView_(std::forward(args)...); - } - template ::value>::type * = nullptr> + template PORTABLE_INLINE_FUNCTION T &operator()(Args... args) const { + static_assert(std::is_same::value, + "DataBox::operator() is only available if the independent " + "variable data transform is Spiner::TransformLinear, because " + "otherwise accessing the underlying data is not doing what you " + "expect. Use DataBox::get_data_value and DataBox::set_data_value " + "instead, because those correctly account for the independent " + "value data transformation."); return dataView_(std::forward(args)...); } From 624d2443d4e740a3b11bb3ef8ee3e5b39e71fdee Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 23 Oct 2024 08:50:11 -0600 Subject: [PATCH 42/48] x is the 'ground truth' representation, so fix the bounds. --- spiner/regular_grid_1d.hpp | 45 +++++++++++++++++++++++++++----------- 1 file changed, 32 insertions(+), 13 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index 032a59601..43922d50e 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -58,9 +58,17 @@ class RegularGrid1D { PORTABLE_INLINE_FUNCTION RegularGrid1D() : 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"); + : xmin_(xmin) + , xmax_(xmax) + , umin_(Transform::forward(xmin)) + , umax_(Transform::forward(xmax)) + , du_((umax_ - umin_) / static_cast(N - 1)) + , inv_du_(1 / du_) + , N_(N) + { + // A transform could be monotonically decreasing, so there's no guarantee + // that umin_ < umax_ + PORTABLE_ALWAYS_REQUIRE(xmin_ < xmax_ && N_ > 0, "Valid grid"); } // Forces x in the interval @@ -72,6 +80,7 @@ class RegularGrid1D { return ix; } + // TODO: make these private to hide u as an internal detail // 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 { @@ -82,13 +91,13 @@ class RegularGrid1D { PORTABLE_INLINE_FUNCTION T x(const int i) const { return Transform::reverse(u(i)); } + // TODO: Delete index since internally we only use index_u? PORTABLE_INLINE_FUNCTION int index(const T x) const { 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 { + PORTABLE_INLINE_FUNCTION void weights(const T &x, int &ix, weights_t &w) const { const T u = Transform::forward(x); ix = index_u(u); const auto floor = static_cast(ix) * du_ + umin_; @@ -107,22 +116,31 @@ class RegularGrid1D { // utitilies PORTABLE_INLINE_FUNCTION bool - operator==(const RegularGrid1D &other) const { - return (umin_ == other.umin_ && umax_ == other.umax_ && du_ == other.du_ && - inv_du_ == other.inv_du_ && N_ == other.N_); + operator==(const RegularGrid1D &other) const { + return (umin_ == other.umin_ && + umax_ == other.umax_ && + du_ == other.du_ && + N_ == other.N_); } PORTABLE_INLINE_FUNCTION bool - operator!=(const RegularGrid1D &other) const { + operator!=(const RegularGrid1D &other) const { return !(*this == other); } + // TODO: umin, umax should be private 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 T min() const { return xmin_; } + PORTABLE_INLINE_FUNCTION T max() const { return xmax_; } PORTABLE_INLINE_FUNCTION size_t nPoints() const { return N_; } PORTABLE_INLINE_FUNCTION bool isnan() const { - return (std::isnan(umin_) || std::isnan(umax_) || std::isnan(du_) || - std::isnan(inv_du_) || std::isnan((T)N_)); + return (std::isnan(xmin_) || + std::isnan(xmax_) || + 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(); } @@ -166,6 +184,7 @@ class RegularGrid1D { #endif private: + T xmin_, xmax_; T umin_, umax_; T du_, inv_du_; size_t N_; From 0486a9d9b921a67de3298c3fc5b8af9d24e90703 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 23 Oct 2024 08:53:33 -0600 Subject: [PATCH 43/48] Cleaning up things that should (probably?) be private. --- spiner/regular_grid_1d.hpp | 41 +++++++++++++++++--------------------- 1 file changed, 18 insertions(+), 23 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index 43922d50e..cbf89b3e5 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -71,30 +71,13 @@ class RegularGrid1D { PORTABLE_ALWAYS_REQUIRE(xmin_ < xmax_ && N_ > 0, "Valid grid"); } - // Forces x in the interval - PORTABLE_INLINE_FUNCTION int bound(int ix) const { -#ifndef SPINER_DISABLE_BOUNDS_CHECKS - if (ix < 0) ix = 0; - if (ix >= (int)N_ - 1) ix = (int)N_ - 2; // Ensures ix+1 exists -#endif - return ix; - } - - // TODO: make these private to hide u as an internal detail - // 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_)); - } - + // TODO: min() and x(0) won't necessarily match. + // max() and x(nPoints-1) won't necessarily match. + // Should we do anything about this? // Translate between x coordinate and index PORTABLE_INLINE_FUNCTION T x(const int i) const { return Transform::reverse(u(i)); } - // TODO: Delete index since internally we only use index_u? - PORTABLE_INLINE_FUNCTION int index(const T x) const { - 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 { @@ -126,9 +109,6 @@ class RegularGrid1D { operator!=(const RegularGrid1D &other) const { return !(*this == other); } - // TODO: umin, umax should be private - PORTABLE_INLINE_FUNCTION T umin() const { return umin_; } - PORTABLE_INLINE_FUNCTION T umax() const { return umax_; } PORTABLE_INLINE_FUNCTION T min() const { return xmin_; } PORTABLE_INLINE_FUNCTION T max() const { return xmax_; } @@ -184,6 +164,21 @@ class RegularGrid1D { #endif private: + // Forces x in the interval + PORTABLE_INLINE_FUNCTION int bound(int ix) const { +#ifndef SPINER_DISABLE_BOUNDS_CHECKS + if (ix < 0) ix = 0; + if (ix >= (int)N_ - 1) ix = (int)N_ - 2; // Ensures ix+1 exists +#endif + return ix; + } + + // 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_)); + } + T xmin_, xmax_; T umin_, umax_; T du_, inv_du_; From aba2fd46b90789e4f527347a1a603fc7c9a7916a Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 23 Oct 2024 09:05:12 -0600 Subject: [PATCH 44/48] More cleanup --- spiner/regular_grid_1d.hpp | 27 ++++++++++++++++----------- test/regular_grid_1d.cpp | 4 ---- 2 files changed, 16 insertions(+), 15 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index cbf89b3e5..246da7a26 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -49,7 +49,7 @@ template ::value, bool>::type = true> class RegularGrid1D { - public: +public: using ValueType = T; static constexpr T rNaN = std::numeric_limits::signaling_NaN(); static constexpr int iNaN = std::numeric_limits::signaling_NaN(); @@ -71,14 +71,6 @@ class RegularGrid1D { PORTABLE_ALWAYS_REQUIRE(xmin_ < xmax_ && N_ > 0, "Valid grid"); } - // TODO: min() and x(0) won't necessarily match. - // max() and x(nPoints-1) won't necessarily match. - // Should we do anything about this? - // Translate between x coordinate and index - PORTABLE_INLINE_FUNCTION T x(const int i) const { - return Transform::reverse(u(i)); - } - // Returns closest index and weights for interpolation PORTABLE_INLINE_FUNCTION void weights(const T &x, int &ix, weights_t &w) const { const T u = Transform::forward(x); @@ -97,7 +89,7 @@ class RegularGrid1D { return w[0] * A(ix) + w[1] * A(ix + 1); } - // utitilies + // (in)equality comparison PORTABLE_INLINE_FUNCTION bool operator==(const RegularGrid1D &other) const { return (umin_ == other.umin_ && @@ -110,6 +102,7 @@ class RegularGrid1D { return !(*this == other); } + // queries PORTABLE_INLINE_FUNCTION T min() const { return xmin_; } PORTABLE_INLINE_FUNCTION T max() const { return xmax_; } PORTABLE_INLINE_FUNCTION size_t nPoints() const { return N_; } @@ -123,7 +116,19 @@ class RegularGrid1D { std::isnan((T)N_)); } PORTABLE_INLINE_FUNCTION bool isWellFormed() const { return !isnan(); } + // TODO: min() and x(0) won't necessarily match. + // max() and x(nPoints-1) won't necessarily match. + // Should we do anything about this? + // The easiest fix: if i == 0 return xmin_ and similar for xmax_ + // 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 index_u(Transform::forward(x)); + } + // HDF #ifdef SPINER_USE_HDF inline herr_t saveHDF(hid_t loc, const std::string &name) const { static_assert( @@ -163,7 +168,7 @@ class RegularGrid1D { } #endif - private: +private: // Forces x in the interval PORTABLE_INLINE_FUNCTION int bound(int ix) const { #ifndef SPINER_DISABLE_BOUNDS_CHECKS diff --git a/test/regular_grid_1d.cpp b/test/regular_grid_1d.cpp index 5beaf5853..3d2e3cbb2 100644 --- a/test/regular_grid_1d.cpp +++ b/test/regular_grid_1d.cpp @@ -51,16 +51,12 @@ TEST_CASE("RegularGrid1D with transformations", // 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)); } From 82fe518caf679a57a90e4e47ec13125a214db7f7 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 23 Oct 2024 10:23:13 -0600 Subject: [PATCH 45/48] Expanded tests for RegularGrid1D. --- test/regular_grid_1d.cpp | 135 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 134 insertions(+), 1 deletion(-) diff --git a/test/regular_grid_1d.cpp b/test/regular_grid_1d.cpp index 3d2e3cbb2..3839094c1 100644 --- a/test/regular_grid_1d.cpp +++ b/test/regular_grid_1d.cpp @@ -22,6 +22,54 @@ using Catch::Matchers::WithinRel; +// test transform: reverse(forward(xlo)) and reverse(forward(xhi)) don't map +// 100% correctly and end up inside the bounds of [xlo, xhi], so long as xlo +// and xhi bracket xref = 0.5 +struct TxNarrow { + template + PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) { + constexpr T xref = 0.5; + constexpr T eps = 10 * std::numeric_limits::epsilon(); + return (static_cast(1) - eps) * (x - xref) + xref; + } + template + PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) { + constexpr T uref = 0.5; + constexpr T eps = 10 * std::numeric_limits::epsilon(); + return (static_cast(1) - eps) * (u - uref) + uref; + } +}; + +// test transform: reverse(forward(xlo)) and reverse(forward(xhi)) don't map +// 100% correctly and end up outside the bounds of [xlo, xhi], so long as xlo +// and xhi bracket xref = 0.5 +struct TxExpand { + template + PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) { + constexpr T xref = 0.5; + constexpr T eps = 10 * std::numeric_limits::epsilon(); + return (static_cast(1) + eps) * (x - xref) + xref; + } + template + PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) { + constexpr T uref = 0.5; + constexpr T eps = 10 * std::numeric_limits::epsilon(); + return (static_cast(1) + eps) * (u - uref) + uref; + } +}; + +// test transform: monotonically decreasing instead of increasing +struct TxDecrease { + template + PORTABLE_FORCEINLINE_FUNCTION static constexpr T forward(const T x) { + return 1 - x; + } + template + PORTABLE_FORCEINLINE_FUNCTION static constexpr T reverse(const T u) { + return 1 - u; + } +}; + TEST_CASE("RegularGrid1D", "[RegularGrid1D]") { SECTION("A regular grid 1d emits appropriate metadata") { constexpr Real min = -1; @@ -34,7 +82,7 @@ TEST_CASE("RegularGrid1D", "[RegularGrid1D]") { } } -TEST_CASE("RegularGrid1D with transformations", +TEST_CASE("RegularGrid1D with production transformations", "[RegularGrid1D][transformations]") { constexpr double min = 1; constexpr double max = 1024; @@ -82,3 +130,88 @@ TEST_CASE("RegularGrid1D with transformations", CHECK_THAT(wlog.second, WithinRel(0.5)); } } + +TEST_CASE("RegularGrid1D with test transformations", + "[RegularGrid1D][transformations]") { + constexpr double min = 0.0; + constexpr double max = 1.0; + constexpr size_t N = 101; + Spiner::RegularGrid1D gn(min, max, N); + Spiner::RegularGrid1D ge(min, max, N); + Spiner::RegularGrid1D gd(min, max, N); + + // Basic tests (narrow) + REQUIRE(gn.min() == min); + REQUIRE(gn.max() == max); + REQUIRE(gn.nPoints() == N); + // TODO: Do we want the bounds to be exact? + CHECK(gn.x(0) == min); + CHECK(gn.x(N-1) == max); + + // Basic tests (expand) + REQUIRE(ge.min() == min); + REQUIRE(ge.max() == max); + REQUIRE(ge.nPoints() == N); + // TODO: Do we want the bounds to be exact? + CHECK(ge.x(0) == min); + CHECK(ge.x(N-1) == max); + + // Basic tests (decrease) + REQUIRE(gd.min() == min); + REQUIRE(gd.max() == max); + REQUIRE(gd.nPoints() == N); + // TODO: Do we want the bounds to be exact? + CHECK(gd.x(0) == min); + CHECK(gd.x(N-1) == max); + + // Check all fixed points (narrow) + for (std::size_t n = 0; n < N; ++n) { + const double xx = 0.01 * double(n); + CHECK_THAT(gn.x(n), WithinRel(xx, 1.0e-12)); + } + + // Check all fixed points (expand) + for (std::size_t n = 0; n < N; ++n) { + const double xx = 0.01 * double(n); + CHECK_THAT(ge.x(n), WithinRel(xx, 1.0e-12)); + } + + // Check all fixed points (decrease) + for (std::size_t n = 0; n < N; ++n) { + const double xx = 0.01 * double(n); + CHECK_THAT(gd.x(n), WithinRel(xx, 1.0e-12)); + } + + // Check weights (narrow) + for (std::size_t n = 1; n < N; ++n) { + const double x = 0.5 * (gn.x(n - 1) + gn.x(n)); + int index; + Spiner::weights_t weights; + gn.weights(x, index, weights); + CHECK(index == n - 1); + CHECK_THAT(weights.first, WithinRel(0.5, 1.0e-12)); + CHECK_THAT(weights.second, WithinRel(0.5, 1.0e-12)); + } + + // Check weights (expand) + for (std::size_t n = 1; n < N; ++n) { + const double x = 0.5 * (ge.x(n - 1) + ge.x(n)); + int index; + Spiner::weights_t weights; + ge.weights(x, index, weights); + CHECK(index == n - 1); + CHECK_THAT(weights.first, WithinRel(0.5, 1.0e-12)); + CHECK_THAT(weights.second, WithinRel(0.5, 1.0e-12)); + } + + // Check weights (decrease) + for (std::size_t n = 1; n < N; ++n) { + const double x = 0.5 * (gd.x(n - 1) + gd.x(n)); + int index; + Spiner::weights_t weights; + gd.weights(x, index, weights); + CHECK(index == n - 1); + CHECK_THAT(weights.first, WithinRel(0.5, 1.0e-12)); + CHECK_THAT(weights.second, WithinRel(0.5, 1.0e-12)); + } +} From ab02e9fcc68176a197b632e4838a5a3ebaae8489 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 23 Oct 2024 10:41:25 -0600 Subject: [PATCH 46/48] Expand TODO comment. --- spiner/regular_grid_1d.hpp | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index 246da7a26..ed4635bdc 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -118,8 +118,18 @@ class RegularGrid1D { PORTABLE_INLINE_FUNCTION bool isWellFormed() const { return !isnan(); } // TODO: min() and x(0) won't necessarily match. // max() and x(nPoints-1) won't necessarily match. - // Should we do anything about this? - // The easiest fix: if i == 0 return xmin_ and similar for xmax_ + // Possible fixes: + // -- If (i == 0) return xmin_. Introduces two "if" statements every + // time you look up at x value. + // -- Allow min() and x(0) to differ. This could be confusing and + // cause issues in corner cases. + // -- Allow min() to be different from xmin passed in by the user, + // which may be unexpected (and possibly undesirable) behavior. + // -- Apply an additional linear transformation calibrated so that the + // bounds are exact. This may muck with any calibrations of the + // transform itself. + // -- Carry an additional array of x values to be used here. This + // introduces a lot of memory, when no memory was used before. // Translate between x coordinate and index PORTABLE_INLINE_FUNCTION T x(const int i) const { return Transform::reverse(u(i)); From 1955937b098834e5cff136a2e0b74ff320296bd2 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 23 Oct 2024 10:46:44 -0600 Subject: [PATCH 47/48] format --- spiner/databox.hpp | 15 ++++++++------- spiner/regular_grid_1d.hpp | 29 ++++++++++------------------- test/regular_grid_1d.cpp | 6 +++--- 3 files changed, 21 insertions(+), 29 deletions(-) diff --git a/spiner/databox.hpp b/spiner/databox.hpp index ce2c340d4..caa752342 100644 --- a/spiner/databox.hpp +++ b/spiner/databox.hpp @@ -167,13 +167,14 @@ class DataBox { // T x = db(n4, n3, n2, n1); template PORTABLE_INLINE_FUNCTION T &operator()(Args... args) const { - static_assert(std::is_same::value, - "DataBox::operator() is only available if the independent " - "variable data transform is Spiner::TransformLinear, because " - "otherwise accessing the underlying data is not doing what you " - "expect. Use DataBox::get_data_value and DataBox::set_data_value " - "instead, because those correctly account for the independent " - "value data transformation."); + static_assert( + std::is_same::value, + "DataBox::operator() is only available if the independent variable " + "data transform is Spiner::TransformLinear, because otherwise " + "accessing the underlying data is not doing what you expect. Use " + "DataBox::get_data_value and DataBox::set_data_value instead, because " + "those correctly account for the independent value data " + "transformation."); return dataView_(std::forward(args)...); } diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index ed4635bdc..b7b371d70 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -49,7 +49,7 @@ template ::value, bool>::type = true> class RegularGrid1D { -public: + public: using ValueType = T; static constexpr T rNaN = std::numeric_limits::signaling_NaN(); static constexpr int iNaN = std::numeric_limits::signaling_NaN(); @@ -58,13 +58,9 @@ class RegularGrid1D { PORTABLE_INLINE_FUNCTION RegularGrid1D() : umin_(rNaN), umax_(rNaN), du_(rNaN), inv_du_(rNaN), N_(iNaN) {} PORTABLE_INLINE_FUNCTION RegularGrid1D(T xmin, T xmax, size_t N) - : xmin_(xmin) - , xmax_(xmax) - , umin_(Transform::forward(xmin)) - , umax_(Transform::forward(xmax)) - , du_((umax_ - umin_) / static_cast(N - 1)) - , inv_du_(1 / du_) - , N_(N) + : xmin_(xmin), xmax_(xmax), umin_(Transform::forward(xmin)), + umax_(Transform::forward(xmax)), + du_((umax_ - umin_) / static_cast(N - 1)), inv_du_(1 / du_), N_(N) { // A transform could be monotonically decreasing, so there's no guarantee // that umin_ < umax_ @@ -72,7 +68,8 @@ class RegularGrid1D { } // Returns closest index and weights for interpolation - PORTABLE_INLINE_FUNCTION void weights(const T &x, int &ix, weights_t &w) const { + PORTABLE_INLINE_FUNCTION void weights(const T &x, int &ix, + weights_t &w) const { const T u = Transform::forward(x); ix = index_u(u); const auto floor = static_cast(ix) * du_ + umin_; @@ -92,9 +89,7 @@ class RegularGrid1D { // (in)equality comparison PORTABLE_INLINE_FUNCTION bool operator==(const RegularGrid1D &other) const { - return (umin_ == other.umin_ && - umax_ == other.umax_ && - du_ == other.du_ && + return (umin_ == other.umin_ && umax_ == other.umax_ && du_ == other.du_ && N_ == other.N_); } PORTABLE_INLINE_FUNCTION bool @@ -107,12 +102,8 @@ class RegularGrid1D { PORTABLE_INLINE_FUNCTION T max() const { return xmax_; } PORTABLE_INLINE_FUNCTION size_t nPoints() const { return N_; } PORTABLE_INLINE_FUNCTION bool isnan() const { - return (std::isnan(xmin_) || - std::isnan(xmax_) || - std::isnan(umin_) || - std::isnan(umax_) || - std::isnan(du_) || - std::isnan(inv_du_) || + return (std::isnan(xmin_) || std::isnan(xmax_) || 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(); } @@ -178,7 +169,7 @@ class RegularGrid1D { } #endif -private: + private: // Forces x in the interval PORTABLE_INLINE_FUNCTION int bound(int ix) const { #ifndef SPINER_DISABLE_BOUNDS_CHECKS diff --git a/test/regular_grid_1d.cpp b/test/regular_grid_1d.cpp index 3839094c1..240519160 100644 --- a/test/regular_grid_1d.cpp +++ b/test/regular_grid_1d.cpp @@ -146,7 +146,7 @@ TEST_CASE("RegularGrid1D with test transformations", REQUIRE(gn.nPoints() == N); // TODO: Do we want the bounds to be exact? CHECK(gn.x(0) == min); - CHECK(gn.x(N-1) == max); + CHECK(gn.x(N - 1) == max); // Basic tests (expand) REQUIRE(ge.min() == min); @@ -154,7 +154,7 @@ TEST_CASE("RegularGrid1D with test transformations", REQUIRE(ge.nPoints() == N); // TODO: Do we want the bounds to be exact? CHECK(ge.x(0) == min); - CHECK(ge.x(N-1) == max); + CHECK(ge.x(N - 1) == max); // Basic tests (decrease) REQUIRE(gd.min() == min); @@ -162,7 +162,7 @@ TEST_CASE("RegularGrid1D with test transformations", REQUIRE(gd.nPoints() == N); // TODO: Do we want the bounds to be exact? CHECK(gd.x(0) == min); - CHECK(gd.x(N-1) == max); + CHECK(gd.x(N - 1) == max); // Check all fixed points (narrow) for (std::size_t n = 0; n < N; ++n) { From fce0f564335fde51350bb213a6ce0445d8359586 Mon Sep 17 00:00:00 2001 From: "Brendan K. Krueger" Date: Wed, 23 Oct 2024 10:55:29 -0600 Subject: [PATCH 48/48] format (again... doing it manually because clang-format here and on GitHub don't always agree) --- spiner/regular_grid_1d.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/spiner/regular_grid_1d.hpp b/spiner/regular_grid_1d.hpp index b7b371d70..fa3772d41 100644 --- a/spiner/regular_grid_1d.hpp +++ b/spiner/regular_grid_1d.hpp @@ -60,8 +60,7 @@ class RegularGrid1D { PORTABLE_INLINE_FUNCTION RegularGrid1D(T xmin, T xmax, size_t N) : xmin_(xmin), xmax_(xmax), umin_(Transform::forward(xmin)), umax_(Transform::forward(xmax)), - du_((umax_ - umin_) / static_cast(N - 1)), inv_du_(1 / du_), N_(N) - { + du_((umax_ - umin_) / static_cast(N - 1)), inv_du_(1 / du_), N_(N) { // A transform could be monotonically decreasing, so there's no guarantee // that umin_ < umax_ PORTABLE_ALWAYS_REQUIRE(xmin_ < xmax_ && N_ > 0, "Valid grid");