diff --git a/.appveyor.yml b/.appveyor.yml index c10a31646f..e0aaf575d1 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -121,7 +121,7 @@ before_build: # - cmd: if "%TARGET_ARCH%"=="x64" "C:\Program Files (x86)\Microsoft Visual Studio 15.9\VC\vcvarsall.bat" amd64 # CMake configure - - cmd: cmake -G "%OPENPMD_CMAKE_GENERATOR%" -DCMAKE_BUILD_TYPE=%CONFIGURATION% -DBUILD_SHARED_LIBS=%SHARED% -DBUILD_TESTING=ON -DCMAKE_INSTALL_PREFIX="%CONDA_INSTALL_LOCN%" -DCMAKE_INSTALL_BINDIR="Library\bin" ".." + - cmd: cmake -G "%OPENPMD_CMAKE_GENERATOR%" -DopenPMD_USE_REGIONS=OFF -DCMAKE_BUILD_TYPE=%CONFIGURATION% -DBUILD_SHARED_LIBS=%SHARED% -DBUILD_TESTING=ON -DCMAKE_INSTALL_PREFIX="%CONDA_INSTALL_LOCN%" -DCMAKE_INSTALL_BINDIR="Library\bin" ".." build_script: - cmd: cmake --build . --config %CONFIGURATION% diff --git a/.github/ci/spack-envs/gcc7_py36_ompi_h5_ad1_ad2/spack.yaml b/.github/ci/spack-envs/gcc7_py36_ompi_h5_ad1_ad2/spack.yaml new file mode 100644 index 0000000000..191de172ed --- /dev/null +++ b/.github/ci/spack-envs/gcc7_py36_ompi_h5_ad1_ad2/spack.yaml @@ -0,0 +1,61 @@ +# This is a Spack environment file. +# +# Activating and installing this environment will provide all dependencies +# that are needed for full-feature development. +# https//spack.readthedocs.io/en/latest/environments.html#anonymous-environments +# +spack: + specs: + - adios + - adios2 + - hdf5 + - openmpi + + packages: + adios: + variants: ~zfp ~sz ~lz4 ~blosc + adios2: + variants: ~zfp ~sz ~png ~dataman ~python ~fortran ~ssc ~shared ~bzip2 + cmake: + externals: + - spec: "cmake" + prefix: /usr + buildable: False + openmpi: + externals: + - spec: "openmpi" + prefix: /usr + buildable: False + perl: + externals: + - spec: "perl" + prefix: /usr + buildable: False + python: + externals: + - spec: "python" + prefix: /usr + buildable: False + all: + target: ['x86_64'] + variants: ~fortran + compiler: [gcc@7.4.0] + + compilers: + - compiler: + environment: {} + extra_rpaths: [] + flags: {} + modules: [] + operating_system: ubuntu18.04 + paths: + cc: /usr/bin/gcc-7 + cxx: /usr/bin/g++-7 + f77: /usr/bin/gfortran + fc: /usr/bin/gfortran + spec: gcc@7.4.0 + target: x86_64 + + config: + build_jobs: 2 + diff --git a/.github/workflows/unix.yml b/.github/workflows/unix.yml new file mode 100644 index 0000000000..b7c0ea895f --- /dev/null +++ b/.github/workflows/unix.yml @@ -0,0 +1,387 @@ +name: unix + +on: [push, pull_request] + +jobs: + clangtidy10_nopy_ompi_h5_ad1_ad2: + runs-on: ubuntu-20.04 + steps: + - uses: actions/checkout@v2 + - name: Spack Cache + uses: actions/cache@v2 + with: {path: /opt/spack, key: clangtidy10_nopy_ompi_h5_ad1_ad2 } + - name: Install + run: | + sudo apt-get update + sudo apt-get install clang clang-tidy gfortran libopenmpi-dev python + sudo .github/workflows/dependencies/install_spack + - name: Build + env: {CC: clang, CXX: clang++, CXXFLAGS: -Wno-deprecated-declarations} + # The main openPMD code base leads to many clang-tidy warnings + # when C++17 is enabled. We cannot enable it yet. + run: | + eval $(spack env activate --sh .github/ci/spack-envs/clangtidy_nopy_ompi_h5_ad1_ad2/) + spack install + + mkdir build && cd build + ../share/openPMD/download_samples.sh && chmod u-w samples/git-sample/*.h5 + cmake -S .. -B . -DCMAKE_BUILD_TYPE=RelWithDebInfo -DCMAKE_CXX_CLANG_TIDY="$(which clang-tidy);-system-headers=0" -DopenPMD_USE_REGIONS=ON -DopenPMD_USE_INVASIVE_TESTS=ON + cmake --build . --parallel 2 2> clang-tidy.log + cat clang-tidy.log + if [[ $(wc -m + +#include +#include +#include +#include +#include +#include + +//////////////////////////////////////////////////////////////////////////////// + +void point_example() { + using namespace openPMD::Regions; + using namespace std; + + cout << "\n" + << "Points (with dimension known at compile time):\n"; + + cout << " Define a point:\n"; + Point x{1, 2}; + cout << " x: " << x << "\n"; + + cout << " Define a point from a vector:\n"; + vector vec{4, 5}; + Point y(vec); + cout << " y: " << y << "\n"; + + cout << " Arithmetic operations:\n"; + auto z = x + 2 * y; + cout << " z: " << z << "\n"; + + cout << " Unit vectors pointing in direction d:\n"; + auto u0 = Point::unit(0); + auto u1 = Point::unit(1); + cout << " u0: " << u0 << "\n" + << " u1: " << u1 << "\n"; + cout << " A points with all elements the same:\n"; + auto p3 = Point::pure(3); + cout << " p3: " << p3 << "\n"; + + cout << " Element-wise operations:\n"; + // using std::abs, std::max; + auto mxy1 = max(abs(x), abs(y)); + // Apply arbitrary functions element-wise + auto mxy2 = fmap([](auto a, auto b) { return max(abs(a), abs(b)); }, x, y); + cout << " mxy1: " << mxy1 << "\n" + << " mxy2: " << mxy2 << "\n"; + + cout << " Reduction operations:\n"; + int mx1 = max_element(x); + // Apply arbitrary reduction operations + int mx2 = fold([](auto r, auto a) { return max(r, a); }, 0, x); + cout << " mx1: " << mx1 << "\n" + << " mx2: " << mx2 << "\n"; +} + +//////////////////////////////////////////////////////////////////////////////// + +void ndpoint_example() { + using namespace openPMD::Regions; + using namespace std; + + cout << "\n" + << "NDPoints (with dimension only known at run time):\n"; + + cout << " Define a point:\n"; + NDPoint x{1, 2}; + cout << " x: " << x << "\n"; + + cout << " Define a point from a vector:\n"; + vector vec{4, 5}; + NDPoint y(vec); + cout << " y: " << y << "\n"; + + // Arithmetic operations + cout << " Arithmetic operations:\n"; + auto z = x + 2 * y; + cout << " z: " << z << "\n"; + + cout << " Unit vectors pointing in direction d:\n"; + auto u0 = NDPoint::unit(2, 0); + auto u1 = NDPoint::unit(2, 1); + cout << " u0: " << u0 << "\n" + << " u1: " << u1 << "\n"; + cout << " A points with all elements the same:\n"; + auto p3 = NDPoint::pure(2, 3); + cout << " p3: " << p3 << "\n"; + + cout << " Element-wise operations:\n"; + // using std::abs, std::max; + auto mxy1 = max(abs(x), abs(y)); + // Apply arbitrary functions element-wise + auto mxy2 = fmap([](auto a, auto b) { return max(abs(a), abs(b)); }, x, y); + cout << " mxy1: " << mxy1 << "\n" + << " mxy2: " << mxy2 << "\n"; + + cout << " Reduction operations:\n"; + int mx1 = max_element(x); + // Apply arbitrary reduction operations + int mx2 = fold([](auto r, auto a) { return max(r, a); }, 0, x); + cout << " mx1: " << mx1 << "\n" + << " mx2: " << mx2 << "\n"; +} + +//////////////////////////////////////////////////////////////////////////////// + +void box_example() { + using namespace openPMD::Regions; + using namespace std; + + cout << "\n" + << "Boxes are spanned between points (inclusive lower, exclusive upper " + "bound):\n"; + + cout << " Define two points:\n"; + Point x{1, 4}, y{2, 5}; + cout << " x:" << x << "\n" + << " y:" << y << "\n"; + + cout << " Define a box between these points:\n"; + Box b(x, y); + cout << " b: " << b << " b.empty: " << b.empty() << "\n"; + cout << " Define an empty box:\n"; + Box be; + cout << " be: " << be << " be.empty: " << be.empty() << "\n"; + + cout << " Boxes can be shifted and scaled (e.g. to change resolution):\n"; + Point offset{1, 2}; + auto b1 = b >> offset; + Point scale{2, 2}; + auto b2 = b * scale; + // Boxes can be grown and shrunk (e.g. to add ghost zones) + auto bg = b.grown(1); + auto bs = b.shrunk(1); + cout << " shifted box: " << b1 << "\n" + << " scaled box: " << b2 << "\n" + << " grown box: " << bg << "\n" + << " shrunk box: " << bs << "\n"; + + cout << " The bounding box containing two boxes:\n"; + auto bb = bounding_box(b, b1); + cout << " bounding box: " << bb << "\n"; + + cout << " Boxes can be intersected:\n"; + auto bi = b & b1; + cout << " intersection: " << bi << "\n"; + + cout << " Set tests:\n"; + cout << " b == b1 (equality): " << (b == b1) << "\n" + << " b <= b1 (is_subset_of): " << (b <= b1) << "\n" + << " b < b1 (is_strict_subset_of): " << (b < b1) << "\n"; +} + +//////////////////////////////////////////////////////////////////////////////// + +void ndbox_example() { + using namespace openPMD::Regions; + using namespace std; + + cout << "\n" + << "Boxes are spanned between points (inclusive lower, exclusive upper " + "bound):\n"; + + cout << " Define two points:\n"; + NDPoint x{1, 4}, y{2, 5}; + cout << " x:" << x << "\n" + << " y:" << y << "\n"; + + cout << " Define a box between these points:\n"; + NDBox b(x, y); + cout << " b: " << b << " b.empty: " << b.empty() << "\n"; + cout << " Define an empty box:\n"; + NDBox be(2); + cout << " be: " << be << " be.empty: " << be.empty() << "\n"; + + cout << " Boxes can be shifted and scaled (e.g. to change resolution):\n"; + NDPoint offset{1, 2}; + auto b1 = b >> offset; + NDPoint scale{2, 2}; + auto b2 = b * scale; + // Boxes can be grown and shrunk (e.g. to add ghost zones) + auto bg = b.grown(1); + auto bs = b.shrunk(1); + cout << " shifted box: " << b1 << "\n" + << " scaled box: " << b2 << "\n" + << " grown box: " << bg << "\n" + << " shrunk box: " << bs << "\n"; + + cout << " The bounding box containing two boxes:\n"; + auto bb = bounding_box(b, b1); + cout << " bounding box: " << bb << "\n"; + + cout << " Boxes can be intersected:\n"; + auto bi = b & b1; + cout << " intersection: " << bi << "\n"; + + cout << " Set tests:\n"; + cout << " b == b1 (equality): " << (b == b1) << "\n" + << " b <= b1 (is_subset_of): " << (b <= b1) << "\n" + << " b < b1 (is_strict_subset_of): " << (b < b1) << "\n"; +} + +//////////////////////////////////////////////////////////////////////////////// + +void region_example() { + using namespace openPMD::Regions; + using namespace std; + + cout << "\n" + << "Regions consist of a set of boxes:\n"; + + cout << " Define two points:\n"; + Point x{1, 4}, y{2, 5}; + cout << " x:" << x << "\n" + << " y:" << y << "\n"; + + cout << " Define a box between these points:\n"; + Box b(x, y); + cout << " b:" << b << "\n"; + + cout << " Define a region consisting of this box:\n"; + Region r(b); + cout << " r: " << r << " r.empty: " << r.empty() << "\n"; + cout << " Define an empty region:\n"; + Region re; + cout << " re: " << re << " re.empty: " << re.empty() << "\n"; + + cout << " Regions can be shifted and scaled (e.g. to change resolution):\n"; + Point offset{1, 2}; + auto r1 = r >> offset; + Point scale{2, 2}; + auto r2 = r * scale; + // Regions can be grown and shrunk (e.g. to add ghost zones) + auto rg = r.grown(1); + auto rs = r.shrunk(1); + cout << " shifted region: " << r1 << "\n" + << " scaled region: " << r2 << "\n" + << " grown region: " << rg << "\n" + << " shrunk region: " << rs << "\n"; + + cout << " The bounding box containing a region:\n"; + auto bb = bounding_box(r); + cout << " bounding box: " << bb << "\n"; + + cout << " Regions can be treated as sets:\n"; + auto ri = r & r1; + auto ru = r | r1; + auto rd = r - r1; + auto rx = r ^ r1; + cout << " intersection: " << ri << "\n"; + cout << " union: " << ru << "\n"; + cout << " difference: " << rd << "\n"; + cout << " symmetric difference: " << rx << "\n"; + + cout << " Set tests:\n"; + cout << " r == r1 (equality): " << (r == r1) << "\n" + << " r <= r1 (is_subset_of): " << (r <= r1) << "\n" + << " r < r1 (is_strict_subset_of): " << (r < r1) << "\n"; + + cout << " Regions can be converted to a list of boxes:\n"; + cout << " rg - r:\n"; + for (const auto &bx : std::vector>(rg - r)) + cout << " " << bx << "\n"; + + // Don't benchmark the debug version; it uses an N^2 algorithm for self-checks + const int n = REGIONS_DEBUG ? 10 : 1000; + cout << " Create a list of " << n + << " 3d boxes and convert it to a region:\n"; + std::mt19937 gen; + std::uniform_int_distribution dist(0, 100); + const auto random_point = [&]() { + return Point::make([&](auto) { return dist(gen); }); + }; + const auto random_box = [&]() { + auto p1 = random_point(); + auto p2 = random_point(); + // Sort the points to avoid creating many empty boxes + return Box(min(p1, p2), max(p1, p2)); + }; + using clock = std::chrono::high_resolution_clock; + using timestamp = std::chrono::time_point; + using second = std::chrono::duration>; + const auto gettime = []() { return clock::now(); }; + const auto elapsed = [](timestamp t0, timestamp t1) { + return std::chrono::duration_cast(t1 - t0).count(); + }; + std::vector> boxlist; + for (int i = 0; i < n; ++i) + boxlist.push_back(random_box()); + timestamp t0 = gettime(); + Region reg(boxlist); + timestamp t1 = gettime(); + cout << " reg.size: " << reg.size() << "\n" + << " reg.nboxes: " << reg.nboxes() << "\n" + << " runtime: " << elapsed(t0, t1) << " sec\n"; + cout << " Grow the region by 1 point:\n"; + timestamp t2 = gettime(); + Region reg1 = reg.grown(1); + timestamp t3 = gettime(); + cout << " reg.size: " << reg1.size() << "\n" + << " reg.nboxes: " << reg1.nboxes() << "\n" + << " runtime: " << elapsed(t2, t3) << " sec\n"; +} + +//////////////////////////////////////////////////////////////////////////////// + +void ndregion_example() { + using namespace openPMD::Regions; + using namespace std; + + cout << "\n" + << "Regions consist of a set of boxes:\n"; + + cout << " Define two points:\n"; + NDPoint x{1, 4}, y{2, 5}; + cout << " x:" << x << "\n" + << " y:" << y << "\n"; + + cout << " Define a box between these points:\n"; + NDBox b(x, y); + cout << " b:" << b << "\n"; + + cout << " Define a region consisting of this box:\n"; + NDRegion r(b); + cout << " r: " << r << " r.empty: " << r.empty() << "\n"; + cout << " Define an empty region:\n"; + NDRegion re(2); + cout << " re: " << re << " re.empty: " << re.empty() << "\n"; + + cout << " Regions can be shifted and scaled (e.g. to change resolution):\n"; + NDPoint offset{1, 2}; + auto r1 = r >> offset; + NDPoint scale{2, 2}; + auto r2 = r * scale; + // Regions can be grown and shrunk (e.g. to add ghost zones) + auto rg = r.grown(1); + auto rs = r.shrunk(1); + cout << " shifted region: " << r1 << "\n" + << " scaled region: " << r2 << "\n" + << " grown region: " << rg << "\n" + << " shrunk region: " << rs << "\n"; + + cout << " The bounding box containing a region:\n"; + auto bb = bounding_box(r); + cout << " bounding box: " << bb << "\n"; + + cout << " Regions can be treated as sets:\n"; + auto ri = r & r1; + auto ru = r | r1; + auto rd = r - r1; + auto rx = r ^ r1; + cout << " intersection: " << ri << "\n"; + cout << " union: " << ru << "\n"; + cout << " difference: " << rd << "\n"; + cout << " symmetric difference: " << rx << "\n"; + + cout << " Set tests:\n"; + cout << " r == r1 (equality): " << (r == r1) << "\n" + << " r <= r1 (is_subset_of): " << (r <= r1) << "\n" + << " r < r1 (is_strict_subset_of): " << (r < r1) << "\n"; + + cout << " Regions can be converted to a list of boxes:\n"; + cout << " rg - r:\n"; + for (const auto &bx : std::vector>(rg - r)) + cout << " " << bx << "\n"; + + // Don't benchmark the debug version; it uses an N^2 algorithm for self-checks + const int n = REGIONS_DEBUG ? 10 : 1000; + cout << " Create a list of " << n + << " 3d boxes and convert it to a region:\n"; + std::mt19937 gen; + std::uniform_int_distribution dist(0, 100); + const auto random_point = [&]() { + return NDPoint::make(3, [&](auto) { return dist(gen); }); + }; + const auto random_box = [&]() { + auto p1 = random_point(); + auto p2 = random_point(); + // Sort the points to avoid creating many empty boxes + return NDBox(min(p1, p2), max(p1, p2)); + }; + using clock = std::chrono::high_resolution_clock; + using timestamp = std::chrono::time_point; + using second = std::chrono::duration>; + const auto gettime = []() { return clock::now(); }; + const auto elapsed = [](timestamp t0, timestamp t1) { + return std::chrono::duration_cast(t1 - t0).count(); + }; + std::vector> boxlist; + for (int i = 0; i < n; ++i) + boxlist.push_back(random_box()); + timestamp t0 = gettime(); + NDRegion reg(3, boxlist); + timestamp t1 = gettime(); + cout << " reg.size: " << reg.size() << "\n" + << " reg.nboxes: " << reg.nboxes() << "\n" + << " runtime: " << elapsed(t0, t1) << " sec\n"; + cout << " Grow the region by 1 point:\n"; + timestamp t2 = gettime(); + NDRegion reg1 = reg.grown(1); + timestamp t3 = gettime(); + cout << " reg.size: " << reg1.size() << "\n" + << " reg.nboxes: " << reg1.nboxes() << "\n" + << " runtime: " << elapsed(t2, t3) << " sec\n"; +} + +//////////////////////////////////////////////////////////////////////////////// + +int main() { + point_example(); + ndpoint_example(); + + box_example(); + ndbox_example(); + + region_example(); + ndregion_example(); + + return 0; +} diff --git a/include/openPMD/regions/Box.hpp b/include/openPMD/regions/Box.hpp new file mode 100644 index 0000000000..bcc7e762c6 --- /dev/null +++ b/include/openPMD/regions/Box.hpp @@ -0,0 +1,434 @@ +#pragma once + +#include "Helpers.hpp" +#include "Point.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace openPMD { +namespace Regions { + +/** A D-dimensional box + * + * A box is described by two points, its lower bound (inclusive) and + * upper bound (exclusive). If the lower bound not less than the upper + * bound, the box is empty. Empty boxes are fine (similar to an empty + * array). + * + * The dimension D needs to be known at compile time. @see NDBox + */ +template class Box; + +template class Box { + bool is_full; + + explicit constexpr Box(bool is_full_) : is_full(is_full_) {} + +public: + constexpr static std::size_t D = 0; + + using value_type = typename Point::value_type; + using size_type = typename Point::size_type; + + /** Create empty box + */ + constexpr Box() : is_full(false) {} + + Box(const Box &) = default; + Box(Box &&) = default; + Box &operator=(const Box &) = default; + Box &operator=(Box &&) = default; + + template Box(const Box &b) : is_full(b.is_full) {} + + /** Create box from lower (inclusive) and upper (exclusive) bound + */ + Box(const Point &, const Point &) : is_full(false) {} + /** Create box holding a single point + */ + explicit Box(const Point &) : is_full(true) {} + + // Predicates + size_type ndims() const { return D; } + bool empty() const { return !is_full; } + Point lower() const { return Point(); } + Point upper() const { return Point(); } + Point shape() const { return Point(); } + size_type size() const { return is_full; } + + // Shift and scale operators + Box &operator>>=(const Point &) { return *this; } + Box &operator<<=(const Point &) { return *this; } + Box &operator*=(const Point &) { return *this; } + Box operator>>(const Point &) const { return Box(*this); } + Box operator<<(const Point &) const { return Box(*this); } + Box operator*(const Point &) const { return Box(*this); } + Box grown(const Point &, const Point &) const { return *this; } + Box grown(const Point &) const { return *this; } + Box grown(const T &) const { return *this; } + Box shrunk(const Point &, const Point &) const { return *this; } + Box shrunk(const Point &) const { return *this; } + Box shrunk(const T &) const { return *this; } + + // Comparison operators + friend bool operator==(const Box &b1, const Box &b2) { + return b1.empty() == b2.empty(); + } + friend bool operator!=(const Box &b1, const Box &b2) { return !(b1 == b2); } + + // Set comparison operators + bool contains(const Point &) const { return !empty(); } + friend bool isdisjoint(const Box &b1, const Box &b2) { + return b1.empty() || b2.empty(); + } + friend bool operator<=(const Box &b1, const Box &b2) { + return !b1.empty() <= !b2.empty(); + } + friend bool operator>=(const Box &b1, const Box &b2) { return b2 <= b1; } + friend bool operator<(const Box &b1, const Box &b2) { + return b1 <= b2 && b1 != b2; + } + friend bool operator>(const Box &b1, const Box &b2) { return b2 < b1; } + bool is_subset_of(const Box &b) const { return *this <= b; } + bool is_superset_of(const Box &b) const { return *this >= b; } + bool is_strict_subset_of(const Box &b) const { return *this < b; } + bool is_strict_superset_of(const Box &b) const { return *this > b; } + + // Set operations + friend Box bounding_box(const Box &b1, const Box &b2) { + return Box(!b1.empty() || !b2.empty()); + } + + friend Box operator&(const Box &b1, const Box &b2) { + return Box(!b1.empty() & !b2.empty()); + } + Box &operator&=(const Box &b) { return *this = *this & b; } + friend Box intersection(const Box &b1, const Box &b2) { return b1 & b2; } + + friend bool operator==(const Box &b, const std::vector &bs) { + // This assumes that the elements of bs are disjoint + return b.empty() == bs.empty(); + } + friend bool operator==(const std::vector &bs, const Box &b) { + return b == bs; + } + friend bool operator!=(const Box &b, const std::vector &bs) { + return !(b == bs); + } + friend bool operator!=(const std::vector &bs, const Box &b) { + return !(bs == b); + } + + friend std::vector operator-(const Box &b1, const Box &b2) { + if (!b1.empty() > b2.empty()) + return std::vector{Box(Point())}; + return std::vector{}; + } + friend std::vector difference(const Box &b1, const Box &b2) { + return b1 - b2; + } + + friend std::vector operator|(const Box &b1, const Box &b2) { + if (!b1.empty() | !b2.empty()) + return std::vector{Box(Point())}; + return std::vector{}; + } + friend std::vector setunion(const Box &b1, const Box &b2) { + return b1 | b2; + } + + friend std::vector operator^(const Box &b1, const Box &b2) { + if (!b1.empty() ^ !b2.empty()) + return std::vector{Box(Point())}; + return std::vector{}; + } + friend std::vector symmetric_difference(const Box &b1, const Box &b2) { + return b1 ^ b2; + } + + /** Output a box + */ + friend std::ostream &operator<<(std::ostream &os, const Box &b) { + return os << "(" << b.is_full << ")"; + } +}; + +//////////////////////////////////////////////////////////////////////////////// + +template class Box { + Point lo, hi; + +public: + using value_type = typename Point::value_type; + using size_type = typename Point::size_type; + + /** Create empty box + */ + Box() = default; + + Box(const Box &) = default; + Box(Box &&) = default; + Box &operator=(const Box &) = default; + Box &operator=(Box &&) = default; + + template Box(const Box &b) : lo(b.lo), hi(b.hi) {} + + /** Create box from lower (inclusive) and upper (exclusive) bound + */ + Box(const Point &lo_, const Point &hi_) : lo(lo_), hi(hi_) {} + /** Create box holding a single point + */ + explicit Box(const Point &p) : lo(p), hi(p + T(1)) {} + + // Predicates + /** Number of dimensions + */ + size_type ndims() const { return D; } + /** Whether the Box is empty + */ + bool empty() const { return any(hi <= lo); } + /** Lower bound (inclusive) + */ + Point lower() const { return lo; } + /** Upper bound (exclusive) + */ + Point upper() const { return hi; } + /** Shape, i.e. the "size" in each direction + */ + Point shape() const { return max(hi - lo, T(0)); } + /** Size, the total number of contained points + */ + size_type size() const { return product(shape()); } + + // Shift and scale operators + /** Shift a Box to the right (upwards). The shift can be negative, which + * shifts left. + */ + Box &operator>>=(const Point &p) { + lo += p; + hi += p; + return *this; + } + /** Shift a Box to the left (downwards). The shift can be negative, which + * shifts right. + */ + Box &operator<<=(const Point &p) { + lo -= p; + hi -= p; + return *this; + } + /** Scale a Box + */ + Box &operator*=(const Point &p) { + lo *= p; + hi *= p; + return *this; + } + Box operator>>(const Point &p) const { return Box(*this) >>= p; } + Box operator<<(const Point &p) const { return Box(*this) <<= p; } + Box operator*(const Point &p) const { return Box(*this) *= p; } + /** Grow a Box by given amounts in each direction. + * + * The growth can be negative, which shrinks the box. If a Box is + * shrunk too much it becomes empty. Growing an empty box always + * results in an empty box. + */ + Box grown(const Point &dlo, const Point &dup) const { + if (empty()) + return Box(); + Box nb(*this); + nb.lo -= dlo; + nb.hi += dup; + return nb; + } + Box grown(const Point &d) const { return grown(d, d); } + Box grown(const T &d) const { return grown(Point::pure(d)); } + /** Shrink a Box by given amounts in each direction. + * + * The shrinkage can be negative, which grows the box. If a Box is + * shrunk too much it becomes empty. Growing an empty box always + * results in an empty box. + */ + Box shrunk(const Point &dlo, const Point &dup) const { + return grown(-dlo, -dup); + } + Box shrunk(const Point &d) const { return shrunk(d, d); } + Box shrunk(const T &d) const { return shrunk(Point::pure(d)); } + + // Comparison operators + /** Compare two boxes + * + * (All empty boxes are equal.) + */ + friend bool operator==(const Box &b1, const Box &b2) { + if (b1.empty() && b2.empty()) + return true; + if (b1.empty() || b2.empty()) + return false; + return all(b1.lo == b2.lo && b1.hi == b2.hi); + } + friend bool operator!=(const Box &b1, const Box &b2) { return !(b1 == b2); } + + // Set comparison operators + /** Check whether a Box contains a given point + */ + bool contains(const Point &p) const { + if (empty()) + return false; + return all(p >= lo && p < hi); + } + /** Check whether two boxes are disjoint, i.e. whether they have no point in + * common + */ + friend bool isdisjoint(const Box &b1, const Box &b2) { + return (b1 & b2).empty(); + } + /** Check whether Box 1 is (completely) contained in Box 2 + */ + friend bool operator<=(const Box &b1, const Box &b2) { + if (b1.empty()) + return true; + if (b2.empty()) + return false; + return all(b1.lo >= b2.lo && b1.hi <= b2.hi); + } + /** Check whether Box 1 (completely) contains Box 2 + */ + friend bool operator>=(const Box &b1, const Box &b2) { return b2 <= b1; } + /** Check whether Box 1 is (completely) contained in Box 2, and Box + * 2 is larger than Box 1 + */ + friend bool operator<(const Box &b1, const Box &b2) { + return b1 <= b2 && b1 != b2; + } + /** Check whether Box 1 (completely) contains in Box 2, and Box + * 1 is larger than Box 2 + */ + friend bool operator>(const Box &b1, const Box &b2) { return b2 < b1; } + /** Check whether a Box is a subset of another Box. This is equivalen to `<=`. + */ + bool is_subset_of(const Box &b) const { return *this <= b; } + /** Check whether a Box is a superset of another Box. This is equivalen to + * `>=`. + */ + bool is_superset_of(const Box &b) const { return *this >= b; } + /** Check whether a Box is a strict subset of another Box. This is equivalent + * to `<`. + */ + bool is_strict_subset_of(const Box &b) const { return *this < b; } + /** Check whether a Box is a strict superset of another Box. This is + * equivalent to `>`. + */ + bool is_strict_superset_of(const Box &b) const { return *this > b; } + + // Set operations + /** Calculate the bounding box of two Boxes. This is the smallest Box that + * contains both Boxes. + */ + friend Box bounding_box(const Box &b1, const Box &b2) { + if (b1.empty()) + return b2; + if (b2.empty()) + return b1; + auto r = Box(min(b1.lo, b2.lo), max(b1.hi, b2.hi)); +// Postcondition +#if REGIONS_DEBUG + assert(b1 <= r && b2 <= r); +#endif + return r; + } + + friend Box operator&(const Box &b1, const Box &b2) { + auto nlo = max(b1.lo, b2.lo); + auto nhi = min(b1.hi, b2.hi); + auto r = Box(nlo, nhi); +// Postcondition +#if REGIONS_DEBUG + assert(r <= b1 && r <= b2); +#endif + return r; + } + /** Calculate the intersection between two Boxes + * + * Other set operations (union, difference, symmetric difference) are not + * supported for Boxes; use Regions instead. + */ + Box &operator&=(const Box &b) { return *this = *this & b; } + /** Calculate the intersection between two Boxes + */ + friend Box intersection(const Box &b1, const Box &b2) { return b1 & b2; } + + /** Output a Box + */ + friend std::ostream &operator<<(std::ostream &os, const Box &b) { + return os << "(" << b.lo << ":" << b.hi << ")"; + } +}; + +} // namespace Regions +} // namespace openPMD + +namespace std { + +template +struct equal_to>; + +template struct equal_to> { + constexpr bool operator()(const openPMD::Regions::Box &x, + const openPMD::Regions::Box &y) const { + return x.empty() == y.empty(); + } +}; + +template +struct equal_to> { + constexpr bool operator()(const openPMD::Regions::Box &x, + const openPMD::Regions::Box &y) const { + if (x.empty() && y.empty()) + return true; + if (x.empty() || y.empty()) + return false; + return openPMD::Regions::helpers::tuple_eq( + make_tuple(x.lower(), x.upper()), make_tuple(y.lower(), y.upper())); + } +}; + +template struct hash> { + constexpr size_t operator()(const openPMD::Regions::Box &x) const { + if (x.empty()) + return size_t(0xc9df21a36550a048ULL); + hash> h; + return openPMD::Regions::helpers::hash_combine(h(x.lower()), h(x.upper())); + } +}; + +template struct less>; + +template struct less> { + constexpr bool operator()(const openPMD::Regions::Box &x, + const openPMD::Regions::Box &y) const { + return !x.empty() < !y.empty(); + } +}; + +template struct less> { + constexpr bool operator()(const openPMD::Regions::Box &x, + const openPMD::Regions::Box &y) const { + if (x.empty() && y.empty()) + return false; + if (x.empty()) + return true; + if (y.empty()) + return false; + return openPMD::Regions::helpers::tuple_lt( + make_tuple(x.lower(), x.upper()), make_tuple(y.lower(), y.upper())); + } +}; + +} // namespace std diff --git a/include/openPMD/regions/Helpers.hpp b/include/openPMD/regions/Helpers.hpp new file mode 100644 index 0000000000..84541e91db --- /dev/null +++ b/include/openPMD/regions/Helpers.hpp @@ -0,0 +1,219 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +namespace openPMD { +namespace Regions { + +#define REGIONS_DEBUG 1 // 0 or 1 + +namespace helpers { + +//////////////////////////////////////////////////////////////////////////////// + +// Combine hash values +template std::size_t hash_combine(std::size_t seed, const T &x) { + std::hash h; + // Taken from Boost + return seed ^ + (h(x) + size_t(0x00e6052366ac4440eULL) + (seed << 6) + (seed >> 2)); +} + +//////////////////////////////////////////////////////////////////////////////// + +// Convert a tuple-like type to a std::array +template +constexpr auto array_push(const Tuple &t, std::index_sequence, + const U &x) { + return std::array{{std::get(t)..., T(x)}}; +} + +// Append an element to a std::array +template +constexpr auto array_push(const std::array &a, const U &e) { + return array_push(a, std::make_index_sequence(), e); +} + +// Construct a std::array from a function +template +constexpr std::array construct_array(const F &f) { + if constexpr (N == 0) + return std::array(); + else + return array_push(construct_array(f), f(N - 1)); +} + +// Convert a std::vector +template +std::vector convert_vector(const std::vector &vec) { + std::vector res(vec.size()); + for (std::size_t i = 0; i < res.size(); ++i) + res[i] = vec[i]; + return res; +} +template +std::vector convert_vector(std::vector &&vec) { + std::vector res(vec.size()); + for (std::size_t i = 0; i < res.size(); ++i) + res[i] = std::move(vec[i]); + return res; +} + +//////////////////////////////////////////////////////////////////////////////// + +// Compare tuples + +namespace impl { + +template struct tuple_eq { + template + constexpr bool operator()(const Tuple1 &x, const Tuple2 &y) const { + using T1 = std::tuple_element_t; + using T2 = std::tuple_element_t; + using T = std::common_type_t; + const std::equal_to eq; + return tuple_eq()(x, y) && + eq(std::get(x), std::get(y)); + } +}; +template <> struct tuple_eq<0> { + template + constexpr bool operator()(const Tuple1 &, const Tuple2 &) const { + return true; + } +}; + +template struct tuple_cmp { + template + constexpr int operator()(const Tuple1 &x, const Tuple2 &y) const { + using T1 = std::tuple_element_t; + using T2 = std::tuple_element_t; + using T = std::common_type_t; + const int cmp = tuple_cmp()(x, y); + if (cmp != 0) + return cmp; + const std::less lt; + if (lt(std::get(x), std::get(y))) + return -1; + if (lt(std::get(y), std::get(x))) + return +1; + return 0; + } +}; +template <> struct tuple_cmp<0> { + template + constexpr int operator()(const Tuple1 &, const Tuple2 &) const { + return 0; + } +}; + +} // namespace impl + +template +constexpr bool tuple_eq(const Tuple1 &x, const Tuple2 &y) { + constexpr std::size_t sz = std::tuple_size_v; + static_assert(std::tuple_size_v == sz); + return impl::tuple_eq()(x, y); +} + +template +constexpr int tuple_cmp(const Tuple1 &x, const Tuple2 &y) { + constexpr std::size_t sz = std::tuple_size_v; + static_assert(std::tuple_size_v == sz); + return impl::tuple_cmp()(x, y); +} + +template +constexpr bool tuple_lt(const Tuple1 &x, const Tuple2 &y) { + return tuple_cmp(x, y) < 0; +} + +//////////////////////////////////////////////////////////////////////////////// + +// Compare vectors + +template +bool vector_eq(const Vector1 &x, const Vector2 &y) { + if (x.size() != y.size()) + return false; + using T1 = typename Vector1::value_type; + using T2 = typename Vector2::value_type; + using T = std::common_type_t; + const std::equal_to eq; + for (std::size_t n = 0; n < x.size(); ++n) + if (!eq(x[n], y[n])) + return false; + return true; +} + +template +bool vector_lt(const Vector1 &x, const Vector2 &y) { + using T1 = typename Vector1::value_type; + using T2 = typename Vector2::value_type; + using T = std::common_type_t; + const std::less lt; + using std::min; + for (std::size_t n = 0; n < min(x.size(), y.size()); ++n) + if (lt(x[n], y[n])) + return true; + else if (lt(y[n], x[n])) + return false; + return x.size() < y.size(); +} + +//////////////////////////////////////////////////////////////////////////////// + +// Reduce over a std::vector, using bisection to reduce cost + +/** Reduce a non-empty range + */ +template ::type> +T reduce1(const Op &op, std::vector &xs) { + assert(!xs.empty()); + for (std::size_t dist = 1; dist < xs.size(); dist *= 2) + for (std::size_t i = 0; i + dist < xs.size(); i += 2 * dist) + xs[i] = op(std::move(xs[i]), std::move(xs[i + dist])); + return std::move(xs[0]); +} + +/** Mapreduce a range with a given neutral element + */ +template +R mapreduce(const F &f, const Op &op, R z, const B &b, const E &e) { + std::vector rs; + for (auto i = b; i != e; ++i) + rs.push_back(f(*i)); + if (rs.empty()) + return z; + return reduce1(op, rs); +} + +/** Mapreduce a range + */ +template ::value_type, + typename R = std::decay_t>> +R mapreduce(const F &f, const Op &op, const I &b, const I &e) { + return mapreduce(f, op, R(), b, e); +} + +/** Mapreduce a container + */ +template >> +R mapreduce(const F &f, const Op &op, const C &c) { + return mapreduce(f, op, std::begin(c), std::end(c)); +} + +} // namespace helpers + +} // namespace Regions +} // namespace openPMD diff --git a/include/openPMD/regions/NDBox.hpp b/include/openPMD/regions/NDBox.hpp new file mode 100644 index 0000000000..d4e2883baf --- /dev/null +++ b/include/openPMD/regions/NDBox.hpp @@ -0,0 +1,527 @@ +#pragma once + +#include "Box.hpp" +#include "Helpers.hpp" +#include "NDPoint.hpp" + +namespace openPMD { +namespace Regions { + +template class NDBox; + +namespace detail { + +// Abstract base helper class + +template class VBox { +public: + using value_type = T; + using size_type = typename VPoint::size_type; + + virtual std::unique_ptr copy() const = 0; + + virtual ~VBox() {} + + virtual size_type ndims() const = 0; + virtual bool empty() const = 0; + virtual NDPoint lower() const = 0; + virtual NDPoint upper() const = 0; + virtual NDPoint shape() const = 0; + virtual size_type size() const = 0; + + virtual VBox &operator>>=(const NDPoint &p) = 0; + virtual VBox &operator<<=(const NDPoint &p) = 0; + virtual VBox &operator*=(const NDPoint &p) = 0; + virtual std::unique_ptr operator>>(const NDPoint &p) const = 0; + virtual std::unique_ptr operator<<(const NDPoint &p) const = 0; + virtual std::unique_ptr operator*(const NDPoint &p) const = 0; + virtual std::unique_ptr grown(const NDPoint &dlo, + const NDPoint &dup) const = 0; + virtual std::unique_ptr grown(const NDPoint &d) const = 0; + virtual std::unique_ptr grown(const T &d) const = 0; + virtual std::unique_ptr shrunk(const NDPoint &dlo, + const NDPoint &dup) const = 0; + virtual std::unique_ptr shrunk(const NDPoint &d) const = 0; + virtual std::unique_ptr shrunk(const T &d) const = 0; + + virtual bool operator==(const VBox &b2) const = 0; + virtual bool operator!=(const VBox &b2) const = 0; + + virtual bool contains(const NDPoint &p) const = 0; + virtual bool isdisjoint1(const VBox &b2) const = 0; + virtual bool operator<=(const VBox &b2) const = 0; + virtual bool operator>=(const VBox &b2) const = 0; + virtual bool operator<(const VBox &b2) const = 0; + virtual bool operator>(const VBox &b2) const = 0; + virtual bool is_subset_of(const VBox &b) const = 0; + virtual bool is_superset_of(const VBox &b) const = 0; + virtual bool is_strict_subset_of(const VBox &b) const = 0; + virtual bool is_strict_superset_of(const VBox &b) const = 0; + + virtual std::unique_ptr bounding_box1(const VBox &b2) const = 0; + + virtual std::unique_ptr operator&(const VBox &b2) const = 0; + virtual VBox &operator&=(const VBox &b) = 0; + virtual std::unique_ptr intersection1(const VBox &b2) const = 0; + + virtual bool equal_to1(const VBox &x) const = 0; + virtual std::size_t hash1() const = 0; + virtual bool less1(const VBox &x) const = 0; + + virtual std::ostream &output(std::ostream &os) const = 0; +}; + +// Helper class wrapping Box + +template class WBox final : public VBox { + Box b; + +public: + using typename VBox::value_type; + using typename VBox::size_type; + + WBox() : b{} {} + + WBox(const WBox &x) = default; + WBox(WBox &&) = default; + WBox &operator=(const WBox &) = default; + WBox &operator=(WBox &&) = default; + + WBox(const Box &b_) : b(b_) {} + WBox(Box &&b_) : b(std::move(b_)) {} + operator Box() const { return b; } + + WBox(const Point &lo_, const Point &hi_) : b(lo_, hi_) {} + explicit WBox(const Point &p) : b(p) {} + + std::unique_ptr> copy() const override { + return std::make_unique(*this); + } + + ~WBox() override {} + + size_type ndims() const override { return b.ndims(); } + bool empty() const override { return b.empty(); } + NDPoint lower() const override { return b.lower(); } + NDPoint upper() const override { return b.upper(); } + NDPoint shape() const override { return b.shape(); } + size_type size() const override { return b.size(); } + + VBox &operator>>=(const NDPoint &p) override { + b >>= Point(p); + return *this; + } + VBox &operator<<=(const NDPoint &p) override { + b <<= Point(p); + return *this; + } + VBox &operator*=(const NDPoint &p) override { + b *= Point(p); + return *this; + } + std::unique_ptr> operator>>(const NDPoint &p) const override { + return std::make_unique(b >> Point(p)); + } + std::unique_ptr> operator<<(const NDPoint &p) const override { + return std::make_unique(b << Point(p)); + } + std::unique_ptr> operator*(const NDPoint &p) const override { + return std::make_unique(b * Point(p)); + } + std::unique_ptr> grown(const NDPoint &dlo, + const NDPoint &dup) const override { + return std::make_unique(b.grown(Point(dlo), Point(dup))); + } + std::unique_ptr> grown(const NDPoint &d) const override { + return std::make_unique(b.grown(Point(d))); + } + std::unique_ptr> grown(const T &d) const override { + return std::make_unique(b.grown(d)); + } + std::unique_ptr> shrunk(const NDPoint &dlo, + const NDPoint &dup) const override { + return std::make_unique(b.shrunk(Point(dlo), Point(dup))); + } + std::unique_ptr> shrunk(const NDPoint &d) const override { + return std::make_unique(b.shrunk(Point(d))); + } + std::unique_ptr> shrunk(const T &d) const override { + return std::make_unique(b.shrunk(d)); + } + + bool operator==(const VBox &b2) const override { + return b == dynamic_cast(b2).b; + } + bool operator!=(const VBox &b2) const override { + return b != dynamic_cast(b2).b; + } + + bool contains(const NDPoint &p) const override { + return b.contains(Point(p)); + } + bool isdisjoint1(const VBox &b2) const override { + return isdisjoint(b, dynamic_cast(b2).b); + } + bool operator<=(const VBox &b2) const override { + return b <= dynamic_cast(b2).b; + } + bool operator>=(const VBox &b2) const override { + return b >= dynamic_cast(b2).b; + } + bool operator<(const VBox &b2) const override { + return b < dynamic_cast(b2).b; + } + bool operator>(const VBox &b2) const override { + return b > dynamic_cast(b2).b; + } + bool is_subset_of(const VBox &b2) const override { + return b.is_subset_of(dynamic_cast(b2).b); + } + bool is_superset_of(const VBox &b2) const override { + return b.is_superset_of(dynamic_cast(b2).b); + } + bool is_strict_subset_of(const VBox &b2) const override { + return b.is_strict_subset_of(dynamic_cast(b2).b); + } + bool is_strict_superset_of(const VBox &b2) const override { + return b.is_strict_superset_of(dynamic_cast(b2).b); + } + + std::unique_ptr> bounding_box1(const VBox &b2) const override { + return std::make_unique( + bounding_box(b, dynamic_cast(b2).b)); + } + + std::unique_ptr> operator&(const VBox &b2) const override { + return std::make_unique(b & dynamic_cast(b2).b); + } + VBox &operator&=(const VBox &b2) override { + b &= dynamic_cast(b2).b; + return *this; + } + std::unique_ptr> intersection1(const VBox &b2) const override { + return std::make_unique( + intersection(b, dynamic_cast(b2).b)); + } + + bool equal_to1(const VBox &x) const override { + return std::equal_to>()(b, dynamic_cast(x).b); + } + std::size_t hash1() const override { return std::hash>()(b); } + bool less1(const VBox &x) const override { + return std::less>()(b, dynamic_cast(x).b); + } + + std::ostream &output(std::ostream &os) const override { return os << b; } +}; + +template +std::unique_ptr> make_VBox(const std::size_t D, const Args &...args) { + switch (D) { + case 0: + return std::make_unique>(args...); + case 1: + return std::make_unique>(args...); + case 2: + return std::make_unique>(args...); + case 3: + return std::make_unique>(args...); + case 4: + return std::make_unique>(args...); + case 5: + return std::make_unique>(args...); + default: + abort(); + } +} + +} // namespace detail + +/** A NDBox + * + * The dimension (number of component) of the Box is only known at + * run-time. @see Box + * + * a NDBox is described by two points, its lower bound (inclusive) and + * upper bound (exclusive). If the lower bound not less than the upper + * bound, the box is empty. Empty boxes are fine (similar to an empty + * array). + */ +template class NDBox { + template using VBox = detail::VBox; + + template friend class NDBox; + friend struct std::equal_to; + friend struct std::hash; + friend struct std::less; + + std::unique_ptr> b; + + NDBox(std::unique_ptr> b_) : b(std::move(b_)) {} + +public: + /** Component type + */ + using value_type = typename VBox::value_type; + /** Return type of Box::size() + */ + using size_type = typename VBox::size_type; + + /** Create an invalid Box + */ + NDBox() : b() {} + /** Create a value-initialized Box with D components + */ + NDBox(const size_type D) : b(detail::make_VBox(D)) {} + + NDBox(const NDBox &x) : b(x.b ? x.b->copy() : nullptr) {} + NDBox(NDBox &&) = default; + NDBox &operator=(const NDBox &x) { + b = x.b ? x.b->copy() : nullptr; + return *this; + } + NDBox &operator=(NDBox &&) = default; + + /** Create an NDBox from a Box + */ + template + NDBox(const Box &b_) : b(std::make_unique>(b_)) {} + template + NDBox(Box &&b_) + : b(std::make_unique>(std::move(b_))) {} + /** Convert an NDBox to a Box + * + * This only works when D == ndims(). + */ + template operator Box() const { + return Box(dynamic_cast &>(*b)); + } + + /** Create an NDBox from lower (inclusive) and upper (exclusive) bound + */ + NDBox(const NDPoint &lo_, const NDPoint &hi_) + : b(detail::make_VBox(lo_.ndims(), lo_, hi_)) {} + /** Create box holding a single point + */ + explicit NDBox(const NDPoint &p) : b(detail::make_VBox(p.ndims(), p)) {} + + /** Check whether an NDBox is valid + * + * A valid NDBox knows its number of dimensions, and its components + * are initialized. An invalid NDBox does not know its number of + * dimensions and holds no data, similar to a null pointer. + * + * Most other member functions must not be called for invalid NDBox. + * + * Note that there is a difference between invalid NDBoxes and empty + * NDBoxes. Invalid NDBoxes are essentially invalid objects that + * cannot be used in a useful manner. Empty NDBoxes are fine, + * similar to empty arrays. + */ + bool has_value() const { return bool(b); } + + /** Number of dimensions + */ + size_type ndims() const { return b->ndims(); } + /** Whether the NDBox box is empty + */ + bool empty() const { return b->empty(); } + /** Lower bound (inclusive) + */ + NDPoint lower() const { return b->lower(); } + /** Upper bound (exclusive) + */ + NDPoint upper() const { return b->upper(); } + /** Shape, i.e. the "size" in each direction + */ + NDPoint shape() const { return b->shape(); } + /** Size, the total number of contained points + */ + size_type size() const { return b->size(); } + + /** Shift a NDBox to the right (upwards). The shift can be negative, which + * shifts left. + */ + NDBox &operator>>=(const NDPoint &p) { + *b >>= p->p; + return *this; + } + /** Shift a NDBox to the left (downwards). The shift can be negative, which + * shifts right. + */ + NDBox &operator<<=(const NDPoint &p) { + *b <<= p->p; + return *this; + } + /** Scale a NDBox + */ + NDBox &operator*=(const NDPoint &p) { + *b *= p->p; + return *this; + } + NDBox operator>>(const NDPoint &p) const { return NDBox(*b >> p); } + NDBox operator<<(const NDPoint &p) const { return NDBox(*b << p); } + NDBox operator*(const NDPoint &p) const { return NDBox(*b * p); } + /** Grow a NDBox by given amounts in each direction. + * + * The growth can be negative, which shrinks the box. If a NDBox is + * shrunk too much it becomes empty. Growing an empty box always + * results in an empty box. + */ + NDBox grown(const NDPoint &dlo, const NDPoint &dup) const { + return NDBox(b->grown(dlo, dup)); + } + NDBox grown(const NDPoint &d) const { return NDBox(b->grown(d)); } + NDBox grown(const T &d) const { return NDBox(b->grown(d)); } + NDBox shrunk(const NDPoint &dlo, const NDPoint &dup) const { + return NDBox(b->shrunk(dlo, dup)); + } + /** Shrink a NDBox by given amounts in each direction. + * + * The shrinkage can be negative, which grows the box. If a NDBox is + * shrunk too much it becomes empty. Growing an empty box always + * results in an empty box. + */ + NDBox shrunk(const NDPoint &d) const { return NDBox(b->shrunk(d)); } + NDBox shrunk(const T &d) const { return NDBox(b->shrunk(d)); } + + /** Compare two boxes + * + * (All empty boxes are equal.) + */ + friend bool operator==(const NDBox &b1, const NDBox &b2) { + return *b1.b == *b2.b; + } + friend bool operator!=(const NDBox &b1, const NDBox &b2) { + return *b1.b != *b2.b; + } + + /** Check whether a NDBox contains a given point + */ + bool contains(const NDPoint &p) const { return b->contains(p); } + /** Check whether two boxes are disjoint, i.e. whether they have no point in + * common + */ + friend bool isdisjoint(const NDBox &b1, const NDBox &b2) { + return b1.b->isdisjoint1(*b2.b); + } + /** Check whether Box 1 is (completely) contained in Box 2 + */ + friend bool operator<=(const NDBox &b1, const NDBox &b2) { + return *b1.b <= *b2.b; + } + /** Check whether Box 1 (completely) contains Box 2 + */ + friend bool operator>=(const NDBox &b1, const NDBox &b2) { + return *b1.b >= *b2.b; + } + /** Check whether Box 1 is (completely) contained in Box 2, and Box + * 2 is larger than Box 1 + */ + friend bool operator<(const NDBox &b1, const NDBox &b2) { + return *b1.b < *b2.b; + } + /** Check whether Box 1 (completely) contains in Box 2, and Box + * 1 is larger than Box 2 + */ + friend bool operator>(const NDBox &b1, const NDBox &b2) { + return *b1.b > *b2.b; + } + /** Check whether a NDBox is a subset of another Box. This is equivalent to + * `<=`. + */ + bool is_subset_of(const NDBox &b2) const { return b->is_subset_of(*b2.b); } + /** Check whether a NDBox is a superset of another Box. This is equivalent to + * `>=`. + */ + bool is_superset_of(const NDBox &b2) const { + return b->is_superset_of(*b2.b); + } + /** Check whether a NDBox is a strict subset of another Box. This is + * equivalent to `<`. + */ + bool is_strict_subset_of(const NDBox &b2) const { + return b->is_strict_subset_of(*b2.b); + } + /** Check whether a NDBox is a strict superset of another Box. This is + * equivalent to `>`. + */ + bool is_strict_superset_of(const NDBox &b2) const { + return b->is_strict_superset_of(*b2.b); + } + + /** Calculate the bounding box of two Boxes. This is the smallest Box that + * contains both Boxes. + */ + friend NDBox bounding_box(const NDBox &b1, const NDBox &b2) { + return b1.b->bounding_box1(*b2.b); + } + + /** Calculate the intersection between two Boxes + * + * Other set operations (union, difference, symmetric difference) are not + * supported for Boxes; use Regions instead. + */ + friend NDBox operator&(const NDBox &b1, const NDBox &b2) { + return *b1.b & *b2.b; + } + NDBox &operator&=(const NDBox &b2) { + *b &= *b2.b; + return *this; + } + /** Calculate the intersection between two Boxes + * + * Other set operations (union, difference, symmetric difference) are not + * supported for Boxes; use Regions instead. + */ + friend NDBox intersection(const NDBox &b1, const NDBox &b2) { + return intersection(*b1.b, *b2.b); + } + + /** Output an NDBox + */ + friend std::ostream &operator<<(std::ostream &os, const NDBox &x) { + if (x.b) + x.b->output(os); + else + os << "(INVALID)"; + return os; + } +}; + +} // namespace Regions +} // namespace openPMD + +namespace std { + +template struct equal_to> { + constexpr bool operator()(const openPMD::Regions::NDBox &x, + const openPMD::Regions::NDBox &y) const { + if (!x.has_value() && !y.has_value()) + return true; + if (!x.has_value() || !y.has_value()) + return false; + return x.b->equal_to1(*y.b); + } +}; + +template struct hash> { + constexpr size_t operator()(const openPMD::Regions::NDBox &x) const { + if (!x.has_value()) + return size_t(0x7e21b87749864bbcULL); + return x.b->hash1(); + } +}; + +template struct less> { + constexpr bool operator()(const openPMD::Regions::NDBox &x, + const openPMD::Regions::NDBox &y) const { + if (!x.has_value()) + return true; + if (!y.has_value()) + return false; + return x.b->less1(*y.b); + } +}; + +} // namespace std diff --git a/include/openPMD/regions/NDPoint.hpp b/include/openPMD/regions/NDPoint.hpp new file mode 100644 index 0000000000..b85bb85009 --- /dev/null +++ b/include/openPMD/regions/NDPoint.hpp @@ -0,0 +1,1238 @@ +#pragma once + +#include "Helpers.hpp" +#include "Point.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace openPMD { +namespace Regions { + +/** Maximum number of dimensions supported by NDPoint + */ +constexpr std::size_t max_ndims = 5; + +namespace detail { + +// Abstract base helper class + +template class VPoint { +public: + using value_type = T; + using size_type = std::size_t; + + virtual std::unique_ptr copy() const = 0; + + virtual ~VPoint() {} + + virtual std::unique_ptr + fmap1(const std::function &f) const = 0; + virtual std::unique_ptr + fmap1(const std::function &f, + const VPoint &x) const = 0; + virtual std::unique_ptr + fmap1(const std::function &f, + const VPoint &x, const VPoint &y) const = 0; + virtual T fold1(const std::function &op, + const T &r) const = 0; + virtual T fold1(const std::function &op, + const T &r, const VPoint &x) const = 0; + + virtual void set_from_vector(const std::vector &vec) = 0; + virtual operator std::vector() const = 0; + + virtual size_type size() const = 0; + + virtual size_type ndims() const = 0; + virtual const T &operator[](const size_type d) const = 0; + virtual T &operator[](const size_type d) = 0; + virtual std::unique_ptr erase(const size_type d) const = 0; + virtual std::unique_ptr insert(const size_type d, + const T &a) const = 0; + virtual std::unique_ptr reversed() const = 0; + + virtual std::unique_ptr operator+() const = 0; + virtual std::unique_ptr operator-() const = 0; + virtual std::unique_ptr operator~() const = 0; + virtual std::unique_ptr> operator!() const = 0; + + virtual std::unique_ptr operator+(const VPoint &x) const = 0; + virtual std::unique_ptr operator-(const VPoint &x) const = 0; + virtual std::unique_ptr operator*(const VPoint &x) const = 0; + virtual std::unique_ptr operator/(const VPoint &x) const = 0; + virtual std::unique_ptr operator%(const VPoint &x) const = 0; + virtual std::unique_ptr operator&(const VPoint &x) const = 0; + virtual std::unique_ptr operator|(const VPoint &x) const = 0; + virtual std::unique_ptr operator^(const VPoint &x) const = 0; + virtual std::unique_ptr> operator&&(const VPoint &x) const = 0; + virtual std::unique_ptr> operator||(const VPoint &x) const = 0; + + virtual std::unique_ptr left_plus(const T &a) const = 0; + virtual std::unique_ptr left_minus(const T &a) const = 0; + virtual std::unique_ptr left_multiplies(const T &a) const = 0; + virtual std::unique_ptr left_divides(const T &a) const = 0; + virtual std::unique_ptr left_modulus(const T &a) const = 0; + virtual std::unique_ptr left_bit_and(const T &a) const = 0; + virtual std::unique_ptr left_bit_or(const T &a) const = 0; + virtual std::unique_ptr left_bit_xor(const T &a) const = 0; + virtual std::unique_ptr> left_logical_and(const T &a) const = 0; + virtual std::unique_ptr> left_logical_or(const T &a) const = 0; + + virtual std::unique_ptr operator+(const T &b) const = 0; + virtual std::unique_ptr operator-(const T &b) const = 0; + virtual std::unique_ptr operator*(const T &b) const = 0; + virtual std::unique_ptr operator/(const T &b) const = 0; + virtual std::unique_ptr operator%(const T &b) const = 0; + virtual std::unique_ptr operator&(const T &b) const = 0; + virtual std::unique_ptr operator|(const T &b) const = 0; + virtual std::unique_ptr operator^(const T &b) const = 0; + virtual std::unique_ptr> operator&&(const T &b) const = 0; + virtual std::unique_ptr> operator||(const T &b) const = 0; + + virtual VPoint &operator+=(const VPoint &x) = 0; + virtual VPoint &operator-=(const VPoint &x) = 0; + virtual VPoint &operator*=(const VPoint &x) = 0; + virtual VPoint &operator/=(const VPoint &x) = 0; + virtual VPoint &operator%=(const VPoint &x) = 0; + virtual VPoint &operator&=(const VPoint &x) = 0; + virtual VPoint &operator|=(const VPoint &x) = 0; + virtual VPoint &operator^=(const VPoint &x) = 0; + + virtual VPoint &operator+=(const T &b) = 0; + virtual VPoint &operator-=(const T &b) = 0; + virtual VPoint &operator*=(const T &b) = 0; + virtual VPoint &operator/=(const T &b) = 0; + virtual VPoint &operator%=(const T &b) = 0; + virtual VPoint &operator&=(const T &b) = 0; + virtual VPoint &operator|=(const T &b) = 0; + virtual VPoint &operator^=(const T &b) = 0; + + virtual std::unique_ptr abs1() const = 0; + virtual std::unique_ptr fabs1() const = 0; + + virtual std::unique_ptr fmax1(const VPoint &x) const = 0; + virtual std::unique_ptr fmin1(const VPoint &x) const = 0; + virtual std::unique_ptr max1(const VPoint &x) const = 0; + virtual std::unique_ptr min1(const VPoint &x) const = 0; + + virtual std::unique_ptr fmax1(const T &b) const = 0; + virtual std::unique_ptr fmin1(const T &b) const = 0; + virtual std::unique_ptr max1(const T &b) const = 0; + virtual std::unique_ptr min1(const T &b) const = 0; + + virtual T all1() const = 0; + virtual T any1() const = 0; + virtual T max_element1() const = 0; + virtual T min_element1() const = 0; + virtual T product1() const = 0; + virtual T sum1() const = 0; + + virtual std::unique_ptr> operator==(const VPoint &x) const = 0; + virtual std::unique_ptr> operator!=(const VPoint &x) const = 0; + virtual std::unique_ptr> operator<(const VPoint &x) const = 0; + virtual std::unique_ptr> operator>(const VPoint &x) const = 0; + virtual std::unique_ptr> operator<=(const VPoint &x) const = 0; + virtual std::unique_ptr> operator>=(const VPoint &x) const = 0; + + virtual std::unique_ptr> operator==(const T &a) const = 0; + virtual std::unique_ptr> operator!=(const T &a) const = 0; + virtual std::unique_ptr> operator<(const T &a) const = 0; + virtual std::unique_ptr> operator>(const T &a) const = 0; + virtual std::unique_ptr> operator<=(const T &a) const = 0; + virtual std::unique_ptr> operator>=(const T &a) const = 0; + + virtual bool equal_to1(const VPoint &x) const = 0; + virtual std::size_t hash1() const = 0; + virtual bool less1(const VPoint &x) const = 0; + + virtual std::ostream &output(std::ostream &os) const = 0; +}; + +// Helper class wrapping Point + +template class WPoint final : public VPoint { + Point p; + +public: + using typename VPoint::value_type; + using typename VPoint::size_type; + + WPoint() : p{} {} + + WPoint(const WPoint &x) = default; + WPoint(WPoint &&) = default; + WPoint &operator=(const WPoint &) = default; + WPoint &operator=(WPoint &&) = default; + + WPoint(const Point &p_) : p(p_) {} + WPoint(Point &&p_) : p(std::move(p_)) {} + operator Point() const { return p; } + + WPoint(const std::array &arr) : p(arr) {} + WPoint(std::array &&arr) : p(std::move(arr)) {} + operator std::array() const { return p; } + + std::unique_ptr> copy() const override { + return std::make_unique(*this); + } + + ~WPoint() override {} + + std::unique_ptr> + fmap1(const std::function &f) const override { + return std::make_unique(fmap(f, p)); + } + std::unique_ptr> + fmap1(const std::function &f, + const VPoint &x) const override { + return std::make_unique( + fmap(f, p, dynamic_cast(x).p)); + } + std::unique_ptr> + fmap1(const std::function &f, + const VPoint &x, const VPoint &y) const override { + return std::make_unique(fmap(f, p, + dynamic_cast(x).p, + dynamic_cast(y).p)); + } + T fold1(const std::function &op, + const T &r) const override { + return fold(op, r, p); + } + T fold1(const std::function &op, + const T &r, const VPoint &x) const override { + return fold(op, r, p, dynamic_cast(x).p); + } + + void set_from_vector(const std::vector &vec) override { + *this = WPoint(vec); + } + operator std::vector() const override { return std::vector(p); } + + /*constexpr*/ size_type size() const override { return p.size(); } + + /*constexpr*/ size_type ndims() const override { return p.ndims(); } + + const T &operator[](const size_type d) const override { return p[d]; } + T &operator[](const size_type d) override { return p[d]; } + + std::unique_ptr> erase(const size_type d + [[maybe_unused]]) const override { + if constexpr (D == 0) + return std::unique_ptr>(); + else + return std::make_unique>(p.erase(d)); + } + std::unique_ptr> insert(const size_type d [[maybe_unused]], + const T &a) const override { + if constexpr (D == max_ndims) + return std::unique_ptr>(); + else + return std::make_unique>(p.insert(d, a)); + } + + std::unique_ptr> reversed() const override { + return std::make_unique(p.reversed()); + } + + std::unique_ptr> operator+() const override { + return std::make_unique(+p); + } + std::unique_ptr> operator-() const override { + return std::make_unique(-p); + } + std::unique_ptr> operator~() const override { + if constexpr (std::is_integral_v) + return std::make_unique(~p); + else + std::abort(); + } + std::unique_ptr> operator!() const override { + return std::make_unique>(!p); + } + + std::unique_ptr> operator+(const VPoint &x) const override { + return std::make_unique(p + dynamic_cast(x).p); + } + std::unique_ptr> operator-(const VPoint &x) const override { + return std::make_unique(p - dynamic_cast(x).p); + } + std::unique_ptr> operator*(const VPoint &x) const override { + return std::make_unique(p * dynamic_cast(x).p); + } + std::unique_ptr> operator/(const VPoint &x) const override { + return std::make_unique(p / dynamic_cast(x).p); + } + std::unique_ptr> operator%(const VPoint &x) const override { + if constexpr (std::is_integral_v) + return std::make_unique(p % dynamic_cast(x).p); + else + std::abort(); + } + std::unique_ptr> operator&(const VPoint &x) const override { + if constexpr (std::is_integral_v) + return std::make_unique(p & dynamic_cast(x).p); + else + std::abort(); + } + std::unique_ptr> operator|(const VPoint &x) const override { + if constexpr (std::is_integral_v) + return std::make_unique(p | dynamic_cast(x).p); + else + std::abort(); + } + std::unique_ptr> operator^(const VPoint &x) const override { + if constexpr (std::is_integral_v) + return std::make_unique(p ^ dynamic_cast(x).p); + else + std::abort(); + } + std::unique_ptr> operator&&(const VPoint &x) const override { + return std::make_unique>(p && + dynamic_cast(x).p); + } + std::unique_ptr> operator||(const VPoint &x) const override { + return std::make_unique>(p || + dynamic_cast(x).p); + } + + std::unique_ptr> left_plus(const T &a) const override { + return std::make_unique(a + p); + } + std::unique_ptr> left_minus(const T &a) const override { + return std::make_unique(a - p); + } + std::unique_ptr> left_multiplies(const T &a) const override { + return std::make_unique(a * p); + } + std::unique_ptr> left_divides(const T &a) const override { + return std::make_unique(a / p); + } + std::unique_ptr> left_modulus(const T &a) const override { + if constexpr (std::is_integral_v) + return std::make_unique(a % p); + else + std::abort(); + } + std::unique_ptr> left_bit_and(const T &a) const override { + if constexpr (std::is_integral_v) + return std::make_unique(a & p); + else + std::abort(); + } + std::unique_ptr> left_bit_or(const T &a) const override { + if constexpr (std::is_integral_v) + return std::make_unique(a | p); + else + std::abort(); + } + std::unique_ptr> left_bit_xor(const T &a) const override { + if constexpr (std::is_integral_v) + return std::make_unique(a ^ p); + else + std::abort(); + } + std::unique_ptr> left_logical_and(const T &a) const override { + return std::make_unique>(a && p); + } + std::unique_ptr> left_logical_or(const T &a) const override { + return std::make_unique>(a || p); + } + + std::unique_ptr> operator+(const T &b) const override { + return std::make_unique(p + b); + } + std::unique_ptr> operator-(const T &b) const override { + return std::make_unique(p - b); + } + std::unique_ptr> operator*(const T &b) const override { + return std::make_unique(p * b); + } + std::unique_ptr> operator/(const T &b) const override { + return std::make_unique(p / b); + } + std::unique_ptr> operator%(const T &b) const override { + if constexpr (std::is_integral_v) + return std::make_unique(p % b); + else + std::abort(); + } + std::unique_ptr> operator&(const T &b) const override { + if constexpr (std::is_integral_v) + return std::make_unique(p & b); + else + std::abort(); + } + std::unique_ptr> operator|(const T &b) const override { + if constexpr (std::is_integral_v) + return std::make_unique(p | b); + else + std::abort(); + } + std::unique_ptr> operator^(const T &b) const override { + if constexpr (std::is_integral_v) + return std::make_unique(p ^ b); + else + std::abort(); + } + std::unique_ptr> operator&&(const T &b) const override { + return std::make_unique>(p && b); + } + std::unique_ptr> operator||(const T &b) const override { + return std::make_unique>(p || b); + } + + VPoint &operator+=(const VPoint &x) override { + p += dynamic_cast(x).p; + return *this; + } + VPoint &operator-=(const VPoint &x) override { + p -= dynamic_cast(x).p; + return *this; + } + VPoint &operator*=(const VPoint &x) override { + p *= dynamic_cast(x).p; + return *this; + } + VPoint &operator/=(const VPoint &x) override { + p /= dynamic_cast(x).p; + return *this; + } + VPoint &operator%=(const VPoint &x) override { + if constexpr (std::is_integral_v) { + p %= dynamic_cast(x).p; + return *this; + } else { + std::abort(); + } + } + VPoint &operator&=(const VPoint &x) override { + if constexpr (std::is_integral_v) { + p &= dynamic_cast(x).p; + return *this; + } else { + std::abort(); + } + } + VPoint &operator|=(const VPoint &x) override { + if constexpr (std::is_integral_v) { + p |= dynamic_cast(x).p; + return *this; + } else { + std::abort(); + } + } + VPoint &operator^=(const VPoint &x) override { + if constexpr (std::is_integral_v) { + p ^= dynamic_cast(x).p; + return *this; + } else { + std::abort(); + } + } + + VPoint &operator+=(const T &b) override { + p += b; + return *this; + } + VPoint &operator-=(const T &b) override { + p -= b; + return *this; + } + VPoint &operator*=(const T &b) override { + p *= b; + return *this; + } + VPoint &operator/=(const T &b) override { + p /= b; + return *this; + } + VPoint &operator%=(const T &b) override { + if constexpr (std::is_integral_v) { + p %= b; + return *this; + } else { + std::abort(); + } + } + VPoint &operator&=(const T &b) override { + if constexpr (std::is_integral_v) { + p &= b; + return *this; + } else { + std::abort(); + } + } + VPoint &operator|=(const T &b) override { + if constexpr (std::is_integral_v) { + p |= b; + return *this; + } else { + std::abort(); + } + } + VPoint &operator^=(const T &b) override { + if constexpr (std::is_integral_v) { + p ^= b; + return *this; + } else { + std::abort(); + } + } + + std::unique_ptr> abs1() const override { + return std::make_unique(abs(p)); + } + std::unique_ptr> fabs1() const override { + return std::make_unique(fabs(p)); + } + + std::unique_ptr> fmax1(const VPoint &x) const override { + return std::make_unique(fmax(p, dynamic_cast(x).p)); + } + std::unique_ptr> fmin1(const VPoint &x) const override { + return std::make_unique(fmin(p, dynamic_cast(x).p)); + } + std::unique_ptr> max1(const VPoint &x) const override { + return std::make_unique(max(p, dynamic_cast(x).p)); + } + std::unique_ptr> min1(const VPoint &x) const override { + return std::make_unique(min(p, dynamic_cast(x).p)); + } + + std::unique_ptr> fmax1(const T &b) const override { + return std::make_unique(fmax(p, b)); + } + std::unique_ptr> fmin1(const T &b) const override { + return std::make_unique(fmin(p, b)); + } + std::unique_ptr> max1(const T &b) const override { + return std::make_unique(max(p, b)); + } + std::unique_ptr> min1(const T &b) const override { + return std::make_unique(min(p, b)); + } + + T all1() const override { return all(p); } + T any1() const override { return any(p); } + T max_element1() const override { return max_element(p); } + T min_element1() const override { return min_element(p); } + T product1() const override { return product(p); } + T sum1() const override { return sum(p); } + + std::unique_ptr> operator==(const VPoint &x) const override { + return std::make_unique>(p == + dynamic_cast(x).p); + } + std::unique_ptr> operator!=(const VPoint &x) const override { + return std::make_unique>(p != + dynamic_cast(x).p); + } + std::unique_ptr> operator<(const VPoint &x) const override { + return std::make_unique>(p < + dynamic_cast(x).p); + } + std::unique_ptr> operator>(const VPoint &x) const override { + return std::make_unique>(p > + dynamic_cast(x).p); + } + std::unique_ptr> operator<=(const VPoint &x) const override { + return std::make_unique>(p <= + dynamic_cast(x).p); + } + std::unique_ptr> operator>=(const VPoint &x) const override { + return std::make_unique>(p >= + dynamic_cast(x).p); + } + + std::unique_ptr> operator==(const T &a) const override { + return std::make_unique>(p == a); + } + std::unique_ptr> operator!=(const T &a) const override { + return std::make_unique>(p != a); + } + std::unique_ptr> operator<(const T &a) const override { + return std::make_unique>(p < a); + } + std::unique_ptr> operator>(const T &a) const override { + return std::make_unique>(p > a); + } + std::unique_ptr> operator<=(const T &a) const override { + return std::make_unique>(p <= a); + } + std::unique_ptr> operator>=(const T &a) const override { + return std::make_unique>(p >= a); + } + + bool equal_to1(const VPoint &x) const override { + return std::equal_to>()(p, dynamic_cast(x).p); + } + std::size_t hash1() const override { return std::hash>()(p); } + bool less1(const VPoint &x) const override { + return std::less>()(p, dynamic_cast(x).p); + } + + std::ostream &output(std::ostream &os) const override { return os << p; } +}; + +template +std::unique_ptr> make_VPoint(const std::size_t D) { + switch (D) { + case 0: + return std::make_unique>(); + case 1: + return std::make_unique>(); + case 2: + return std::make_unique>(); + case 3: + return std::make_unique>(); + case 4: + return std::make_unique>(); + case 5: + return std::make_unique>(); + default: + abort(); + } +} + +template +std::unique_ptr> make_VPoint(const std::size_t D, const F &f) { + switch (D) { + case 0: + return std::make_unique>(Point::make(f)); + case 1: + return std::make_unique>(Point::make(f)); + case 2: + return std::make_unique>(Point::make(f)); + case 3: + return std::make_unique>(Point::make(f)); + case 4: + return std::make_unique>(Point::make(f)); + case 5: + return std::make_unique>(Point::make(f)); + default: + abort(); + } +} + +} // namespace detail + +/** A Point + * + * The dimension (number of component) of the Point is only known at + * run-time. @see Point + * + * Points can represent either Points or distances. Points are + * fixed-size vectors that support arithmetic operations. + */ +template class NDPoint { + template using VPoint = detail::VPoint; + + template friend class NDPoint; + friend struct std::equal_to; + friend struct std::hash; + friend struct std::less; + + std::unique_ptr> p; + + NDPoint(std::unique_ptr> p_) : p(std::move(p_)) {} + +public: + /** Component type + */ + using value_type = typename VPoint::value_type; + /** Return type of Point::size() + */ + using size_type = typename VPoint::size_type; + + /** Create an invalid Point + */ + NDPoint() : p() {} + /** Create a value-initialized Point with D components + */ + NDPoint(const size_type D) : p(detail::make_VPoint(D)) {} + + NDPoint(const NDPoint &x) : p(x.p ? x.p->copy() : nullptr) {} + NDPoint(NDPoint &&) = default; + NDPoint &operator=(const NDPoint &x) { + p = x.p ? x.p->copy() : nullptr; + return *this; + } + NDPoint &operator=(NDPoint &&) = default; + + /** Create an NDPoint from a Point + */ + template + NDPoint(const Point &p_) + : p(std::make_unique>(p_)) {} + template + NDPoint(Point &&p_) + : p(std::make_unique>(std::move(p_))) {} + /** Convert an NDPoint to a Point + * + * This only works when D == ndims(). + */ + template operator Point() const { + return Point(dynamic_cast &>(*p)); + } + + /** Create an NDPoint from a std::array + */ + template + NDPoint(const std::array &arr) + : p(std::make_unique>(arr)) {} + template + NDPoint(std::array &&arr) + : p(std::make_unique>(std::move(arr))) {} + /** Convert an NDPoint to a std::array + * + * This only works when D == ndims(). + */ + template operator std::array() const { + return std::array(dynamic_cast>(&p)); + } + /** Create an NDPoint from an initializer list + * + * Example: NDPoint{1, 2} + */ + NDPoint(std::initializer_list lst) : NDPoint(std::vector(lst)) {} + +private: + template + static std::vector convert_vector(const std::vector &vec) { + std::vector res(vec.size()); + for (std::size_t n = 0; n < res.size(); ++n) + res[n] = T(vec[n]); + return res; + } + +public: + template + NDPoint(const NDPoint &x) + : p(x.p ? NDPoint(convert_vector(std::vector(x))) : nullptr) {} + + /** Create an NDPoint from a generating function + * + * Example: NDPoint::make(3, [](int i) { return i+1; }) + * Result: NDPoint{1, 2, 3} + */ + template static NDPoint make(const size_type D, const F &f) { + return NDPoint(detail::make_VPoint(D, f)); + } + /** Create an NDPoint with each component set to the same value `a` + */ + static NDPoint pure(const size_type D, const T &val) { + return make(D, [&](size_type) { return val; }); + } + /** Create a unit NDPoint, where component dir is one, and all other + * components are zero + */ + static NDPoint unit(const size_type D, const size_type dir) { + return make(D, [&](size_type d) { return d == dir; }); + } + /** Create an NDPoint with components set to the natural number + * sequence [0, ..., D-1] + */ + static NDPoint iota(const size_type D) { + return make(D, [&](size_type d) { return d; }); + } + + /** Check whether an NDPoint is valid + * + * A valid NDPoint knows its number of dimensions, and its + * components are initialized. An invalid NDPoint does not know its + * number of dimensions and holds no data, similar to a null + * Pointer. + * + * Most other member functions cannot be called for invalid + * NDPoints. + */ + bool has_value() const { return bool(p); } + + /** Map a function over all components of one or several Points + * + * Example: + * NDPoint pi, pj; + * NDPoint pk = fmap([](auto i, auto j) { return i+j; }, pi, pj); + * This calculates the component-wise sum of pi and pj, i.e. pi + pj . + */ + template friend NDPoint fmap(const F &f, const NDPoint &x) { + return NDPoint(x.p->fmap1(f)); + } + template + friend NDPoint fmap(const F &f, const NDPoint &x, const NDPoint &y) { + return NDPoint(x.p->fmap1(f, *y.p)); + } + template + friend NDPoint fmap(const F &f, const NDPoint &x, const NDPoint &y, + const NDPoint &z) { + return NDPoint(x.p->fmap1(f, *y.p, *z.p)); + } + /** Reduce over all components of one or several NDPoints + * + * Example: + * NDPoint p; + * int s = fold([](auto r, auto i) { return r+i; }, 0, p); + * This calculates the sum of all components ("horizonal sum") of p, + * same as sum(p). + */ + template + friend T fold(const Op &op, const T &r, const NDPoint &x) { + return x.p->fold1(op, r); + } + template + friend T fold(const Op &op, const T &r, const NDPoint &x, const NDPoint &y) { + return x.p->fold1(op, r, *y.p); + } + + /** Create an NDPoint from a std::vector + */ + NDPoint(const std::vector &vec) : NDPoint(vec.size()) { + p->set_from_vector(vec); + } + /** Convert an NDPoint to a std::vector + */ + operator std::vector() const { return std::vector(p); } + + /** Number of comopnents (same as number of dimensions) + */ + size_type size() const { return p->size(); } + + /** Number of dimensions (same as number of comopnents) + */ + size_type ndims() const { return p->ndims(); } + + /** Get a component of a Point + */ + const T &operator[](const size_type d) const { return (*p)[d]; } + /** Get a component of a Point + */ + T &operator[](const size_type d) { return (*p)[d]; } + /** Remove a component from a point + * + * This reduces the dimension of a point by one. + */ + constexpr NDPoint erase(size_type dir) const { + return NDPoint(p->erase(dir)); + } + /** Add a component to a point + * + * This increases the dimension of a point by one. + */ + constexpr NDPoint insert(size_type dir, const T &a) const { + return NDPoint(p->insert(dir, a)); + } + /** Reverse the components of a point + */ + constexpr NDPoint reversed() const { return NDPoint(p->reversed()); } + + /** Apply unary plus operator element-wise + */ + friend NDPoint operator+(const NDPoint &x) { return NDPoint(+*x.p); } + /** Element-wise negate + */ + friend NDPoint operator-(const NDPoint &x) { return NDPoint(-*x.p); } + /** Element-wise bitwise not + */ + template > * = nullptr> + friend NDPoint operator~(const NDPoint &x) { + return NDPoint(~*x.p); + } + /** Element-wise logical not + */ + friend NDPoint operator!(const NDPoint &x) { + return NDPoint(!*x.p); + } + + /** Add element-wise + */ + friend NDPoint operator+(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p + *y.p); + } + /** Subtract element-wise + */ + friend NDPoint operator-(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p - *y.p); + } + /** Multiply element-wise + */ + friend NDPoint operator*(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p * *y.p); + } + /** Divide element-wise + */ + friend NDPoint operator/(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p / *y.p); + } + /** Element-wise modulo + */ + template > * = nullptr> + friend NDPoint operator%(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p % *y.p); + } + /** Element-wise bitwise and + */ + template > * = nullptr> + friend NDPoint operator&(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p & *y.p); + } + /** Element-wise bitwise or + */ + template > * = nullptr> + friend NDPoint operator|(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p | *y.p); + } + /** Element-wise bitwise exclusive or + */ + template > * = nullptr> + friend NDPoint operator^(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p ^ *y.p); + } + /** Element-wise logical and + */ + friend NDPoint operator&&(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p && *y.p); + } + /** Element-wise logical or + */ + friend NDPoint operator||(const NDPoint &x, const NDPoint &y) { + return NDPoint(*x.p || *y.p); + } + + friend NDPoint operator+(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_plus(a)); + } + friend NDPoint operator-(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_minus(a)); + } + friend NDPoint operator*(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_multiplies(a)); + } + friend NDPoint operator/(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_divides(a)); + } + template > * = nullptr> + friend NDPoint operator%(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_modulus(a)); + } + template > * = nullptr> + friend NDPoint operator&(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_bit_and(a)); + } + template > * = nullptr> + friend NDPoint operator|(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_bit_or(a)); + } + template > * = nullptr> + friend NDPoint operator^(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_bit_xor(a)); + } + friend NDPoint operator&&(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_logical_and(a)); + } + friend NDPoint operator||(const T &a, const NDPoint &y) { + return NDPoint(y.p->left_logical_or(a)); + } + + friend NDPoint operator+(const NDPoint &x, const T &b) { + return NDPoint(*x.p + b); + } + friend NDPoint operator-(const NDPoint &x, const T &b) { + return NDPoint(*x.p - b); + } + friend NDPoint operator*(const NDPoint &x, const T &b) { + return NDPoint(*x.p * b); + } + friend NDPoint operator/(const NDPoint &x, const T &b) { + return NDPoint(*x.p / b); + } + template > * = nullptr> + friend NDPoint operator%(const NDPoint &x, const T &b) { + return NDPoint(*x.p % b); + } + template > * = nullptr> + friend NDPoint operator&(const NDPoint &x, const T &b) { + return NDPoint(*x.p & b); + } + template > * = nullptr> + friend NDPoint operator|(const NDPoint &x, const T &b) { + return NDPoint(*x.p | b); + } + template > * = nullptr> + friend NDPoint operator^(const NDPoint &x, const T &b) { + return NDPoint(*x.p ^ b); + } + friend NDPoint operator&&(const NDPoint &x, const T &b) { + return NDPoint(*x.p && b); + } + friend NDPoint operator||(const NDPoint &x, const T &b) { + return NDPoint(*x.p || b); + } + + NDPoint &operator+=(const NDPoint &x) { + *p += *x.p; + return *this; + } + NDPoint &operator-=(const NDPoint &x) { + *p -= *x.p; + return *this; + } + NDPoint &operator*=(const NDPoint &x) { + *p *= *x.p; + return *this; + } + NDPoint &operator/=(const NDPoint &x) { + *p /= *x.p; + return *this; + } + template > * = nullptr> + NDPoint &operator%=(const NDPoint &x) { + *p %= *x.p; + return *this; + } + template > * = nullptr> + NDPoint &operator&=(const NDPoint &x) { + *p &= *x.p; + return *this; + } + template > * = nullptr> + NDPoint &operator|=(const NDPoint &x) { + *p |= *x.p; + return *this; + } + template > * = nullptr> + NDPoint &operator^=(const NDPoint &x) { + *p ^= *x.p; + return *this; + } + + NDPoint &operator+=(const T &a) { + *p += a; + return *this; + } + NDPoint &operator-=(const T &a) { + *p -= a; + return *this; + } + NDPoint &operator*=(const T &a) { + *p *= a; + return *this; + } + NDPoint &operator/=(const T &a) { + *p /= a; + return *this; + } + template > * = nullptr> + NDPoint &operator%=(const T &a) { + *p %= a; + return *this; + } + template > * = nullptr> + NDPoint &operator&=(const T &a) { + *p &= a; + return *this; + } + template > * = nullptr> + NDPoint &operator|=(const T &a) { + *p |= a; + return *this; + } + template > * = nullptr> + NDPoint &operator^=(const T &a) { + *p ^= a; + return *this; + } + + /** Element-wise absolute value + */ + friend NDPoint abs(const NDPoint &x) { return x.p->abs1(); } + /** Element-wise absolute value + */ + friend NDPoint fabs(const NDPoint &x) { return x.p->fabs1(); } + + /** Element-wise maximum of two points + */ + friend NDPoint fmax(const NDPoint &x, const NDPoint &y) { + return x.p->fmax1(*y.p); + } + /** Element-wise minimum of two points + */ + friend NDPoint fmin(const NDPoint &x, const NDPoint &y) { + return x.p->fmin1(*y.p); + } + /** Element-wise maximum of two points + */ + friend NDPoint max(const NDPoint &x, const NDPoint &y) { + return x.p->max1(*y.p); + } + /** Element-wise minimum of two points + */ + friend NDPoint min(const NDPoint &x, const NDPoint &y) { + return x.p->min1(*y.p); + } + + friend NDPoint fmax(const T &a, const NDPoint &y) { return y.p->fmax1(a); } + friend NDPoint fmin(const T &a, const NDPoint &y) { return y.p->fmin1(a); } + friend NDPoint max(const T &a, const NDPoint &y) { return y.p->max1(a); } + friend NDPoint min(const T &a, const NDPoint &y) { return y.p->min1(a); } + + friend NDPoint fmax(const NDPoint &x, const T &b) { return x.p->fmax1(b); } + friend NDPoint fmin(const NDPoint &x, const T &b) { return x.p->fmin1(b); } + friend NDPoint max(const NDPoint &x, const T &b) { return x.p->max1(b); } + friend NDPoint min(const NDPoint &x, const T &b) { return x.p->min1(b); } + + /** Return true if all elements are true + */ + friend T all(const NDPoint &x) { return x.p->all1(); } + /** Return true if any element is true + */ + friend T any(const NDPoint &x) { return x.p->any1(); } + /** Return maximum element + */ + friend T max_element(const NDPoint &x) { return x.p->max_element1(); } + /** Return minimum element + */ + friend T min_element(const NDPoint &x) { return x.p->min_element1(); } + /** Product of all elements + */ + friend T product(const NDPoint &x) { return x.p->product1(); } + /** Sum of all elements + */ + friend T sum(const NDPoint &x) { return x.p->sum1(); } + + /** Pointwise comparison + */ + friend NDPoint operator==(const NDPoint &x, const NDPoint &y) { + return *x.p == *y.p; + } + /** Pointwise comparison + */ + friend NDPoint operator!=(const NDPoint &x, const NDPoint &y) { + return *x.p != *y.p; + } + /** Pointwise comparison + */ + friend NDPoint operator<(const NDPoint &x, const NDPoint &y) { + return *x.p < *y.p; + } + /** Pointwise comparison + */ + friend NDPoint operator>(const NDPoint &x, const NDPoint &y) { + return *x.p > *y.p; + } + /** Pointwise comparison + */ + friend NDPoint operator<=(const NDPoint &x, const NDPoint &y) { + return *x.p <= *y.p; + } + /** Pointwise comparison + */ + friend NDPoint operator>=(const NDPoint &x, const NDPoint &y) { + return *x.p >= *y.p; + } + + friend NDPoint operator==(const NDPoint &x, const T &b) { + return *x.p == b; + } + friend NDPoint operator!=(const NDPoint &x, const T &b) { + return *x.p != b; + } + friend NDPoint operator<(const NDPoint &x, const T &b) { + return *x.p < b; + } + friend NDPoint operator>(const NDPoint &x, const T &b) { + return *x.p > b; + } + friend NDPoint operator<=(const NDPoint &x, const T &b) { + return *x.p <= b; + } + friend NDPoint operator>=(const NDPoint &x, const T &b) { + return *x.p >= b; + } + + /** Output an NDoint + */ + friend std::ostream &operator<<(std::ostream &os, const NDPoint &x) { + if (x.p) + x.p->output(os); + else + os << "[INVALID]"; + return os; + } +}; + +} // namespace Regions +} // namespace openPMD + +namespace std { + +/** Specialization of `equal_to` for NDPoint + */ +template struct equal_to> { + constexpr bool operator()(const openPMD::Regions::NDPoint &x, + const openPMD::Regions::NDPoint &y) const { + if (!x.has_value() && !y.has_value()) + return true; + if (!x.has_value() || !y.has_value()) + return false; + return x.p->equal_to1(*y.p); + } +}; + +/** Specialization of `hash` for NDPoint + */ +template struct hash> { + constexpr size_t operator()(const openPMD::Regions::NDPoint &x) const { + if (!x.has_value()) + return size_t(0xf458b18eca40aef1ULL); + return x.p->hash1(); + } +}; + +/** Specialization of `less` for NDPoint + */ +template struct less> { + constexpr bool operator()(const openPMD::Regions::NDPoint &x, + const openPMD::Regions::NDPoint &y) const { + if (!x.has_value()) + return true; + if (!y.has_value()) + return false; + return x.p->less1(*y.p); + } +}; + +} // namespace std diff --git a/include/openPMD/regions/NDRegion.hpp b/include/openPMD/regions/NDRegion.hpp new file mode 100644 index 0000000000..78b963d793 --- /dev/null +++ b/include/openPMD/regions/NDRegion.hpp @@ -0,0 +1,664 @@ +#pragma once + +#include "Helpers.hpp" +#include "NDBox.hpp" +#include "NDPoint.hpp" +#include "Region.hpp" + +namespace openPMD { +namespace Regions { + +template class NDRegion; + +namespace detail { + +// Abstract base helper class + +template class VRegion { +public: + using value_type = T; + using size_type = typename VPoint::size_type; + + virtual std::unique_ptr copy() const = 0; + + virtual operator std::vector>() const = 0; + + virtual ~VRegion() {} + + virtual size_type ndims() const = 0; + virtual bool empty() const = 0; + + virtual size_type size() const = 0; + + virtual size_type nboxes() const = 0; + + virtual std::unique_ptr operator>>(const NDPoint &d) const = 0; + virtual std::unique_ptr operator<<(const NDPoint &d) const = 0; + virtual VRegion &operator>>=(const NDPoint &d) = 0; + virtual VRegion &operator<<=(const NDPoint &d) = 0; + virtual std::unique_ptr operator*(const NDPoint &d) const = 0; + virtual VRegion &operator*=(const NDPoint &d) = 0; + virtual std::unique_ptr grown(const NDPoint &dlo, + const NDPoint &dup) const = 0; + virtual std::unique_ptr grown(const NDPoint &d) const = 0; + virtual std::unique_ptr grown(const T &n) const = 0; + virtual std::unique_ptr shrunk(const NDPoint &dlo, + const NDPoint &dup) const = 0; + virtual std::unique_ptr shrunk(const NDPoint &d) const = 0; + virtual std::unique_ptr shrunk(const T &n) const = 0; + + virtual bool operator==(const VRegion ®ion2) const = 0; + virtual bool operator!=(const VRegion ®ion2) const = 0; + virtual bool operator==(const NDBox &box2) const = 0; + virtual bool operator!=(const NDBox &box2) const = 0; + + virtual NDBox bounding_box1() const = 0; + + virtual std::unique_ptr operator&(const VRegion ®ion2) const = 0; + virtual std::unique_ptr operator|(const VRegion ®ion2) const = 0; + virtual std::unique_ptr operator^(const VRegion ®ion2) const = 0; + virtual std::unique_ptr operator-(const VRegion ®ion2) const = 0; + + virtual VRegion &operator&=(const VRegion ®ion) = 0; + virtual VRegion &operator|=(const VRegion ®ion) = 0; + virtual VRegion &operator^=(const VRegion ®ion) = 0; + virtual VRegion &operator-=(const VRegion ®ion) = 0; + + virtual std::unique_ptr + intersection1(const VRegion ®ion2) const = 0; + virtual std::unique_ptr setunion1(const VRegion ®ion2) const = 0; + virtual std::unique_ptr + symmetric_difference1(const VRegion ®ion2) const = 0; + virtual std::unique_ptr + difference1(const VRegion ®ion2) const = 0; + + virtual bool contains(const NDPoint &p) const = 0; + virtual bool isdisjoint1(const VRegion ®ion2) const = 0; + + virtual bool operator<=(const VRegion ®ion2) const = 0; + virtual bool operator>=(const VRegion ®ion2) const = 0; + virtual bool operator<(const VRegion ®ion2) const = 0; + virtual bool operator>(const VRegion ®ion2) const = 0; + + virtual bool is_subset_of(const VRegion ®ion) const = 0; + virtual bool is_superset_of(const VRegion ®ion) const = 0; + virtual bool is_strict_subset_of(const VRegion ®ion) const = 0; + virtual bool is_strict_superset_of(const VRegion ®ion) const = 0; + + virtual bool equal_to1(const VRegion &x) const = 0; + virtual bool less1(const VRegion &x) const = 0; + + virtual std::ostream &output(std::ostream &os) const = 0; +}; + +// Helper class wrapping Region + +template class WRegion final : public VRegion { + Region r; + +public: + using typename VRegion::value_type; + using typename VRegion::size_type; + + WRegion() : r{} {} + + WRegion(const WRegion &x) = default; + WRegion(WRegion &&) = default; + WRegion &operator=(const WRegion &) = default; + WRegion &operator=(WRegion &&) = default; + + WRegion(const Region &r_) : r(r_) {} + WRegion(Region &&r_) : r(std::move(r_)) {} + operator Region() const { return r; } + + WRegion(const Point &p) : r(p) {} + WRegion(const Box &b) : r(b) {} + WRegion(const std::vector> &boxes) : r(boxes) {} + + WRegion(const NDPoint &p) : r(Point(p)) {} + WRegion(const NDBox &b) : r(Box(b)) {} + WRegion(const std::vector> &boxes) + : r(helpers::convert_vector>(boxes)) {} + + operator std::vector>() const { return std::vector>(r); } + operator std::vector>() const override { + return helpers::convert_vector>(std::vector>(r)); + } + + std::unique_ptr> copy() const override { + return std::make_unique(*this); + } + + ~WRegion() override {} + + size_type ndims() const override { return r.ndims(); } + bool empty() const override { return r.empty(); } + + size_type size() const override { return r.size(); } + + size_type nboxes() const override { return r.nboxes(); } + + std::unique_ptr> operator>>(const NDPoint &d) const override { + return std::make_unique(r >> Point(d)); + } + std::unique_ptr> operator<<(const NDPoint &d) const override { + return std::make_unique(r << Point(d)); + } + VRegion &operator>>=(const NDPoint &d) override { + r >>= Point(d); + return *this; + } + VRegion &operator<<=(const NDPoint &d) override { + r <<= Point(d); + return *this; + } + std::unique_ptr> operator*(const NDPoint &d) const override { + return std::make_unique(r * Point(d)); + } + VRegion &operator*=(const NDPoint &d) override { + r *= Point(d); + return *this; + } + std::unique_ptr> grown(const NDPoint &dlo, + const NDPoint &dup) const override { + return std::make_unique( + r.grown(Point(dlo), Point(dup))); + } + std::unique_ptr> grown(const NDPoint &d) const override { + return std::make_unique(r.grown(Point(d))); + } + std::unique_ptr> grown(const T &n) const override { + return std::make_unique(r.grown(n)); + } + std::unique_ptr> shrunk(const NDPoint &dlo, + const NDPoint &dup) const override { + return std::make_unique( + r.shrunk(Point(dlo), Point(dup))); + } + std::unique_ptr> shrunk(const NDPoint &d) const override { + return std::make_unique(r.shrunk(Point(d))); + } + std::unique_ptr> shrunk(const T &n) const override { + return std::make_unique(r.shrunk(n)); + } + + bool operator==(const VRegion ®ion2) const override { + return r == dynamic_cast(region2); + } + bool operator!=(const VRegion ®ion2) const override { + return r != dynamic_cast(region2); + } + bool operator==(const NDBox &box2) const override { + return r == Box(box2); + } + bool operator!=(const NDBox &box2) const override { + return r != Box(box2); + } + + NDBox bounding_box1() const override { return NDBox(bounding_box(r)); } + + std::unique_ptr> + operator&(const VRegion ®ion2) const override { + return std::make_unique(r & + dynamic_cast(region2)); + } + std::unique_ptr> + operator|(const VRegion ®ion2) const override { + return std::make_unique(r | + dynamic_cast(region2)); + } + std::unique_ptr> + operator^(const VRegion ®ion2) const override { + return std::make_unique(r ^ + dynamic_cast(region2)); + } + std::unique_ptr> + operator-(const VRegion ®ion2) const override { + return std::make_unique(r - + dynamic_cast(region2)); + } + + VRegion &operator&=(const VRegion ®ion) override { + r &= dynamic_cast(region); + return *this; + } + VRegion &operator|=(const VRegion ®ion) override { + r |= dynamic_cast(region); + return *this; + } + VRegion &operator^=(const VRegion ®ion) override { + r ^= dynamic_cast(region); + return *this; + } + VRegion &operator-=(const VRegion ®ion) override { + r -= dynamic_cast(region); + return *this; + } + + std::unique_ptr> + intersection1(const VRegion ®ion2) const override { + return std::make_unique( + intersection(r, dynamic_cast(region2))); + } + std::unique_ptr> + setunion1(const VRegion ®ion2) const override { + return std::make_unique( + setunion(r, dynamic_cast(region2))); + } + std::unique_ptr> + symmetric_difference1(const VRegion ®ion2) const override { + return std::make_unique( + symmetric_difference(r, dynamic_cast(region2))); + } + std::unique_ptr> + difference1(const VRegion ®ion2) const override { + return std::make_unique( + difference(r, dynamic_cast(region2))); + } + + bool contains(const NDPoint &p) const override { + return r.contains(Point(p)); + } + bool isdisjoint1(const VRegion ®ion2) const override { + return isdisjoint(r, dynamic_cast(region2)); + } + + bool operator<=(const VRegion ®ion2) const override { + return r <= dynamic_cast(region2); + } + bool operator>=(const VRegion ®ion2) const override { + return r >= dynamic_cast(region2); + } + bool operator<(const VRegion ®ion2) const override { + return r < dynamic_cast(region2); + } + bool operator>(const VRegion ®ion2) const override { + return r > dynamic_cast(region2); + } + + bool is_subset_of(const VRegion ®ion) const override { + return r.is_subset_of(dynamic_cast(region)); + } + bool is_superset_of(const VRegion ®ion) const override { + return r.is_superset_of(dynamic_cast(region)); + } + bool is_strict_subset_of(const VRegion ®ion) const override { + return r.is_strict_subset_of(dynamic_cast(region)); + } + bool is_strict_superset_of(const VRegion ®ion) const override { + return r.is_strict_superset_of(dynamic_cast(region)); + } + + bool equal_to1(const VRegion &x) const override { + return std::equal_to>()(r, dynamic_cast(x).r); + } + bool less1(const VRegion &x) const override { + return std::less>()(r, dynamic_cast(x).r); + } + + std::ostream &output(std::ostream &os) const override { return os << r; } +}; + +template +std::unique_ptr> make_VRegion(const std::size_t D, + const Args &...args) { + switch (D) { + case 0: + return std::make_unique>(args...); + case 1: + return std::make_unique>(args...); + case 2: + return std::make_unique>(args...); + case 3: + return std::make_unique>(args...); + case 4: + return std::make_unique>(args...); + case 5: + return std::make_unique>(args...); + default: + abort(); + } +} + +} // namespace detail + +/** A Region + * + * The dimension of the NDRegion is only known at run-time. @see + * Region + * + * A region is an arbitrarily shaped set of points. The internal + * representation is likely based on boxes, and is thus most efficient + * if the region has many axis-aligned boundaries. + */ +template class NDRegion { + template using VRegion = detail::VRegion; + + template friend class NDRegion; + friend struct std::equal_to; + friend struct std::less; + + std::unique_ptr> r; + + NDRegion(std::unique_ptr> r_) : r(std::move(r_)) {} + +public: + /** Component type + */ + using value_type = typename VRegion::value_type; + /** Return type of Region::size() + */ + using size_type = typename VRegion::size_type; + + /** Create an invalid Region + */ + NDRegion() : r() {} + /** Create an empty Region in D dimensions + */ + NDRegion(const size_type D) : r(detail::make_VRegion(D)) {} + + NDRegion(const NDRegion &x) : r(x.r ? x.r->copy() : nullptr) {} + NDRegion(NDRegion &&) = default; + NDRegion &operator=(const NDRegion &x) { + r = x.r ? x.r->copy() : nullptr; + return *this; + } + NDRegion &operator=(NDRegion &&) = default; + + /** Create an NDRegion from a Region + */ + template + NDRegion(const Region &r_) + : r(std::make_unique>(r_)) {} + template + NDRegion(Region &&r_) + : r(std::make_unique>(std::move(r_))) {} + /** Convert an NDRegion to a Region + * + * This only works when D == ndims(). + */ + template operator Region() const { + return Region(dynamic_cast &>(*r)); + } + + template + NDRegion(const Point &p) + : r(std::make_unique>(detail::WRegion(p))) {} + template + NDRegion(const Box &b) + : r(std::make_unique>(detail::WRegion(b))) {} + template + NDRegion(const std::vector> &boxes) + : r(std::make_unique>(detail::WRegion(boxes))) {} + template operator std::vector>() const { + return std::vector>( + dynamic_cast &>(*r)); + } + + /** Create Region containing a single NDPoint + */ + NDRegion(const NDPoint &p) : r(detail::make_VRegion(p.ndims(), p)) {} + /** Create Region containina a single NDBox + */ + NDRegion(const NDBox &b) : r(detail::make_VRegion(b.ndims(), b)) {} + /** Create a Region from a vector of NDBoxes + */ + NDRegion(const size_type D, const std::vector> &boxes) + : r(detail::make_VRegion(D, boxes)) {} + operator std::vector>() const { return std::vector>(*r); } + + /** Check whether an NDRegion is valid + * + * A valid NDRegion knows its number of dimensions, and its components + * are initialized. An invalid NDRegion does not know its number of + * dimensions and holds no data, similar to a null pointer. + * + * Most other member functions cannot not be called for invalid + * NDRegions. + */ + bool has_value() const { return bool(r); } + + /** Number of dimensions + */ + size_type ndims() const { return r->ndims(); } + /** Whether the Region is empty + */ + bool empty() const { return r->empty(); } + + /** Size, the total number of contained points + */ + size_type size() const { return r->size(); } + + /** A measure of the number of vertices defining the Region + */ + size_type nboxes() const { return r->nboxes(); } + + /** Shift an NDRegion to the right (upwards). The shift can be negative, which + * shifts left. + */ + NDRegion operator>>(const NDPoint &d) const { return NDRegion(*r >> d); } + /** Shift an NDRegion to the left (downwards). The shift can be negative, + * which shifts right. + */ + NDRegion operator<<(const NDPoint &d) const { return NDRegion(*r << d); } + NDRegion &operator>>=(const NDPoint &d) { + *r >>= d; + return *this; + } + NDRegion &operator<<=(const NDPoint &d) { + *r <<= d; + return *this; + } + /** Scale an NDRegion + */ + NDRegion operator*(const NDPoint &d) const { return NDRegion(*r * d); } + NDRegion &operator*=(const NDPoint &d) { + *r *= d; + return *this; + } + + /** Grow an NDRegion by given amounts in each direction. + * + * The growth can be negative, which shrinks the NDRegion. If an NDRegion is + * shrunk too much it becomes empty. Growing an empty NDRegion always + * results in an empty NDRegion. + */ + NDRegion grown(const NDPoint &dlo, const NDPoint &dup) const { + return NDRegion(r->grown(dlo, dup)); + } + NDRegion grown(const NDPoint &d) const { return NDRegion(r->grown(d)); } + NDRegion grown(const T &n) const { return NDRegion(r->grown(n)); } + /** Shrink an NDRegion by given amounts in each direction. + * + * The shrinkage can be negative, which shrinks the NDRegion. If a Region is + * shrunk too much it becomes empty. Growing an empty NDRegion always + * results in an empty NDRegion. + */ + NDRegion shrunk(const NDPoint &dlo, const NDPoint &dup) const { + return NDRegion(r->shrunk(dlo, dup)); + } + NDRegion shrunk(const NDPoint &d) const { return NDRegion(r->shrunk(d)); } + NDRegion shrunk(const T &n) const { return NDRegion(r->shrunk(n)); } + + /** Compare two Regions + */ + friend bool operator==(const NDRegion ®ion1, const NDRegion ®ion2) { + return *region1.r == *region2.r; + } + /** Compare two Regions + */ + friend bool operator!=(const NDRegion ®ion1, const NDRegion ®ion2) { + return *region1.r != *region2.r; + } + friend bool operator==(const NDRegion ®ion1, const NDBox &box2) { + return *region1.r == box2; + } + friend bool operator!=(const NDRegion ®ion1, const NDBox &box2) { + return *region1.r != box2; + } + friend bool operator==(const NDBox &box1, const NDRegion ®ion2) { + return *region2.r == box1; + } + friend bool operator!=(const NDBox &box1, const NDRegion ®ion2) { + return *region2.r != box1; + } + + /** Calculate the bounding box of an NDRegion. This is the smallest NDBox that + * contains the NDRegion. + */ + friend NDBox bounding_box(const NDRegion ®ion) { + return region.r->bounding_box1(); + } + + /** Set intersection of two NDRegions + */ + friend NDRegion operator&(const NDRegion ®ion1, const NDRegion ®ion2) { + return NDRegion(*region1.r & *region2.r); + } + /** Set union of two NDRegions + */ + friend NDRegion operator|(const NDRegion ®ion1, const NDRegion ®ion2) { + return NDRegion(*region1.r | *region2.r); + } + /** Symmetric difference of two NDRegions + */ + friend NDRegion operator^(const NDRegion ®ion1, const NDRegion ®ion2) { + return NDRegion(*region1.r ^ *region2.r); + } + /** Set difference of two NDRegions + */ + friend NDRegion operator-(const NDRegion ®ion1, const NDRegion ®ion2) { + return NDRegion(*region1.r - *region2.r); + } + + NDRegion &operator&=(const NDRegion ®ion) { + *r &= *region.r; + return *this; + } + NDRegion &operator|=(const NDRegion ®ion) { + *r |= *region.r; + return *this; + } + NDRegion &operator^=(const NDRegion ®ion) { + *r ^= *region.r; + return *this; + } + NDRegion &operator-=(const NDRegion ®ion) { + *r -= *region.r; + return *this; + } + + /** Set intersection of two NDRegions + */ + friend NDRegion intersection(const NDRegion ®ion1, + const NDRegion ®ion2) { + return NDRegion(intersection(*region1.r, *region2.r)); + } + /** Set union of two NDRegions + */ + friend NDRegion setunion(const NDRegion ®ion1, const NDRegion ®ion2) { + return NDRegion(setunion(*region1.r, *region2.r)); + } + /** Symmetric difference of two NDRegions + */ + friend NDRegion symmetric_difference(const NDRegion ®ion1, + const NDRegion ®ion2) { + return NDRegion(symmetric_difference(*region1.r, *region2.r)); + } + /** Set difference of two NDRegions + */ + friend NDRegion difference(const NDRegion ®ion1, const NDRegion ®ion2) { + return NDRegion(difference(*region1.r, *region2.r)); + } + + /** Whether an NDRegion contains a Point + */ + bool contains(const NDPoint &p) const { return r->contains(p); } + /** Whether two NDRegions are disjoint, i.e. whether they have no Point in + * common + */ + friend bool isdisjoint(const NDRegion ®ion1, const NDRegion ®ion2) { + return region1.r->isdisjoint1(*region2.r); + } + + /** is subset of + */ + friend bool operator<=(const NDRegion ®ion1, const NDRegion ®ion2) { + return *region1.r <= *region2.r; + } + /** is superset of + */ + friend bool operator>=(const NDRegion ®ion1, const NDRegion ®ion2) { + return *region1.r >= *region2.r; + } + /** is strict subset of + */ + friend bool operator<(const NDRegion ®ion1, const NDRegion ®ion2) { + return *region1.r < *region2.r; + } + /** is strict superset of + */ + friend bool operator>(const NDRegion ®ion1, const NDRegion ®ion2) { + return *region1.r > *region2.r; + } + + /** is subset of + */ + bool is_subset_of(const NDRegion ®ion) const { + return r->is_subset_of(*region.r); + } + /** is superset of + */ + bool is_superset_of(const NDRegion ®ion) const { + return r->is_superset_of(*region.r); + } + /** is strict subset of + */ + bool is_strict_subset_of(const NDRegion ®ion) const { + return r->is_strict_subset_of(*region.r); + } + /** is strict superset of + */ + bool is_strict_superset_of(const NDRegion ®ion) const { + return r->is_strict_superset_of(*region.r); + } + + /** Output an NDegion + */ + friend std::ostream &operator<<(std::ostream &os, const NDRegion &x) { + if (x.r) + x.r->output(os); + else + os << "{INVALID}"; + return os; + } +}; + +} // namespace Regions +} // namespace openPMD + +namespace std { + +template struct equal_to> { + constexpr bool operator()(const openPMD::Regions::NDRegion &x, + const openPMD::Regions::NDRegion &y) const { + if (!x.has_value() && !y.has_value()) + return true; + if (!x.has_value() || !y.has_value()) + return false; + return x.r->equal_to1(*y.r); + } +}; + +template struct less> { + constexpr bool operator()(const openPMD::Regions::NDRegion &x, + const openPMD::Regions::NDRegion &y) const { + if (!x.has_value()) + return true; + if (!y.has_value()) + return false; + return x.r->less1(*y.r); + } +}; + +} // namespace std diff --git a/include/openPMD/regions/Point.hpp b/include/openPMD/regions/Point.hpp new file mode 100644 index 0000000000..70e2dc449d --- /dev/null +++ b/include/openPMD/regions/Point.hpp @@ -0,0 +1,656 @@ +#pragma once + +#include "Helpers.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace openPMD { +namespace Regions { + +/** A D-dimensional point + * + * The dimension D needs to be known at compile time. @see NDPoint + * + * Points can represent either points or distances. Points are + * fixed-size vectors that support arithmetic operations. + */ +template class Point { + std::array elts; + + friend struct std::equal_to; + friend struct std::hash; + friend struct std::less; + +public: + /** Component type + */ + using value_type = T; + /** Return type of Point::size() + */ + using size_type = std::size_t; + + /** Create a value-initialized Point + * + * For most types, this initializes all components to zero. + */ + constexpr Point() : elts{} {} // value initialization + + Point(const Point &) = default; + Point(Point &&) = default; + Point &operator=(const Point &) = default; + Point &operator=(Point &&) = default; + + /** Loop over the natural number sequence [0, ..., D-1], and evaluate f for + * each number + */ + template static constexpr void loop(const F &f) { + for (std::size_t d = 0; d < D; ++d) + f(d); + } + /** Create a new Point by applying a function to the natural number sequence + * [0, ..., D-1] + */ + template static constexpr Point make(const F &f) { + return helpers::construct_array(f); + } + + /** Create a Point with each component set to the same value `a` + */ + static constexpr Point pure(const T &a) { + return make([&](size_type) { return a; }); + } + + /** Create a unit Point, where component dir is one, and all other components + * are zero + */ + static constexpr Point unit(size_type dir) { + return make([&](size_type d) { return d == dir; }); + } + + /** Create a Point with components set to the natural number + * sequence [0, ..., D-1] + */ + static constexpr Point iota() { + return Point::make([](size_type d) { return d; }); + } + + /** Map a function over all components of one or several Points + * + * Example: + * Point pi, pj; + * Point pk = fmap([](auto i, auto j) { return i+j; }, pi, pj); + * This calculates the component-wise sum of pi and pj, i.e. pi + pj . + */ + template >>> + friend constexpr Point fmap(const F &f, const Point &x, + const Point &...args) { + return Point::make( + [&](size_type d) { return f(x.elts[d], args.elts[d]...); }); + } + /** Map a function over all components of one or several Points, returning + * void + * + * Example: + * Point pi; + * fmap_([](auto& i) { return i*=2; }, pi); + * This doubles each component of pi, the same as pi *= 2 . + */ + template + friend constexpr void fmap_(const F &f, const Point &x, + const Point &...args) { + loop([&](size_type d) { f(x.elts[d], args.elts[d]...); }); + } + + /** Reduce over all components of one or several Points + * + * Example: + * Point pi; + * int s = fold([](auto r, auto i) { return r+i; }, 0, pi); + * This calculates the sum of all components ("horizonal sum") of pi, + * same as sum(pi). + */ + template >>> + friend constexpr R fold(const Op &op, R r, const Point &x, + const Point &...args) { + loop([&](int d) { r = op(r, x[d], args[d]...); }); + return r; + } + + /** Create a Point from a Point with different component type + */ + template + constexpr Point(const Point &x) + : elts(fmap([](const U &a) { return T(a); }, x)) {} + + /** Create a Point from Pointers to first and one past the last element + */ + constexpr Point(const T *begin, const T *end [[maybe_unused]]) + : elts((assert(begin + D == end), + make([&](size_type d) { return begin[d]; }))) {} + /** Create a Point from an initializer list + * + * Example: Point{1,2,3} + */ + constexpr Point(std::initializer_list lst) + : Point(lst.begin(), lst.end()) {} + /** Create a Point from a C-style array + */ + template * = nullptr> + constexpr Point(const T (&arr)[DD]) : elts(&arr[0], &arr[D]) {} + /** Create a Point from a std::array + */ + constexpr Point(const std::array &arr) : elts(arr) {} + /** Create a Point from a std::vector + */ + template > * = nullptr> + constexpr Point(const std::vector &vec) + : Point(&*vec.begin(), &*vec.end()) {} + constexpr Point(const std::vector &vec) + : Point((assert(vec.size() == D), + make([&](size_type d) { return vec[d]; }))) {} + + /** Convert a Point to a std::array + */ + constexpr operator std::array() const { return elts; } + /** Convert a Point to a std::vector + */ + explicit constexpr operator std::vector() const { + return std::vector(elts.begin(), elts.end()); + } + + /** Number of components (same as number of dimensions) + */ + constexpr size_type size() const { return D; } + + /** Number of dimensions (same as number of components) + */ + constexpr size_type ndims() const { return D; } + + /** Get a component of a Point + */ + constexpr const T &operator[](const size_type d) const { return elts[d]; } + /** Get a component of a Point + */ + constexpr T &operator[](const size_type d) { return elts[d]; } + + /** Remove a component from a Point + * + * This reduces the dimension of a Point by one. + */ + constexpr Point 0 ? D - 1 : 0)> erase(size_type dir) const { + assert(dir >= 0 && dir < D); + return Point 0 ? D - 1 : 0)>::make( + [&](size_type d) { return d < dir ? (*this)[d] : (*this)[d + 1]; }); + } + /** Add a component to a Point + * + * This increases the dimension of a Point by one. + */ + constexpr Point insert(size_type dir, const T &a) const { + assert(dir >= 0 && dir <= D); + return Point::make([&](size_type d) { + return d < dir ? (*this)[d] : d == dir ? a : (*this)[d - 1]; + }); + } + + /** Reverse the components of a Point + */ + constexpr Point reversed() const { + return Point::make([&](size_type d) { return (*this)[D - 1 - d]; }); + } + + /** Apply unary plus operator element-wise + */ + friend constexpr Point operator+(const Point &x) { + return fmap([](const T &a) { return +a; }, x); + } + /** Element-wise negate + */ + friend constexpr Point operator-(const Point &x) { + return fmap([](const T &a) { return -a; }, x); + } + /** Element-wise bitwise not + */ + template > * = nullptr> + friend constexpr Point operator~(const Point &x) { + if constexpr (std::is_same_v) + // Special case to avoid compiler warnings + return fmap([](const T &) { return true; }, x); + else + return fmap([](const T &a) { return ~a; }, x); + } + /** Element-wise logical not + */ + friend constexpr Point operator!(const Point &x) { + return fmap([](const T &a) { return !a; }, x); + } + + /** Add element-wise + */ + friend constexpr Point operator+(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a + b; }, x, y); + } + /** Subtract element-wise + */ + friend constexpr Point operator-(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a - b; }, x, y); + } + /** Multiply element-wise + */ + friend constexpr Point operator*(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a * b; }, x, y); + } + /** Divide element-wise + */ + friend constexpr Point operator/(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a / b; }, x, y); + } + /** Element-wise modulo + */ + template > * = nullptr> + friend constexpr Point operator%(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a % b; }, x, y); + } + /** Element-wise bitwise and + */ + template > * = nullptr> + friend constexpr Point operator&(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a & b; }, x, y); + } + /** Element-wise bitwise or + */ + template > * = nullptr> + friend constexpr Point operator|(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a | b; }, x, y); + } + /** Element-wise bitwise exclusive or + */ + template > * = nullptr> + friend constexpr Point operator^(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a ^ b; }, x, y); + } + /** Element-wise logical and + */ + friend constexpr Point operator&&(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a && b; }, x, y); + } + /** Element-wise logical or + */ + friend constexpr Point operator||(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a || b; }, x, y); + } + + friend constexpr Point operator+(const T &a, const Point &y) { + return fmap([&](const T &b) { return a + b; }, y); + } + friend constexpr Point operator-(const T &a, const Point &y) { + return fmap([&](const T &b) { return a - b; }, y); + } + friend constexpr Point operator*(const T &a, const Point &y) { + return fmap([&](const T &b) { return a * b; }, y); + } + friend constexpr Point operator/(const T &a, const Point &y) { + return fmap([&](const T &b) { return a / b; }, y); + } + template > * = nullptr> + friend constexpr Point operator%(const T &a, const Point &y) { + return fmap([&](const T &b) { return a % b; }, y); + } + template > * = nullptr> + friend constexpr Point operator&(const T &a, const Point &y) { + return fmap([&](const T &b) { return a & b; }, y); + } + template > * = nullptr> + friend constexpr Point operator|(const T &a, const Point &y) { + return fmap([&](const T &b) { return a | b; }, y); + } + template > * = nullptr> + friend constexpr Point operator^(const T &a, const Point &y) { + return fmap([&](const T &b) { return a ^ b; }, y); + } + friend constexpr Point operator&&(const T &a, const Point &y) { + return fmap([&](const T &b) { return a && b; }, y); + } + friend constexpr Point operator||(const T &a, const Point &y) { + return fmap([&](const T &b) { return a || b; }, y); + } + + friend constexpr Point operator+(const Point &x, const T &b) { + return fmap([&](const T &a) { return a + b; }, x); + } + friend constexpr Point operator-(const Point &x, const T &b) { + return fmap([&](const T &a) { return a - b; }, x); + } + friend constexpr Point operator*(const Point &x, const T &b) { + return fmap([&](const T &a) { return a * b; }, x); + } + friend constexpr Point operator/(const Point &x, const T &b) { + return fmap([&](const T &a) { return a / b; }, x); + } + template > * = nullptr> + friend constexpr Point operator%(const Point &x, const T &b) { + return fmap([&](const T &a) { return a % b; }, x); + } + template > * = nullptr> + friend constexpr Point operator&(const Point &x, const T &b) { + return fmap([&](const T &a) { return a & b; }, x); + } + template > * = nullptr> + friend constexpr Point operator|(const Point &x, const T &b) { + return fmap([&](const T &a) { return a | b; }, x); + } + template > * = nullptr> + friend constexpr Point operator^(const Point &x, const T &b) { + return fmap([&](const T &a) { return a ^ b; }, x); + } + friend constexpr Point operator&&(const Point &x, const T &b) { + return fmap([&](const T &a) { return a && b; }, x); + } + friend constexpr Point operator||(const Point &x, const T &b) { + return fmap([&](const T &a) { return a || b; }, x); + } + + constexpr Point &operator+=(const Point &x) { return *this = *this + x; } + constexpr Point &operator-=(const Point &x) { return *this = *this - x; } + constexpr Point &operator*=(const Point &x) { return *this = *this * x; } + constexpr Point &operator/=(const Point &x) { return *this = *this / x; } + template > * = nullptr> + constexpr Point &operator%=(const Point &x) { + return *this = *this % x; + } + template > * = nullptr> + constexpr Point &operator&=(const Point &x) { + return *this = *this & x; + } + template > * = nullptr> + constexpr Point &operator|=(const Point &x) { + return *this = *this | x; + } + template > * = nullptr> + constexpr Point &operator^=(const Point &x) { + return *this = *this ^ x; + } + + constexpr Point &operator+=(const T &a) { return *this = *this + a; } + constexpr Point &operator-=(const T &a) { return *this = *this - a; } + constexpr Point &operator*=(const T &a) { return *this = *this * a; } + constexpr Point &operator/=(const T &a) { return *this = *this / a; } + template > * = nullptr> + constexpr Point &operator%=(const T &a) { + return *this = *this % a; + } + template > * = nullptr> + constexpr Point &operator&=(const T &a) { + return *this = *this & a; + } + template > * = nullptr> + constexpr Point &operator|=(const T &a) { + return *this = *this | a; + } + template > * = nullptr> + constexpr Point &operator^=(const T &a) { + return *this = *this ^ a; + } + + /** Element-wise absolute value + */ + friend constexpr Point abs(const Point &x) { + using std::abs; + return fmap( + [&](const T &a) { + if constexpr (std::is_same_v) + return a; + else + return abs(a); + }, + x); + } + /** Element-wise absolute value + */ + friend constexpr Point fabs(const Point &x) { + using std::fabs; + return fmap([&](const T &a) { return fabs(a); }, x); + } + + /** Element-wise maximum of two points + */ + friend constexpr Point fmax(const Point &x, const Point &y) { + using std::fmax; + return fmap([](const T &a, const T &b) { return fmax(a, b); }, x, y); + } + /** Element-wise minimum of two points + */ + friend constexpr Point fmin(const Point &x, const Point &y) { + using std::fmin; + return fmap([](const T &a, const T &b) { return fmin(a, b); }, x, y); + } + /** Element-wise maximum of two points + */ + friend constexpr Point max(const Point &x, const Point &y) { + using std::max; + return fmap([](const T &a, const T &b) { return max(a, b); }, x, y); + } + /** Element-wise minimum of two points + */ + friend constexpr Point min(const Point &x, const Point &y) { + using std::min; + return fmap([](const T &a, const T &b) { return min(a, b); }, x, y); + } + + friend constexpr Point fmax(const T &a, const Point &y) { + using std::fmax; + return fmap([&](const T &b) { return fmax(a, b); }, y); + } + friend constexpr Point fmin(const T &a, const Point &y) { + using std::fmin; + return fmap([&](const T &b) { return fmin(a, b); }, y); + } + friend constexpr Point max(const T &a, const Point &y) { + using std::max; + return fmap([&](const T &b) { return max(a, b); }, y); + } + friend constexpr Point min(const T &a, const Point &y) { + using std::min; + return fmap([&](const T &b) { return min(a, b); }, y); + } + + friend constexpr Point fmax(const Point &x, const T &b) { + using std::fmax; + return fmap([&](const T &a) { return fmax(a, b); }, x); + } + friend constexpr Point fmin(const Point &x, const T &b) { + using std::fmin; + return fmap([&](const T &a) { return fmin(a, b); }, x); + } + friend constexpr Point max(const Point &x, const T &b) { + using std::max; + return fmap([&](const T &a) { return max(a, b); }, x); + } + friend constexpr Point min(const Point &x, const T &b) { + using std::min; + return fmap([&](const T &a) { return min(a, b); }, x); + } + + /** Return true if all elements are true + */ + friend constexpr bool all(const Point &x) { + return fold([](bool r, const T &a) { return r && a; }, true, x); + } + /** Return true if any element is true + */ + friend constexpr bool any(const Point &x) { + return fold([](bool r, const T &a) { return r || a; }, false, x); + } + /** Return maximum element + */ + friend constexpr T max_element(const Point &x) { + using std::max; + const T neutral = std::is_floating_point_v + ? -std::numeric_limits::infinity() + : std::numeric_limits::lowest(); + return fold([](const T &r, const T &a) { return max(r, a); }, neutral, x); + } + /** Return minimum element + */ + friend constexpr T min_element(const Point &x) { + using std::min; + const T neutral = std::is_floating_point_v + ? std::numeric_limits::infinity() + : std::numeric_limits::max(); + return fold([](const T &r, const T &a) { return min(r, a); }, neutral, x); + } + /** Product of all elements + */ + friend constexpr T product(const Point &x) { + return fold([](const T &r, const T &a) { return r * a; }, T(1), x); + } + /** Sum of all elements + */ + friend constexpr T sum(const Point &x) { + return fold([](const T &r, const T &a) { return r + a; }, T(0), x); + } + + /** Pointwise comparison + */ + friend constexpr Point operator==(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a == b; }, x, y); + } + /** Pointwise comparison + */ + friend constexpr Point operator!=(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a != b; }, x, y); + } + /** Pointwise comparison + */ + friend constexpr Point operator<(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a < b; }, x, y); + } + /** Pointwise comparison + */ + friend constexpr Point operator>(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a > b; }, x, y); + } + /** Pointwise comparison + */ + friend constexpr Point operator<=(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a <= b; }, x, y); + } + /** Pointwise comparison + */ + friend constexpr Point operator>=(const Point &x, const Point &y) { + return fmap([](const T &a, const T &b) { return a >= b; }, x, y); + } + + friend constexpr Point operator==(const Point &x, const T &b) { + return fmap([&](const T &a) { return a == b; }, x); + } + friend constexpr Point operator!=(const Point &x, const T &b) { + return fmap([&](const T &a) { return a != b; }, x); + } + friend constexpr Point operator<(const Point &x, const T &b) { + return fmap([&](const T &a) { return a < b; }, x); + } + friend constexpr Point operator>(const Point &x, const T &b) { + return fmap([&](const T &a) { return a > b; }, x); + } + friend constexpr Point operator<=(const Point &x, const T &b) { + return fmap([&](const T &a) { return a <= b; }, x); + } + friend constexpr Point operator>=(const Point &x, const T &b) { + return fmap([&](const T &a) { return a >= b; }, x); + } + + /** Output a Point + */ + friend std::ostream &operator<<(std::ostream &os, const Point &x) { + os << "["; + for (std::size_t d = 0; d < D; ++d) { + if (d != 0) + os << ","; + os << x[d]; + } + os << "]"; + return os; + } +}; + +} // namespace Regions +} // namespace openPMD + +namespace std { + +/** Specialization of `equal_to` for Point + */ +template +struct equal_to> { + constexpr bool operator()(const openPMD::Regions::Point &x, + const openPMD::Regions::Point &y) const { + const equal_to> eq; + return eq(x.elts, y.elts); + } +}; + +/** Specialization of `hash` for Point + */ +template +struct hash> { + constexpr size_t operator()(const openPMD::Regions::Point &x) const { + const hash h; + return fold( + [&](size_t r, const T &b) { + return openPMD::Regions::helpers::hash_combine(r, h(b)); + }, + size_t(0xb22da17173243869ULL), x); + } +}; + +/** Specialization of `less` for Point + */ +template +struct less> { + constexpr bool operator()(const openPMD::Regions::Point &x, + const openPMD::Regions::Point &y) const { + const less> lt; + return lt(x.elts, y.elts); + } +}; + +} // namespace std diff --git a/include/openPMD/regions/Region.hpp b/include/openPMD/regions/Region.hpp new file mode 100644 index 0000000000..a7ef816308 --- /dev/null +++ b/include/openPMD/regions/Region.hpp @@ -0,0 +1,1297 @@ +#pragma once + +#include "Box.hpp" +#include "Helpers.hpp" +#include "Point.hpp" + +#include +#include +#include +#include +#include +#include +#include + +namespace openPMD { +namespace Regions { + +/** A D-dimensional region + * + * A region is an arbitrarily shaped set of points. The internal + * representation is likely based on boxes, and is thus most efficient + * if the region has many axis-aligned boundaries. + * + * The dimension D needs to be known at compile time. @see NDRegion + */ +template class Region; + +//////////////////////////////////////////////////////////////////////////////// + +template class Region { +public: + constexpr static std::size_t D = 0; + +private: + bool is_full; + + friend class std::equal_to>; + friend class std::less>; + friend class Region; + + // This should be private, but GCC 9 seems to ignore the friend + // declaration above +public: + explicit constexpr Region(bool is_full_) : is_full(is_full_) {} + +private: +public: + using value_type = typename Point::value_type; + using size_type = typename Point::size_type; + + /** Invariant + */ + constexpr bool invariant() const { return true; } + + void check_invariant() const { +#if REGIONS_DEBUG + assert(invariant()); +#endif + } + + /** Create empty region + */ + constexpr Region() : is_full(false) {} + + Region(const Region &) = default; + Region(Region &&) = default; + Region &operator=(const Region &) = default; + Region &operator=(Region &&) = default; + + Region(const Point &) : is_full(true) {} + Region(const Box &b) : is_full(!b.empty()) {} + + template + Region(const Region ®ion) : is_full(region.is_full) {} + + Region(const std::vector> &bs) { + is_full = false; + for (const auto &b : bs) + is_full |= !b.empty(); + } + operator std::vector>() const { + if (empty()) + return std::vector>{}; + return std::vector>{Box(Point())}; + } + + // Predicates + size_type ndims() const { return D; } + bool empty() const { return !is_full; } + size_type size() const { return is_full; } + size_type nboxes() const { return is_full; } + + // Comparison operators + friend bool operator==(const Region ®ion1, const Region ®ion2) { + return region1.empty() == region2.empty(); + } + friend bool operator!=(const Region ®ion1, const Region ®ion2) { + return !(region1 == region2); + } + friend bool operator==(const Region ®ion1, const Box &box2) { + return region1 == Region(box2); + } + friend bool operator!=(const Region ®ion1, const Box &box2) { + return region1 != Region(box2); + } + friend bool operator==(const Box &box1, const Region ®ion2) { + return Region(box1) == region2; + } + friend bool operator!=(const Box &box1, const Region ®ion2) { + return Region(box1) != region2; + } + + // Shift and scale operators + Region operator>>(const Point &) const { return *this; } + Region operator<<(const Point &) const { return *this; } + Region &operator>>=(const Point &d) { return *this = *this >> d; } + Region &operator<<=(const Point &d) { return *this = *this << d; } + Region operator*(const Point &) const { return *this; } + Region &operator*=(const Point &s) { return *this = *this * s; } + Region grown(const Point &, const Point &) const { return *this; } + Region grown(const Point &d) const { return grown(d, d); } + Region grown(T n) const { return grown(Point::pure(n)); } + Region shrunk(const Point &, const Point &) const { + return *this; + } + Region shrunk(const Point &d) const { return shrunk(d, d); } + Region shrunk(T n) const { return shrunk(Point::pure(n)); } + + // Set operations + friend Box bounding_box(const Region ®ion) { + if (region.empty()) + return Box(); + return Box(Point()); + } + + friend Region operator&(const Region ®ion1, const Region ®ion2) { + return Region(!region1.empty() & !region2.empty()); + } + friend Region operator|(const Region ®ion1, const Region ®ion2) { + return Region(!region1.empty() | !region2.empty()); + } + friend Region operator^(const Region ®ion1, const Region ®ion2) { + return Region(!region1.empty() ^ !region2.empty()); + } + friend Region operator-(const Region ®ion1, const Region ®ion2) { + return Region(!region1.empty() & region2.empty()); + } + + Region &operator&=(const Region ®ion) { return *this = *this & region; } + Region &operator|=(const Region ®ion) { return *this = *this | region; } + Region &operator^=(const Region ®ion) { return *this = *this ^ region; } + Region &operator-=(const Region ®ion) { return *this = *this - region; } + + friend Region intersection(const Region ®ion1, const Region ®ion2) { + return region1 & region2; + } + friend Region setunion(const Region ®ion1, const Region ®ion2) { + return region1 | region2; + } + friend Region symmetric_difference(const Region ®ion1, + const Region ®ion2) { + return region1 ^ region2; + } + friend Region difference(const Region ®ion1, const Region ®ion2) { + return region1 - region2; + } + + // Set comparison operators + bool contains(const Point &p) const { + return !isdisjoint(*this, Region(p)); + } + friend bool isdisjoint(const Region ®ion1, const Region ®ion2) { + return (region1 & region2).empty(); + } + + // Comparison operators + friend bool operator<=(const Region ®ion1, const Region ®ion2) { + return (region1 - region2).empty(); + } + friend bool operator>=(const Region ®ion1, const Region ®ion2) { + return region2 <= region1; + } + friend bool operator<(const Region ®ion1, const Region ®ion2) { + return region1 != region2 && region1 <= region2; + } + friend bool operator>(const Region ®ion1, const Region ®ion2) { + return region2 < region1; + } + + bool is_subset_of(const Region ®ion) const { return *this <= region; } + bool is_superset_of(const Region ®ion) const { return *this >= region; } + bool is_strict_subset_of(const Region ®ion) const { + return *this < region; + } + bool is_strict_superset_of(const Region ®ion) const { + return *this > region; + } + + friend std::ostream &operator<<(std::ostream &os, const Region ®ion) { + os << "{"; + if (!region.empty()) + os << "(1)"; + os << "}"; + return os; + } +}; + +//////////////////////////////////////////////////////////////////////////////// + +template class Region { +public: + constexpr static std::size_t D = 1; + +private: + using Subregion = Region; // This is essentially a bool + using Subregions = std::vector; + Subregions subregions; + + friend class std::equal_to>; + friend class std::less>; + +public: + using value_type = typename Point::value_type; + using size_type = typename Point::size_type; + + /** Invariant + */ + bool invariant() const { + const size_type nboxes = subregions.size(); + for (size_type i = 1; i < nboxes; ++i) + // The subregions must be ordered + if (!(subregions[i - 1] < subregions[i])) { + assert(false); + return false; + } + // There must be an even number of subregions + if (subregions.size() % 2 != 0) { + assert(false); + return false; + } + return true; + } + + void check_invariant() const { +#if REGIONS_DEBUG + assert(invariant()); +#endif + } + + /** Create empty region + */ + Region() = default; + + Region(const Region &) = default; + Region(Region &&) = default; + Region &operator=(const Region &) = default; + Region &operator=(Region &&) = default; + + Region(const Point &p) : Region(Box(p)) {} + Region(const Box &b) { + if (b.empty()) + return; + subregions = {b.lower()[0], b.upper()[0]}; + check_invariant(); + } + + template + Region(const Region ®ion) : subregions(region.subregions.size()) { + std::transform(region.subregions.begin(), region.subregions.end(), + subregions.begin(), + [](const auto &subregion) { return Subregion(subregion); }); + } + +private: + static std::vector subregions_from_bounds(std::vector lbnds, + std::vector ubnds) { + const std::size_t nboxes = lbnds.size(); + assert(ubnds.size() == nboxes); + std::vector subregions; + if (nboxes == 0) + return subregions; + std::sort(lbnds.begin(), lbnds.end()); + std::sort(ubnds.begin(), ubnds.end()); + std::size_t nactive = 0; + std::size_t lpos = 0, upos = 0; + while (lpos < nboxes) { + const auto lbnd = lbnds[lpos]; + assert(upos < nboxes); + const auto ubnd = ubnds[upos]; + // Process lower bounds before upper bounds + if (lbnd <= ubnd) { + if (nactive == 0) + subregions.push_back(lbnd); + ++nactive; + ++lpos; + } else { + assert(nactive > 0); + --nactive; + if (nactive == 0) + subregions.push_back(ubnd); + ++upos; + } + } + assert(nactive > 0); + assert(upos < nboxes); + assert(upos + nactive == nboxes); + subregions.push_back(ubnds[nboxes - 1]); +#if REGIONS_DEBUG + // TODO: turn this into a test case + { + Region reg; + for (std::size_t i = 0; i < nboxes; ++i) + reg |= Region(Box(Point{lbnds[i]}, Point{ubnds[i]})); + assert(subregions == reg.subregions); + } +#endif + return subregions; + } + +public: + Region(const std::vector> &boxes) { + std::vector lbnds, ubnds; + lbnds.reserve(boxes.size()); + ubnds.reserve(boxes.size()); + for (const auto &box : boxes) { + lbnds.push_back(box.lower()[0]); + ubnds.push_back(box.upper()[0]); + } + subregions = subregions_from_bounds(std::move(lbnds), std::move(ubnds)); + check_invariant(); + } + operator std::vector>() const { + std::vector> res; + res.reserve(subregions.size() / 2); + auto iter = subregions.begin(); + const auto end = subregions.end(); + while (iter != end) { + const auto pos1 = *iter++; + const auto pos2 = *iter++; + res.emplace_back(Box(Point{pos1}, Point{pos2})); + } +#if REGIONS_DEBUG + // TODO: turn this into a test case + assert(std::is_sorted(res.begin(), res.end())); + { + Region reg; + for (const auto &b : res) { + assert(isdisjoint(Region(b), reg)); + reg |= b; + } + if (!(reg == *this)) { + std::cout << "subregions=["; + for (const auto &sr : subregions) + std::cout << sr << ","; + std::cout << "]\n"; + std::cout << "reg=["; + for (const auto &r : reg.subregions) + std::cout << r << ","; + std::cout << "]\n"; + } + assert(reg == *this); + } +#endif + return res; + } + +private: + template + friend void traverse_subregions(const F &f, const Region ®ion) { + Subregion decoded_subregion; + for (const auto &pos : region.subregions) { + decoded_subregion ^= Subregion(Point()); + f(pos, decoded_subregion); + } + assert(decoded_subregion.empty()); + } + + template + friend void traverse_subregions(const F &f, const Region ®ion1, + const Region ®ion2) { + Subregion decoded_subregion1, decoded_subregion2; + + auto iter1 = region1.subregions.begin(); + auto iter2 = region2.subregions.begin(); + const auto end1 = region1.subregions.end(); + const auto end2 = region2.subregions.end(); + while (iter1 != end1 && iter2 != end2) { + const T next_pos1 = *iter1; + const T next_pos2 = *iter2; + using std::min; + const T pos = min(next_pos1, next_pos2); + const bool active1 = next_pos1 == pos; + const bool active2 = next_pos2 == pos; + decoded_subregion1 ^= Subregion(active1); + decoded_subregion2 ^= Subregion(active2); + f(pos, decoded_subregion1, decoded_subregion2); + if (active1) + ++iter1; + if (active2) + ++iter2; + } + for (; iter1 != end1; ++iter1) { + const T pos = *iter1; + decoded_subregion1 ^= Subregion(Point()); + f(pos, decoded_subregion1, Subregion()); + } + for (; iter2 != end2; ++iter2) { + const T pos = *iter2; + decoded_subregion2 ^= Subregion(Point()); + f(pos, Subregion(), decoded_subregion2); + } + assert(decoded_subregion1.empty()); + assert(decoded_subregion2.empty()); + } + + template + friend Region unary_operator(const F &op, const Region ®ion1) { + Region res; + Subregion old_decoded_subregion; + traverse_subregions( + [&](const T pos, const Subregion &decoded_subregion1) { + auto decoded_subregion_res = op(decoded_subregion1); + auto subregion = decoded_subregion_res ^ old_decoded_subregion; + if (!subregion.empty()) + res.subregions.push_back(pos); + old_decoded_subregion = decoded_subregion_res; + }, + region1); + assert(old_decoded_subregion.empty()); + assert(res.invariant()); + return res; + } + + template + friend Region binary_operator(const F &op, const Region ®ion1, + const Region ®ion2) { + Region res; + Subregion old_decoded_subregion; + traverse_subregions( + [&](const T pos, const Subregion &decoded_subregion1, + const Subregion &decoded_subregion2) { + auto decoded_subregion_res = + op(decoded_subregion1, decoded_subregion2); + auto subregion = decoded_subregion_res ^ old_decoded_subregion; + if (!subregion.empty()) + res.subregions.push_back(pos); + old_decoded_subregion = std::move(decoded_subregion_res); + }, + region1, region2); + assert(old_decoded_subregion.empty()); + assert(res.invariant()); + return res; + } + +public: + // Predicates + size_type ndims() const { return D; } + bool empty() const { return subregions.empty(); } + + size_type size() const { + size_type total_size = 0; + auto iter = subregions.begin(); + const auto end = subregions.end(); + while (iter != end) { + const auto pos0 = *iter++; + const auto pos1 = *iter++; + total_size += pos1 - pos0; + } + return total_size; + } + + size_type nboxes() const { return subregions.size(); } + + // Comparison operators + friend bool operator==(const Region ®ion1, const Region ®ion2) { + return region1.subregions == region2.subregions; + } + friend bool operator!=(const Region ®ion1, const Region ®ion2) { + return !(region1 == region2); + } + friend bool operator==(const Region ®ion1, const Box &box2) { + return region1 == Region(box2); + } + friend bool operator!=(const Region ®ion1, const Box &box2) { + return region1 != Region(box2); + } + friend bool operator==(const Box &box1, const Region ®ion2) { + return Region(box1) == region2; + } + friend bool operator!=(const Box &box1, const Region ®ion2) { + return Region(box1) != region2; + } + + // Shift and scale operators + Region operator>>(const Point &d) const { + Region nr; + nr.subregions.reserve(subregions.size()); + for (const auto &pos : subregions) + nr.subregions.push_back(pos + d[0]); + check_invariant(); + return nr; + } + Region operator<<(const Point &d) const { return *this >> -d; } + Region &operator>>=(const Point &d) { return *this = *this >> d; } + Region &operator<<=(const Point &d) { return *this = *this << d; } + Region operator*(const Point &s) const { + if (s[0] == 0) + return empty() ? Region() : Region(Point()); + Region nr; + nr.subregions.reserve(subregions.size()); + for (const auto &pos : subregions) + nr.subregions.push_back(pos * s[0]); + if (s[0] < 0) + std::reverse(nr.subregions.begin(), nr.subregions.end()); + check_invariant(); + return nr; + } + Region &operator*=(const Point &s) { return *this = *this * s; } + +private: + Region grown_(const Point &dlo, const Point &dup) const { + // Cannot shrink + assert(all(dlo + dup >= Point())); + return helpers::mapreduce( + [&](const Box &b) { return Region(b.grown(dlo, dup)); }, + [](const Region &x, const Region &y) { return x | y; }, + std::vector>(*this)); + } + Region shrunk_(const Point &dlo, const Point &dup) const { + // Cannot grow + assert(all(dlo + dup >= Point())); + auto world = bounding_box(*this).grown(1); + return Region(world.grown(dup, dlo)) - + (Region(world) - *this).grown(dup, dlo); + } + +public: + Region grown(const Point &dlo, const Point &dup) const { + const Region ®ion0 = *this; + const bool need_shrink = any(dlo + dup < Point()); + Region shrunk_region; + if (need_shrink) { + // Shrink in some directions + const Point ndlo = fmap( + [](auto lo, auto up) { return lo + up < 0 ? lo : T(0); }, dlo, dup); + const Point ndup = fmap( + [](auto lo, auto up) { return lo + up < 0 ? up : T(0); }, dlo, dup); + shrunk_region = region0.shrunk_(-ndlo, -ndup); + } + const Region ®ion1 = need_shrink ? shrunk_region : *this; + const bool need_grow = any(dlo + dup > Point()); + Region grown_region; + if (need_grow) { + // Grow in some direction + const Point ndlo = fmap( + [](auto lo, auto up) { return lo + up > 0 ? lo : T(0); }, dlo, dup); + const Point ndup = fmap( + [](auto lo, auto up) { return lo + up > 0 ? up : T(0); }, dlo, dup); + grown_region = region1.grown_(ndlo, ndup); + } + const Region ®ion2 = need_grow ? grown_region : region1; + return region2; + } + Region grown(const Point &d) const { return grown(d, d); } + Region grown(const T &n) const { return grown(Point::pure(n)); } + Region shrunk(const Point &dlo, const Point &dup) const { + return grown(-dlo, -dup); + } + Region shrunk(const Point &d) const { return shrunk(d, d); } + Region shrunk(T n) const { return shrunk(Point::pure(n)); } + + // Set operations + friend Box bounding_box(const Region ®ion) { + if (region.empty()) + return Box(); + return Box(Point{*region.subregions.begin()}, + Point{*(region.subregions.end() - 1)}); + } + + friend Region operator&(const Region ®ion1, const Region ®ion2) { + return binary_operator([](const Subregion &set1, + const Subregion &set2) { return set1 & set2; }, + region1, region2); + } + friend Region operator|(const Region ®ion1, const Region ®ion2) { + return binary_operator([](const Subregion &set1, + const Subregion &set2) { return set1 | set2; }, + region1, region2); + } + friend Region operator^(const Region ®ion1, const Region ®ion2) { + return binary_operator([](const Subregion &set1, + const Subregion &set2) { return set1 ^ set2; }, + region1, region2); + } + friend Region operator-(const Region ®ion1, const Region ®ion2) { + return binary_operator([](const Subregion &set1, + const Subregion &set2) { return set1 - set2; }, + region1, region2); + } + + Region &operator&=(const Region ®ion) { return *this = *this & region; } + Region &operator|=(const Region ®ion) { return *this = *this | region; } + Region &operator^=(const Region ®ion) { return *this = *this ^ region; } + Region &operator-=(const Region ®ion) { return *this = *this - region; } + + friend Region intersection(const Region ®ion1, const Region ®ion2) { + return region1 & region2; + } + friend Region setunion(const Region ®ion1, const Region ®ion2) { + return region1 | region2; + } + friend Region symmetric_difference(const Region ®ion1, + const Region ®ion2) { + return region1 ^ region2; + } + friend Region difference(const Region ®ion1, const Region ®ion2) { + return region1 - region2; + } + + // Set comparison operators + bool contains(const Point &p) const { + return !isdisjoint(*this, Region(p)); + } + friend bool isdisjoint(const Region ®ion1, const Region ®ion2) { + return (region1 & region2).empty(); + } + + // Comparison operators + friend bool operator<=(const Region ®ion1, const Region ®ion2) { + return (region1 - region2).empty(); + } + friend bool operator>=(const Region ®ion1, const Region ®ion2) { + return region2 <= region1; + } + friend bool operator<(const Region ®ion1, const Region ®ion2) { + return region1 != region2 && region1 <= region2; + } + friend bool operator>(const Region ®ion1, const Region ®ion2) { + return region2 < region1; + } + + bool is_subset_of(const Region ®ion) const { return *this <= region; } + bool is_superset_of(const Region ®ion) const { return *this >= region; } + bool is_strict_subset_of(const Region ®ion) const { + return *this < region; + } + bool is_strict_superset_of(const Region ®ion) const { + return *this > region; + } + + friend std::ostream &operator<<(std::ostream &os, const Region ®ion) { + os << "{"; + const std::vector> boxes(region); + for (std::size_t i = 0; i < boxes.size(); ++i) { + if (i > 0) + os << ","; + os << boxes[i]; + } + os << "}"; + return os; + } +}; + +//////////////////////////////////////////////////////////////////////////////// + +template class Region { + using Subregion = Region; + using Subregions = std::vector>; + Subregions subregions; + + friend class std::equal_to>; + friend class std::less>; + +public: + using value_type = typename Point::value_type; + using size_type = typename Point::size_type; + + /** Invariant + */ + bool invariant() const { + for (const auto &pos_subregion : subregions) { + const auto &subregion = pos_subregion.second; + // The subregions must not be empty, and their invariants must hold + if (subregion.empty() || !subregion.invariant()) { + assert(false); + return false; + } + } + return true; + } + + /** Check invariant when in debug mode; abort the code if invariant is not + * met. + */ + void check_invariant() const { +#if REGIONS_DEBUG + assert(invariant()); +#endif + } + + /** Create empty Region + */ + Region() = default; + + Region(const Region &) = default; + Region(Region &&) = default; + Region &operator=(const Region &) = default; + Region &operator=(Region &&) = default; + + /** Create Region containing a single Point + */ + Region(const Point &p) : Region(Box(p)) {} + /** Create Region containina a single Box + */ + Region(const Box &b) { + if (b.empty()) + return; + Box subbox(b.lower().erase(D - 1), b.upper().erase(D - 1)); + subregions = {{b.lower()[D - 1], Subregion(subbox)}, + {b.upper()[D - 1], Subregion(subbox)}}; + check_invariant(); + } + + template + Region(const Region ®ion) : subregions(region.subregions.size()) { + std::transform(region.subregions.begin(), region.subregions.end(), + subregions.begin(), [&](const auto &pos_subregion) { + T pos(pos_subregion.first); + Subregion subregion(pos_subregion.second); + subregions.emplace_back( + std::make_pair(pos, std::move(subregion))); + }); + } + +private: + static Region region_from_boxes( + const typename std::vector>::const_iterator &begin, + const typename std::vector>::const_iterator &end) { + auto sz = end - begin; + if (sz == 0) + return Region(); + if (sz == 1) + return Region(*begin); + const auto mid = begin + sz / 2; + return region_from_boxes(begin, mid) | region_from_boxes(mid, end); + } + +public: + /** Create a Region from a vector of Boxes + */ + Region(const std::vector> &boxes) { + *this = region_from_boxes(boxes.begin(), boxes.end()); +#if REGIONS_DEBUG + // TODO: turn this into a test case + { + Region reg; + for (const auto &box : boxes) + reg |= Region(box); + assert(*this == reg); + } +#endif + check_invariant(); + } + + /** Convert a Region to a vector of non-overlapping Boxes + * + * This creates a canonical, reasonably efficient representation. + */ + operator std::vector>() const { + std::vector> res; + std::map, T> old_subboxes; + traverse_subregions( + [&](const T pos, const Subregion &subregion) { + // Convert subregion to boxes + const std::vector> subboxes1(subregion); + + auto iter0 = old_subboxes.begin(); + auto iter1 = subboxes1.begin(); + const auto end0 = old_subboxes.end(); + const auto end1 = subboxes1.end(); +#if REGIONS_DEBUG + assert(is_sorted(iter1, end1)); +#endif + std::map, T> subboxes; + while (iter0 != end0 || iter1 != end1) { + bool active0 = iter0 != end0; + bool active1 = iter1 != end1; + Box dummy; + const Box &subbox0 = active0 ? iter0->first : dummy; + const Box &subbox1 = active1 ? *iter1 : dummy; + // When both subboxes are active, keep only the first (as determined + // by less<>) + std::equal_to eq; + std::less lt; + if (active0 && active1) { + active0 = !lt(subbox0, subbox1); + active1 = !lt(subbox1, subbox0); + } + + if (active0 && active1 && eq(subbox0, subbox1)) { + // The current bbox continues unchanged -- keep it + const T old_pos = iter0->second; + subboxes[subbox1] = old_pos; + } else { + if (active0) { + // The current box changed; finalize it + const T old_pos = iter0->second; + res.push_back(Box(subbox0.lower().insert(D - 1, old_pos), + subbox0.upper().insert(D - 1, pos))); + } + if (active1) + // There is a new box; add it + subboxes[subbox1] = pos; + } + + if (active0) + ++iter0; + if (active1) + ++iter1; + } + old_subboxes = std::move(subboxes); + }, + *this); + assert(old_subboxes.empty()); +#if REGIONS_DEBUG + // TODO: turn this into a test case + assert(std::is_sorted(res.begin(), res.end())); + { + Region reg; + for (const auto &b : res) { + assert(isdisjoint(Region(b), reg)); + reg |= b; + } + assert(reg == *this); + } +#endif + return res; + } + +private: + template + friend void traverse_subregions(const F &f, const Region ®ion) { + Subregion decoded_subregion; + for (const auto &pos_subregion : region.subregions) { + const T pos = pos_subregion.first; + const auto &subregion = pos_subregion.second; + decoded_subregion ^= subregion; + f(pos, decoded_subregion); + } + assert(decoded_subregion.empty()); + } + + template + friend void traverse_subregions(const F &f, const Region ®ion1, + const Region ®ion2) { + Subregion decoded_subregion1, decoded_subregion2; + + auto iter1 = region1.subregions.begin(); + auto iter2 = region2.subregions.begin(); + const auto end1 = region1.subregions.end(); + const auto end2 = region2.subregions.end(); + while (iter1 != end1 || iter2 != end2) { + const T next_pos1 = + iter1 != end1 ? iter1->first : std::numeric_limits::max(); + const T next_pos2 = + iter2 != end2 ? iter2->first : std::numeric_limits::max(); + using std::min; + const T pos = min(next_pos1, next_pos2); + const bool active1 = next_pos1 == pos; + const bool active2 = next_pos2 == pos; + Subregion dummy; + const Subregion &subregion1 = active1 ? iter1->second : dummy; + const Subregion &subregion2 = active2 ? iter2->second : dummy; + if (active1) + decoded_subregion1 ^= subregion1; + if (active2) + decoded_subregion2 ^= subregion2; + + f(pos, decoded_subregion1, decoded_subregion2); + + if (active1) + ++iter1; + if (active2) + ++iter2; + } + assert(decoded_subregion1.empty()); + assert(decoded_subregion2.empty()); + } + + template + friend Region unary_operator(const F &op, const Region ®ion1) { + Region res; + Subregion old_decoded_subregion; + traverse_subregions( + [&](const T pos, const Subregion &decoded_subregion1) { + auto decoded_subregion_res = op(decoded_subregion1); + auto subregion = decoded_subregion_res ^ old_decoded_subregion; + if (!subregion.empty()) + res.subregions.emplace_back( + std::make_pair(pos, std::move(subregion))); + old_decoded_subregion = std::move(decoded_subregion_res); + }, + region1); + assert(old_decoded_subregion.empty()); + res.check_invariant(); + return res; + } + + template + friend Region binary_operator(const F &op, const Region ®ion1, + const Region ®ion2) { + Region res; + Subregion old_decoded_subregion; + traverse_subregions( + [&](const T pos, const Subregion &decoded_subregion1, + const Subregion &decoded_subregion2) { + auto decoded_subregion_res = + op(decoded_subregion1, decoded_subregion2); + auto subregion = decoded_subregion_res ^ old_decoded_subregion; + if (!subregion.empty()) + res.subregions.emplace_back( + std::make_pair(pos, std::move(subregion))); + old_decoded_subregion = std::move(decoded_subregion_res); + }, + region1, region2); + assert(old_decoded_subregion.empty()); + res.check_invariant(); + return res; + } + +public: + // Predicates + /** Number of dimensions + */ + size_type ndims() const { return D; } + /** Whether the Region is empty + */ + bool empty() const { return subregions.empty(); } + + /** Size, the total number of contained points + */ + size_type size() const { + size_type total_size = 0; + T old_pos = std::numeric_limits::min(); // location of last subregion + size_type old_subregion_size = 0; // number of points in the last subregion + traverse_subregions( + [&](const T pos, const Subregion &subregion) { + const size_type subregion_size = subregion.size(); + total_size += old_subregion_size == 0 + ? 0 + : (pos - old_pos) * old_subregion_size; + old_pos = pos; + old_subregion_size = subregion_size; + }, + *this); + assert(old_subregion_size == 0); + return total_size; + } + + /** A measure of the number of vertices defining the Region + */ + size_type nboxes() const { + size_type sz = 0; + for (const auto &pos_subregion : subregions) { + const auto &subregion = pos_subregion.second; + sz += subregion.nboxes(); + } + return sz; + } + + // Shift and scale operators + /** Shift a Region to the right (upwards). The shift can be negative, which + * shifts left. + */ + Region operator>>(const Point &d) const { + Region nr; + nr.subregions.reserve(subregions.size()); + const T dx = d[D - 1]; + auto subd = d.erase(D - 1); + for (const auto &pos_subregion : subregions) { + const T pos = pos_subregion.first; + const auto &subregion = pos_subregion.second; + nr.subregions.emplace_back(std::make_pair(pos + dx, subregion >> subd)); + } + nr.check_invariant(); + return nr; + } + /** Shift a Region to the left (downwards). The shift can be negative, which + * shifts right. + */ + Region operator<<(const Point &d) const { return *this >> -d; } + Region &operator>>=(const Point &d) { return *this = *this >> d; } + Region &operator<<=(const Point &d) { return *this = *this << d; } + /** Scale a Region + */ + Region operator*(const Point &s) const { + if (s[D - 1] == 0) + return empty() ? Region() : Region(Point()); + Region nr; + nr.subregions.reserve(subregions.size()); + const T ds = s[D - 1]; + auto subs = s.erase(D - 1); + for (const auto &pos_subregion : subregions) { + const T pos = pos_subregion.first; + const auto &subregion = pos_subregion.second; + nr.subregions.emplace_back(std::make_pair(pos * ds, subregion * subs)); + } + if (ds < 0) { + std::reverse(nr.subregions.begin(), nr.subregions.end()); + std::transform(nr.subregions.begin(), nr.subregions.end(), + nr.subregions.begin(), [](const auto &pos_subregion) { + const auto &pos = pos_subregion.first; + const auto &subregion = pos_subregion.second; + return std::make_pair(pos + 1, subregion); + }); + } + nr.check_invariant(); + return nr; + } + Region &operator*=(const Point &s) { return *this = *this * s; } + +private: + Region grown_(const Point &dlo, const Point &dup) const { + // Cannot shrink + assert(all(dlo + dup >= Point())); + return helpers::mapreduce( + [&](const Box &b) { return Region(b.grown(dlo, dup)); }, + [](const Region &x, const Region &y) { return x | y; }, + std::vector>(*this)); + } + Region shrunk_(const Point &dlo, const Point &dup) const { + // Cannot grow + assert(all(dlo + dup >= Point())); + auto world = bounding_box(*this).grown(1); + return Region(world.grown(dup, dlo)) - + (Region(world) - *this).grown(dup, dlo); + } + +public: + /** Grow a Region by given amounts in each direction. + * + * The growth can be negative, which shrinks the Region. If a Region is + * shrunk too much it becomes empty. Growing an empty Region always + * results in an empty Region. + */ + Region grown(const Point &dlo, const Point &dup) const { + const Region ®ion0 = *this; + const bool need_shrink = any(dlo + dup < Point()); + Region shrunk_region; + if (need_shrink) { + // Shrink in some directions + const Point ndlo = fmap( + [](auto lo, auto up) { return lo + up < 0 ? lo : T(0); }, dlo, dup); + const Point ndup = fmap( + [](auto lo, auto up) { return lo + up < 0 ? up : T(0); }, dlo, dup); + shrunk_region = region0.shrunk_(-ndlo, -ndup); + } + const Region ®ion1 = need_shrink ? shrunk_region : *this; + const bool need_grow = any(dlo + dup > Point()); + Region grown_region; + if (need_grow) { + // Grow in some direction + const Point ndlo = fmap( + [](auto lo, auto up) { return lo + up > 0 ? lo : T(0); }, dlo, dup); + const Point ndup = fmap( + [](auto lo, auto up) { return lo + up > 0 ? up : T(0); }, dlo, dup); + grown_region = region1.grown_(ndlo, ndup); + } + const Region ®ion2 = need_grow ? grown_region : region1; + return region2; + } + Region grown(const Point &d) const { return grown(d, d); } + Region grown(const T &n) const { return grown(Point::pure(n)); } + /** Shrink a Region by given amounts in each direction. + * + * The shrinkage can be negative, which shrinks the Region. If a Region is + * shrunk too much it becomes empty. Growing an empty Region always + * results in an empty Region. + */ + Region shrunk(const Point &dlo, const Point &dup) const { + return grown(-dlo, -dup); + } + Region shrunk(const Point &d) const { return shrunk(d, d); } + Region shrunk(const T &n) const { return shrunk(Point::pure(n)); } + + // Comparison operators + /** Compare two Regions + */ + friend bool operator==(const Region ®ion1, const Region ®ion2) { + return region1.subregions == region2.subregions; + } + /** Compare two Regions + */ + friend bool operator!=(const Region ®ion1, const Region ®ion2) { + return !(region1 == region2); + } + friend bool operator==(const Region ®ion1, const Box &box2) { + return region1 == Region(box2); + } + friend bool operator!=(const Region ®ion1, const Box &box2) { + return region1 != Region(box2); + } + friend bool operator==(const Box &box1, const Region ®ion2) { + return Region(box1) == region2; + } + friend bool operator!=(const Box &box1, const Region ®ion2) { + return Region(box1) != region2; + } + + // Set operations + /** Calculate the bounding box of a REgion. This is the smallest Box that + * contains the Region. + */ + friend Box bounding_box(const Region ®ion) { + if (region.empty()) + return Box(); + auto pmin = Point::pure(std::numeric_limits::max()); + auto pmax = Point::pure(std::numeric_limits::min()); + for (const auto &pos_subregion : region.subregions) { + const auto &subregion = pos_subregion.second; + auto subbox = bounding_box(subregion); + using std::max, std::min; + pmin = min(pmin, subbox.lower()); + pmax = max(pmax, subbox.upper()); + } + const T xmin = region.subregions.begin()->first; + const T xmax = (region.subregions.end() - 1)->first; + return Box(pmin.insert(D - 1, xmin), pmax.insert(D - 1, xmax)); + } + + /** Set intersection of two Regions + */ + friend Region operator&(const Region ®ion1, const Region ®ion2) { + return binary_operator([](const Subregion &set1, + const Subregion &set2) { return set1 & set2; }, + region1, region2); + } + /** Set union of two Regions + */ + friend Region operator|(const Region ®ion1, const Region ®ion2) { + return binary_operator([](const Subregion &set1, + const Subregion &set2) { return set1 | set2; }, + region1, region2); + } + /** Symmetric difference of two Regions + */ + friend Region operator^(const Region ®ion1, const Region ®ion2) { + // TODO: If region2 is much smaller than region1, direct insertion may be + // faster + return binary_operator([](const Subregion &set1, + const Subregion &set2) { return set1 ^ set2; }, + region1, region2); + } + /** Set difference of two Regions + */ + friend Region operator-(const Region ®ion1, const Region ®ion2) { + return binary_operator([](const Subregion &set1, + const Subregion &set2) { return set1 - set2; }, + region1, region2); + } + + Region &operator&=(const Region ®ion) { return *this = *this & region; } + Region &operator|=(const Region ®ion) { return *this = *this | region; } + Region &operator^=(const Region ®ion) { return *this = *this ^ region; } + Region &operator-=(const Region ®ion) { return *this = *this - region; } + + /** Set intersection of two Regions + */ + friend Region intersection(const Region ®ion1, const Region ®ion2) { + return region1 & region2; + } + /** Set union of two Regions + */ + friend Region setunion(const Region ®ion1, const Region ®ion2) { + return region1 | region2; + } + /** Symmetric difference of two Regions + */ + friend Region symmetric_difference(const Region ®ion1, + const Region ®ion2) { + return region1 ^ region2; + } + /** Set difference of two Regions + */ + friend Region difference(const Region ®ion1, const Region ®ion2) { + return region1 - region2; + } + + // Set comparison operators + /** Whether a Region contains a Point + */ + bool contains(const Point &p) const { + return !isdisjoint(*this, Region(p)); + } + /** Whether two Regions are disjoint, i.e. whether they have no Point in + * common + */ + friend bool isdisjoint(const Region ®ion1, const Region ®ion2) { + return (region1 & region2).empty(); + } + + // Comparison operators + /** is subset of + */ + friend bool operator<=(const Region ®ion1, const Region ®ion2) { + return (region1 - region2).empty(); + } + /** is superset of + */ + friend bool operator>=(const Region ®ion1, const Region ®ion2) { + return region2 <= region1; + } + /** is strict subset of + */ + friend bool operator<(const Region ®ion1, const Region ®ion2) { + return region1 != region2 && region1 <= region2; + } + /** is strict superset of + */ + friend bool operator>(const Region ®ion1, const Region ®ion2) { + return region2 < region1; + } + + /** is subset of + */ + bool is_subset_of(const Region ®ion) const { return *this <= region; } + /** is superset of + */ + bool is_superset_of(const Region ®ion) const { return *this >= region; } + /** is strict subset of + */ + bool is_strict_subset_of(const Region ®ion) const { + return *this < region; + } + /** is strict superset of + */ + bool is_strict_superset_of(const Region ®ion) const { + return *this > region; + } + + /** Output a Region + */ + friend std::ostream &operator<<(std::ostream &os, const Region ®ion) { + os << "{"; + const std::vector> boxes(region); + for (std::size_t i = 0; i < boxes.size(); ++i) { + if (i > 0) + os << ","; + os << boxes[i]; + } + os << "}"; + return os; + } +}; + +} // namespace Regions +} // namespace openPMD + +namespace std { + +template +struct equal_to>; + +template struct equal_to> { + constexpr bool operator()(const openPMD::Regions::Region &x, + const openPMD::Regions::Region &y) const { + return x.is_full == y.is_full; + } +}; + +template +struct equal_to> { + constexpr bool operator()(const openPMD::Regions::Region &x, + const openPMD::Regions::Region &y) const { + return openPMD::Regions::helpers::vector_eq(x.subregions, y.subregions); + } +}; + +template +struct less>; + +template struct less> { + constexpr bool operator()(const openPMD::Regions::Region &x, + const openPMD::Regions::Region &y) const { + return x.is_full < y.is_full; + } +}; + +template +struct less> { + constexpr bool operator()(const openPMD::Regions::Region &x, + const openPMD::Regions::Region &y) const { + return openPMD::Regions::helpers::vector_lt(x.subregions, y.subregions); + } +}; + +} // namespace std diff --git a/include/openPMD/regions/Regions.hpp b/include/openPMD/regions/Regions.hpp new file mode 100644 index 0000000000..fe4c9bfb2c --- /dev/null +++ b/include/openPMD/regions/Regions.hpp @@ -0,0 +1,9 @@ +#pragma once + +#include "Box.hpp" +#include "Point.hpp" +#include "Region.hpp" + +#include "NDBox.hpp" +#include "NDPoint.hpp" +#include "NDRegion.hpp" diff --git a/test/BoxTest.cpp b/test/BoxTest.cpp new file mode 100644 index 0000000000..c739a12fa7 --- /dev/null +++ b/test/BoxTest.cpp @@ -0,0 +1,226 @@ +#include +#include + +#include + +#include +#include +#include +#include +#include + +using namespace openPMD::Regions; + +template void test_Box(const B &box) { + const std::size_t D = box.ndims(); + using T = typename B::value_type; + const auto p = box.lower(); + typedef std::decay_t P; + const std::equal_to

eqp; + const std::equal_to eqb; + const std::less ltb; + CHECK(box.empty()); + + std::mt19937 gen; + std::uniform_int_distribution dist0(0, 9); + std::uniform_int_distribution dist(-1000, 1000); + const auto rand = [&]() { return dist(gen); }; + const auto randp = [&]() { return fmap([&](auto) { return rand(); }, p); }; + const auto randb = [&]() { + if (D == 0) { + if (dist0(gen) < 5) + return B(p, p); // empty box + else + return B(p); + } + if (dist0(gen) == 0) + return box; + while (1) { + auto lo = randp(); + auto hi = randp(); + auto nb = B(lo, hi); + if (!nb.empty()) + return nb; + } + }; + + for (int iter = 0; iter < 100; ++iter) { + + B N(box); + CHECK(N.ndims() == D); + CHECK(N.empty()); + for (std::size_t d = 0; d < D; ++d) + CHECK(N.lower()[d] >= N.upper()[d]); + + const B X = randb(); + const B Y = randb(); + const B Z = randb(); + + const P n = fmap([](auto) { return T(0); }, p); + CHECK(eqp(n + n, n)); + const P x = randp(); + const P y = randp(); + + const T a = rand(); + + CHECK(N.empty()); + if (D > 0) { + CHECK(X.empty() == all(X.shape() == 0)); + CHECK(Y.empty() == all(Y.shape() == 0)); + CHECK(Z.empty() == all(Z.shape() == 0)); + } + + CHECK(X.empty() == (X.size() == 0)); + CHECK(Y.empty() == (Y.size() == 0)); + CHECK(Z.empty() == (Z.size() == 0)); + + if (D > 0) { + CHECK(X.empty() == all(fmap([](auto lo, auto up) { return lo >= up; }, + X.lower(), X.upper()))); + CHECK(Y.empty() == all(fmap([](auto lo, auto up) { return lo >= up; }, + Y.lower(), Y.upper()))); + CHECK(Z.empty() == all(fmap([](auto lo, auto up) { return lo >= up; }, + Z.lower(), Z.upper()))); + } + + CHECK(eqb(N, N)); + CHECK(eqb(X, X)); + CHECK(!ltb(N, N)); + CHECK(!ltb(X, X)); + if (X.empty()) { + CHECK(eqb(N, X)); + CHECK(!ltb(N, X)); + } else { + CHECK(!eqb(N, X)); + CHECK(ltb(N, X)); + } + if (eqb(X, Y)) + CHECK(ltb(X, Y) + ltb(Y, X) == 0); + else + CHECK(ltb(X, Y) + ltb(Y, X) == 1); + if (ltb(X, Y) && ltb(Y, Z)) + CHECK(ltb(X, Z)); + if (!ltb(Y, X) && !ltb(Z, Y)) + CHECK(!ltb(Z, X)); + + CHECK(((X >> x) << x) == X); + CHECK((X >> x) == (X << -x)); + CHECK((X >> x) == (X << -x)); + CHECK((X >> (x + y)) == ((X >> x) >> y)); + + CHECK((X * x) * y == X * (x * y)); + CHECK((X >> x) * y == (X * y) >> (x * y)); + + CHECK((X.grown(1) == X) == (D == 0 || X.empty())); + if (all(x >= 0 && y >= 0)) + CHECK(X.grown(x).grown(y) == X.grown(x + y)); + else + CHECK((X.grown(x).grown(y).empty() || + X.grown(x).grown(y) == X.grown(x + y)) == true); + if (all(x >= 0)) + CHECK(X.grown(x).grown(-x) == X); + else + CHECK((X.grown(x).grown(-x).empty() || X.grown(x).grown(-x) == X) == + true); + CHECK(X.grown(x) == X.grown(x, x)); + CHECK(X.grown(a) == X.grown(fmap([&](int) { return a; }, x))); + + CHECK(X.shrunk(x, y) == X.grown(-x, -y)); + CHECK(X.shrunk(x) == X.shrunk(x, x)); + CHECK(X.shrunk(a) == X.shrunk(fmap([&](int) { return a; }, x))); + + CHECK(N == N); + CHECK(X == X); + CHECK((N == X) == X.empty()); + CHECK(!(N != N)); + CHECK(!(X != X)); + CHECK((N != X) != (N == X)); + + CHECK(X.contains(X.lower()) == !X.empty()); + CHECK(X.contains(X.upper() - 1) == !X.empty()); + CHECK(X.grown(1).contains(X.upper()) == !X.empty()); + CHECK(isdisjoint(X, X) == X.empty()); + + // a <= b means "a implies b" for booleans (yes, it looks wrong) + CHECK((X < Y) <= (X <= Y)); + CHECK((X > Y) <= (X >= Y)); + CHECK((X <= Y) <= (X.empty() || !isdisjoint(X, Y))); + CHECK((X >= Y) <= (Y.empty() || !isdisjoint(X, Y))); + CHECK(!(X < Y && Y < X)); + CHECK((X <= Y && X >= Y) == (X == Y)); + CHECK((X < X.grown(1)) == (D > 0 && !X.empty())); + CHECK((X.shrunk(1) < X) == (D > 0 && !X.empty())); + + CHECK(N <= N); + CHECK(!(N < N)); + CHECK(N <= X); + CHECK((N < X) == !X.empty()); + + const auto BXY = bounding_box(X, Y); + CHECK(bounding_box(N, X) == X); + CHECK(bounding_box(X, N) == X); + CHECK(bounding_box(X, Y) == bounding_box(Y, X)); + CHECK(bounding_box(bounding_box(X, Y), Z) == + bounding_box(X, bounding_box(Y, Z))); + + CHECK(X <= BXY); + CHECK(Y <= BXY); + CHECK((X.grown(1) <= BXY && Y.grown(1) <= BXY) == (D == 0 || BXY.empty())); + CHECK( + eqb(bounding_box(X.grown(abs(x)), Y.grown(abs(x))), BXY.grown(abs(x)))); + CHECK(eqb(bounding_box(X >> x, Y >> x), BXY >> x)); + CHECK(eqb(bounding_box(X * x, Y * x), BXY * x)); + CHECK( + eqb(bounding_box(X.grown(abs(x)), Y.grown(abs(x))), BXY.grown(abs(x)))); + + const B E = bounding_box(bounding_box(X, Y), Z).grown(10); + + CHECK((N & X) == N); + CHECK((X & N) == N); + CHECK((E & X) == X); + CHECK((X & E) == X); + CHECK((X & Y) == (Y & X)); + CHECK(((X & Y) & Z) == (X & (Y & Z))); + + const B IXY = X & Y; + CHECK((IXY <= X && IXY <= Y) == true); + CHECK((IXY.grown(1) <= X && IXY.grown(1) <= Y) == (D == 0 || IXY.empty())); + + } // for iter +} + +TEST_CASE("Box", "[regions]") { + test_Box(Box()); +} +TEST_CASE("Box", "[regions]") { + test_Box(Box()); +} +TEST_CASE("Box", "[regions]") { + test_Box(Box()); +} +TEST_CASE("Box", "[regions]") { + test_Box(Box()); +} + +TEST_CASE("Box", "[regions]") { test_Box(Box()); } +TEST_CASE("Box", "[regions]") { test_Box(Box()); } +TEST_CASE("Box", "[regions]") { test_Box(Box()); } +TEST_CASE("Box", "[regions]") { test_Box(Box()); } + +TEST_CASE("NDBox(0)", "[regions]") { + test_Box(NDBox(0)); +} +TEST_CASE("NDBox(1)", "[regions]") { + test_Box(NDBox(1)); +} +TEST_CASE("NDBox(2)", "[regions]") { + test_Box(NDBox(2)); +} +TEST_CASE("NDBox(3)", "[regions]") { + test_Box(NDBox(3)); +} + +TEST_CASE("NDBox(0)", "[regions]") { test_Box(NDBox(0)); } +TEST_CASE("NDBox(1)", "[regions]") { test_Box(NDBox(1)); } +TEST_CASE("NDBox(2)", "[regions]") { test_Box(NDBox(2)); } +TEST_CASE("NDBox(3)", "[regions]") { test_Box(NDBox(3)); } diff --git a/test/PointTest.cpp b/test/PointTest.cpp new file mode 100644 index 0000000000..28bb697229 --- /dev/null +++ b/test/PointTest.cpp @@ -0,0 +1,549 @@ +#include +#include + +#include + +#include +#include +#include +#include +#include +#include + +using namespace openPMD::Regions; + +template +auto maxabs(const T &xs) { + if (xs.size() == 0) + return VT(0); + using std::abs, std::max_element; + return max_element(abs(xs)); +} +template > * = nullptr> +auto maxabs(const T &x) { + using std::abs; + return abs(x); +} +template bool is_approx(const T &xs, const T &ys) { + using std::max; + auto m = maxabs(xs - ys); + if (m == 0) + return true; + m /= max(maxabs(xs), maxabs(ys)); + const auto eps = std::numeric_limits::epsilon(); + return m <= 100 * eps; +} + +template bool eq_helper(const T &x1, const T &x2) { + std::equal_to eq; + return eq(x1, x2); +} +template bool eq_helper(const T1 &, const T2 &) { + return false; +} + +template void test_Point_bool(const P &p) { + const std::size_t D = p.ndims(); + using T = typename P::value_type; + const std::equal_to

eq; + const std::less

lt; + + std::mt19937 gen; + std::uniform_int_distribution dist(0, 1); + const auto rand = [&]() { return dist(gen); }; + const auto randp = [&]() { return fmap([&](auto) { return rand(); }, p); }; + + for (int iter = 0; iter < 100; ++iter) { + + P n = p; + CHECK(n.ndims() == D); + for (std::size_t d = 0; d < D; ++d) + CHECK(n[d] == 0); + CHECK(n.size() == D); + + const P x = randp(); + const P y = randp(); + const P z = randp(); + + CHECK(eq(n, n)); + CHECK(eq(x, x)); + CHECK(!lt(n, n)); + CHECK(!lt(x, x)); + if (all(x == n)) { + CHECK(eq(n, x)); + CHECK(!lt(n, x)); + } else { + CHECK(!eq(n, x)); + CHECK(lt(n, x)); + } + + CHECK(!any(n)); + CHECK(all(!n)); + + CHECK(eq(n & x, n)); + CHECK(eq(0 & x, n)); + CHECK(eq(x & n, n)); + CHECK(eq(x & 0, n)); + + CHECK(eq((!n & x), x)); + CHECK(eq((T(1) & x), x)); + CHECK(eq((x & !n), x)); + CHECK(eq((x & T(1)), x)); + + CHECK(eq((n | x), x)); + CHECK(eq((0 | x), x)); + CHECK(eq((x | n), x)); + CHECK(eq((x | 0), x)); + + CHECK(eq((!n | x), !n)); + CHECK(eq((T(1) | x), !n)); + CHECK(eq((x | !n), !n)); + CHECK(eq((x | T(1)), !n)); + + CHECK(eq((x & y), (y & x))); + CHECK(eq((x | y), (y | x))); + + CHECK(eq(((x & y) & z), (x & (y & z)))); + CHECK(eq(((x | y) | z), (x | (y | z)))); + + CHECK(eq((x & (y | z)), ((y & x) | (x & z)))); + CHECK(eq((x | (y & z)), ((y | x) & (x | z)))); + + CHECK(eq((x & y), !(!x | !y))); + CHECK(eq((x | y), !(!x & !y))); + + CHECK(eq((n ^ x), x)); + CHECK(eq((0 ^ x), x)); + CHECK(eq((x ^ n), x)); + CHECK(eq((x ^ 0), x)); + + CHECK(eq((!n ^ x), !x)); + CHECK(eq((T(1) ^ x), !x)); + CHECK(eq((x ^ !n), !x)); + CHECK(eq((x ^ T(1)), !x)); + + CHECK(eq((x ^ x), n)); + + CHECK(eq((x ^ y), (y ^ x))); + CHECK(eq(((x ^ y) ^ z), (x ^ (y ^ z)))); + + CHECK(eq(!(!x), x)); + + CHECK(eq((n && x), n)); + CHECK(eq((0 && x), n)); + CHECK(eq((x && n), n)); + CHECK(eq((x && 0), n)); + + CHECK(eq((!n && x), x)); + CHECK(eq((!T(0) && x), x)); + CHECK(eq((x && !n), x)); + CHECK(eq((x && !T(0)), x)); + + CHECK(eq((n || x), x)); + CHECK(eq((0 || x), x)); + CHECK(eq((x || n), x)); + CHECK(eq((x || 0), x)); + + CHECK(eq((!n || x), !n)); + CHECK(eq((!T(0) || x), !n)); + CHECK(eq((x || !n), !n)); + CHECK(eq((x || !T(0)), !n)); + + CHECK(eq((x && y), (y && x))); + CHECK(eq((x || y), (y || x))); + + CHECK(eq(((x && y) && z), (x && (y && z)))); + CHECK(eq(((x || y) || z), (x || (y || z)))); + + CHECK(eq((x && (y || z)), ((y && x) || (x && z)))); + CHECK(eq((x || (y && z)), ((y || x) && (x || z)))); + + CHECK(eq((x && y), !(!x || !y))); + CHECK(eq((x || y), !(!x && !y))); + + P t; + t = x; + t &= y; + CHECK(eq(t, (x & y))); + t = x; + t |= y; + CHECK(eq(t, (x | y))); + t = x; + t ^= y; + CHECK(eq(t, (x ^ y))); + + } // for iter +} + +template void test_Point_int(const P &p) { + const std::size_t D = p.ndims(); + using T = typename P::value_type; + const std::equal_to

eq; + + std::mt19937 gen; + std::uniform_int_distribution dist(-1000, 1000); + const auto rand = [&]() { return dist(gen); }; + const auto randp = [&]() { return fmap([&](auto) { return rand(); }, p); }; + + for (int iter = 0; iter < 100; ++iter) { + + P n(p); + CHECK(n.size() == D); + for (std::size_t d = 0; d < D; ++d) + CHECK(n[d] == 0); + + const P x = randp(); + const P y = randp(); + const P z = randp(); + + const T a = rand(); + const T b = rand(); + + CHECK(eq(fmap([](auto i) { return i; }, x), x)); + CHECK(eq(fmap([](auto i) { return i + 1; }, + fmap([](auto i) { return 2 * i; }, x)), + fmap([](auto i) { return 2 * i + 1; }, x))); + + CHECK(eq(fmap([](auto i, auto j) { return 2 * i + j; }, x, y), 2 * x + y)); + CHECK(eq( + fmap([](auto i, auto j, auto k) { return 3 * i + 2 * j + k; }, x, y, z), + 3 * x + 2 * y + z)); + + CHECK(fold([](auto i, auto j) { return i + j; }, 0, x) == sum(x)); + CHECK(fold([](auto i, auto j, auto k) { return i + j + k; }, 0, x, y) == + sum(x + y)); + + CHECK(sum(n) == 0); + CHECK(sum(n + 1) == T(D)); + CHECK(product(n) == (D == 0 ? 1 : 0)); + CHECK(product(n + 1) == 1); + CHECK(min_element(n) == (D == 0 ? std::numeric_limits::max() : 0)); + CHECK(max_element(n) == (D == 0 ? std::numeric_limits::min() : 0)); + CHECK(min_element(n + 1) == (D == 0 ? std::numeric_limits::max() : 1)); + CHECK(max_element(n + 1) == (D == 0 ? std::numeric_limits::min() : 1)); + + CHECK(eq(+x, x)); + CHECK(eq(n + x, x)); + CHECK(eq(T(0) + x, x)); + CHECK(eq(x + n, x)); + CHECK(eq(x + T(0), x)); + + CHECK(eq(x + y, y + x)); + + CHECK(eq((x + y) + z, x + (y + z))); + + CHECK(eq(-x, -T(1) * x)); + CHECK(eq(-(-x), x)); + CHECK(eq(x - x, n)); + + CHECK(eq(a * n, n)); + CHECK(eq(n * a, n)); + CHECK(eq(T(0) * x, n)); + CHECK(eq(x * T(0), n)); + CHECK(eq(T(1) * x, x)); + CHECK(eq(x * T(1), x)); + + CHECK(eq(a * x, x * a)); + + CHECK(eq(a * x + b * x, (a + b) * x)); + CHECK(eq(a * (x + y), a * x + a * y)); + + CHECK(eq(x * (y + z), x * y + x * z)); + + if (all(y != 0)) { + CHECK(eq(x * y / y, x)); + CHECK(eq(x / y * y + x % y, x)); + } + + CHECK(eq(~(~x), x)); + + CHECK(eq((n & x), n)); + CHECK(eq((0 & x), n)); + CHECK(eq((x & n), n)); + CHECK(eq((x & 0), n)); + + CHECK(eq((~n & x), x)); + CHECK(eq((~T(0) & x), x)); + CHECK(eq((x & ~n), x)); + CHECK(eq((x & ~T(0)), x)); + + CHECK(eq((n | x), x)); + CHECK(eq((0 | x), x)); + CHECK(eq((x | n), x)); + CHECK(eq((x | 0), x)); + + CHECK(eq((~n | x), ~n)); + CHECK(eq((~T(0) | x), ~n)); + CHECK(eq((x | ~n), ~n)); + CHECK(eq((x | ~T(0)), ~n)); + + CHECK(eq((x & y), (y & x))); + CHECK(eq((x | y), (y | x))); + + CHECK(eq(((x & y) & z), (x & (y & z)))); + CHECK(eq(((x | y) | z), (x | (y | z)))); + + CHECK(eq((x & (y | z)), ((y & x) | (x & z)))); + CHECK(eq((x | (y & z)), ((y | x) & (x | z)))); + + CHECK(eq((x & y), ~(~x | ~y))); + CHECK(eq((x | y), ~(~x & ~y))); + + CHECK(eq((n ^ x), x)); + CHECK(eq((0 ^ x), x)); + CHECK(eq((x ^ n), x)); + CHECK(eq((x ^ 0), x)); + + CHECK(eq((~n ^ x), ~x)); + CHECK(eq((~T(0) ^ x), ~x)); + CHECK(eq((x ^ ~n), ~x)); + CHECK(eq((x ^ ~T(0)), ~x)); + + CHECK(eq((x ^ x), n)); + + CHECK(eq((x ^ y), (y ^ x))); + CHECK(eq(((x ^ y) ^ z), (x ^ (y ^ z)))); + + P t; + t = x; + t += y; + CHECK(eq(t, x + y)); + t = x; + t -= y; + CHECK(eq(t, x - y)); + t = x; + t *= y; + CHECK(eq(t, x * y)); + if (all(y != 0)) { + t = x; + t /= y; + CHECK(eq(t, x / y)); + t = x; + t %= y; + CHECK(eq(t, x % y)); + } + t = x; + t &= y; + CHECK(eq(t, (x & y))); + t = x; + t |= y; + CHECK(eq(t, (x | y))); + t = x; + t ^= y; + CHECK(eq(t, (x ^ y))); + + } // for iter +} + +template void test_Point_float(const P &p) { + const std::size_t D = p.ndims(); + using T = typename P::value_type; + const std::equal_to

eq; + const std::less

lt; + + std::mt19937 gen; + std::uniform_real_distribution dist(-1.0, 1.0); + const auto rand = [&]() { return dist(gen); }; + const auto randp = [&]() { return fmap([&](auto) { return rand(); }, p); }; + + for (int iter = 0; iter < 100; ++iter) { + + P n(p); + CHECK(n.size() == D); + for (std::size_t d = 0; d < D; ++d) + CHECK(n[d] == 0); + + const P x = randp(); + const P y = randp(); + const P z = randp(); + + const T a = rand(); + const T b = rand(); + + CHECK(eq(x, x)); + CHECK(!lt(x, x)); + if (eq(x, y)) + CHECK(lt(x, y) + lt(y, x) == 0); + else + CHECK(lt(x, y) + lt(y, x) == 1); + if (lt(x, y) && lt(y, z)) + CHECK(lt(x, z)); + if (!lt(y, x) && !lt(z, y)) + CHECK(!lt(z, x)); + + // remove-insert is no-op + if (D > 0) { + for (std::size_t d = 0; d < D; ++d) { + const auto a1 = x[d]; + const auto x1 = x.erase(d); + CHECK(x1.ndims() == D - 1); + const auto x2 = x1.insert(d, a1); + CHECK(x2.ndims() == D); + CHECK(eq_helper(x2, x)); + } + } + // insert-remove is no-op + for (std::size_t d = 0; d <= D; ++d) { + const auto x1 = x.insert(d, a); + CHECK(x1.ndims() == D + 1); + CHECK(x1[d] == a); + CHECK(eq(x1.erase(d), x)); + } + + CHECK(eq(x.reversed().reversed(), x)); + + CHECK(eq(fmap([](auto i) { return i; }, x), x)); + CHECK(eq(fmap([](auto i) { return i + 1; }, + fmap([](auto i) { return 2 * i; }, x)), + fmap([](auto i) { return 2 * i + 1; }, x))); + + CHECK(eq(fmap([](auto i, auto j) { return 2 * i + j; }, x, y), 2 * x + y)); + CHECK(eq( + fmap([](auto i, auto j, auto k) { return 3 * i + 2 * j + k; }, x, y, z), + 3 * x + 2 * y + z)); + + CHECK(fold([](auto i, auto j) { return i + j; }, T(0), x) == sum(x)); + CHECK(is_approx( + fold([](auto i, auto j, auto k) { return i + j + k; }, T(0), x, y), + sum(x + y))); + + CHECK(sum(n) == 0); + CHECK(sum(n + 1) == D); + CHECK(product(n) == (D == 0 ? 1 : 0)); + CHECK(product(n + 1) == 1); + const T inf_zero = D == 0 ? T(1) / 0 : 0; + const T neg_inf_zero = D == 0 ? -T(1) / 0 : 0; + const T inf_one = D == 0 ? T(1) / 0 : 1; + const T neg_inf_one = D == 0 ? -T(1) / 0 : 1; + CHECK(min_element(n) == inf_zero); + CHECK(max_element(n) == neg_inf_zero); + CHECK(min_element(n + 1) == inf_one); + CHECK(max_element(n + 1) == neg_inf_one); + + CHECK(eq(+x, x)); + CHECK(eq(n + x, x)); + CHECK(eq(T(0) + x, x)); + CHECK(eq(x + n, x)); + CHECK(eq(x + T(0), x)); + + CHECK(eq(x + y, y + x)); + + CHECK(is_approx((x + y) + z, x + (y + z))); + + CHECK(eq(-x, -T(1) * x)); + CHECK(eq(-(-x), x)); + CHECK(eq(x - x, n)); + + CHECK(eq(a * n, n)); + CHECK(eq(n * a, n)); + CHECK(eq(T(0) * x, n)); + CHECK(eq(x * T(0), n)); + CHECK(eq(T(1) * x, x)); + CHECK(eq(x * T(1), x)); + + CHECK(eq(a * x, x * a)); + + if (all(x != 0)) { + CHECK(eq(x / x, n + 1)); + CHECK(is_approx(1 / (1 / x), x)); + CHECK(is_approx(a / x, a * (1 / x))); + } + if (a != 0) { + CHECK(is_approx(x / a, x * (1 / a))); + } + + CHECK(is_approx(a * x + b * x, (a + b) * x)); + CHECK(is_approx(a * (x + y), a * x + a * y)); + + CHECK(is_approx(x * (y + z), x * y + x * z)); + + if (all(y != 0)) { + CHECK(is_approx(x * y / y, x)); + } + + P t; + t = x; + t += y; + CHECK(eq(t, x + y)); + t = x; + t -= y; + CHECK(eq(t, x - y)); + t = x; + t *= y; + CHECK(eq(t, x * y)); + t = x; + t /= y; + CHECK(eq(t, x / y)); + + } // for iter +} + +TEST_CASE("Point", "[regions]") { test_Point_bool(Point()); } +TEST_CASE("Point", "[regions]") { test_Point_bool(Point()); } +TEST_CASE("Point", "[regions]") { test_Point_bool(Point()); } +TEST_CASE("Point", "[regions]") { test_Point_bool(Point()); } + +TEST_CASE("Point", "[regions]") { + test_Point_int(Point()); +} +TEST_CASE("Point", "[regions]") { + test_Point_int(Point()); +} +TEST_CASE("Point", "[regions]") { + test_Point_int(Point()); +} +TEST_CASE("Point", "[regions]") { + test_Point_int(Point()); +} + +TEST_CASE("Point", "[regions]") { + test_Point_float(Point()); +} +TEST_CASE("Point", "[regions]") { + test_Point_float(Point()); +} +TEST_CASE("Point", "[regions]") { + test_Point_float(Point()); +} +TEST_CASE("Point", "[regions]") { + test_Point_float(Point()); +} + +TEST_CASE("NDPoint(0)", "[regions]") { + test_Point_bool(NDPoint(0)); +} +TEST_CASE("NDPoint(1)", "[regions]") { + test_Point_bool(NDPoint(1)); +} +TEST_CASE("NDPoint(2)", "[regions]") { + test_Point_bool(NDPoint(2)); +} +TEST_CASE("NDPoint(3)", "[regions]") { + test_Point_bool(NDPoint(3)); +} + +TEST_CASE("NDPoint(0)", "[regions]") { + test_Point_int(NDPoint(0)); +} +TEST_CASE("NDPoint(1)", "[regions]") { + test_Point_int(NDPoint(1)); +} +TEST_CASE("NDPoint(2)", "[regions]") { + test_Point_int(NDPoint(2)); +} +TEST_CASE("NDPoint(3)", "[regions]") { + test_Point_int(NDPoint(3)); +} + +TEST_CASE("NDPoint(0)", "[regions]") { + test_Point_float(NDPoint(0)); +} +TEST_CASE("NDPoint(1)", "[regions]") { + test_Point_float(NDPoint(1)); +} +TEST_CASE("NDPoint(2)", "[regions]") { + test_Point_float(NDPoint(2)); +} +TEST_CASE("NDPoint(3)", "[regions]") { + test_Point_float(NDPoint(3)); +} diff --git a/test/RegionTest.cpp b/test/RegionTest.cpp new file mode 100644 index 0000000000..141df92618 --- /dev/null +++ b/test/RegionTest.cpp @@ -0,0 +1,201 @@ +#include +#include + +#include + +#include +#include +#include +#include +#include + +using namespace openPMD::Regions; + +template void test_Region(const R &r) { + const std::size_t D = r.ndims(); + const auto b = bounding_box(r); + const auto p = b.lower(); + using B = std::decay_t; + using P = std::decay_t; + const std::equal_to eqr; + const std::less ltr; + CHECK(r.empty()); + CHECK(b.empty()); + + std::mt19937 gen; + std::uniform_int_distribution dist0(0, 9); + std::uniform_int_distribution dist(-1000, 1000); + const auto rand = [&]() { return dist(gen); }; + const auto randp = [&]() { return fmap([&](auto) { return rand(); }, p); }; + const auto randb = [&]() { + if (D == 0) { + if (dist0(gen) < 5) + return B(p, p); // empty box + else + return B(p); + } + if (dist0(gen) == 0) + return b; + while (1) { + auto lo = randp(); + auto hi = randp(); + auto nb = B(lo, hi); + if (!nb.empty()) + return nb; + } + }; + const auto randr = [&]() { + if (D == 0) { + if (dist0(gen) < 5) + return R(B(p, p)); // empty region + else + return R(B(p)); + } + const int nboxes = dist0(gen) / 2; + R nr(B(p, p)); // empty region + for (int n = 0; n < nboxes; ++n) + nr |= randb(); + return nr; + }; + + for (int iter = 0; iter < 100; ++iter) { + + R N = r; + CHECK(N.empty()); + R X = randr(); + R Y = randr(); + R Z = randr(); + + P n = p; + P x = randp(); + P y = randp(); + + CHECK(eqr(X, X)); + CHECK(!ltr(X, X)); + if (eqr(X, Y)) + CHECK(ltr(X, Y) + ltr(Y, X) == 0); + else + CHECK(ltr(X, Y) + ltr(Y, X) == 1); + if (ltr(X, Y) && ltr(Y, Z)) + CHECK(ltr(X, Z)); + if (!ltr(Y, X) && !ltr(Z, Y)) + CHECK(!ltr(Z, X)); + + CHECK(N == N); + CHECK(X == X); + CHECK((N != X) != X.empty()); + + CHECK((X >> n) == X); + CHECK((X >> x) == (X << -x)); + CHECK(((X >> x) << x) == X); + CHECK(((X >> x) >> y) == (X >> (x + y))); + + CHECK((X * abs(x)) * abs(y) == X * (abs(x) * abs(y))); + CHECK((X >> x) * y == (X * y) >> (x * y)); + + if (X.empty()) { + CHECK(X.grown(abs(x)).empty()); + CHECK(X.shrunk(abs(x)).empty()); + } else if (all(abs(x) == n)) { + CHECK(X.grown(abs(x)) == X); + CHECK(X.shrunk(abs(x)) == X); + } else { + CHECK(X.grown(abs(x)) > X); + CHECK(X.shrunk(abs(x)) < X); + } + + CHECK(X.grown(abs(x)).shrunk(abs(x)) >= X); + CHECK(X.shrunk(abs(x)).grown(abs(x)) <= X); + + CHECK((X.grown(x) >> y) == (X >> y).grown(x)); + + CHECK((X * abs(y)).grown(abs(x) * abs(y)) == X.grown(abs(x)) * abs(y)); + + const B E = bounding_box(bounding_box(bounding_box(X), bounding_box(Y)), + bounding_box(Z)) + .grown(10); + + CHECK((N & X) == N); + CHECK((X & N) == N); + CHECK((E & X) == X); + CHECK((X & E) == X); + CHECK((X & Y) == (Y & X)); + CHECK(((X & Y) & Z) == (X & (Y & Z))); + + CHECK((N | X) == X); + CHECK((E | X) == E); + CHECK((X | E) == E); + CHECK((X | Y) == (Y | X)); + CHECK(((X | Y) | Z) == (X | (Y | Z))); + + CHECK(E - (X & Y) == ((E - X) | (E - Y))); + CHECK(E - (X | Y) == ((E - X) & (E - Y))); + + CHECK((N ^ X) == X); + CHECK((X ^ N) == X); + CHECK((X ^ X) == N); + CHECK((X ^ Y) == (Y ^ X)); + CHECK(((X ^ Y) ^ Z) == (X ^ (Y ^ Z))); + + const R IXY = X & Y; + CHECK((IXY <= X && IXY <= Y) == true); + CHECK((IXY.grown(1) <= X && IXY.grown(1) <= Y) == (D == 0 || IXY.empty())); + + const R UXY = X | Y; + CHECK((X <= UXY && Y <= UXY) == true); + + const R DXY = X - Y; + CHECK((DXY <= X || !isdisjoint(DXY, Y)) == true); + + const R SXY = X ^ Y; + CHECK((SXY <= UXY && isdisjoint(SXY, IXY)) == true); + + CHECK(IXY <= UXY); + CHECK((IXY | SXY) == UXY); + + } // for iter +} + +TEST_CASE("Region", "[regions]") { + test_Region(Region()); +} +TEST_CASE("Region", "[regions]") { + test_Region(Region()); +} +TEST_CASE("Region", "[regions]") { + test_Region(Region()); +} +TEST_CASE("Region", "[regions]") { + test_Region(Region()); +} + +TEST_CASE("Region", "[regions]") { test_Region(Region()); } +TEST_CASE("Region", "[regions]") { test_Region(Region()); } +TEST_CASE("Region", "[regions]") { test_Region(Region()); } +TEST_CASE("Region", "[regions]") { test_Region(Region()); } + +TEST_CASE("NDRegion(0)", "[regions]") { + test_Region(NDRegion(0)); +} +TEST_CASE("NDRegion(1)", "[regions]") { + test_Region(NDRegion(1)); +} +TEST_CASE("NDRegion(2)", "[regions]") { + test_Region(NDRegion(2)); +} +TEST_CASE("NDRegion(3)", "[regions]") { + test_Region(NDRegion(3)); +} + +TEST_CASE("NDRegion(0)", "[regions]") { + test_Region(NDRegion(0)); +} +TEST_CASE("NDRegion(1)", "[regions]") { + test_Region(NDRegion(1)); +} +TEST_CASE("NDRegion(2)", "[regions]") { + test_Region(NDRegion(2)); +} +TEST_CASE("NDRegion(3)", "[regions]") { + test_Region(NDRegion(3)); +}