aboutsummaryrefslogtreecommitdiff
path: root/Carpet/CarpetLib/src/bboxset2.hh
diff options
context:
space:
mode:
Diffstat (limited to 'Carpet/CarpetLib/src/bboxset2.hh')
-rw-r--r--Carpet/CarpetLib/src/bboxset2.hh1369
1 files changed, 1369 insertions, 0 deletions
diff --git a/Carpet/CarpetLib/src/bboxset2.hh b/Carpet/CarpetLib/src/bboxset2.hh
new file mode 100644
index 000000000..bac295377
--- /dev/null
+++ b/Carpet/CarpetLib/src/bboxset2.hh
@@ -0,0 +1,1369 @@
+#ifndef BBOXSET2_HH
+#define BBOXSET2_HH
+
+#include <algorithm>
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <limits>
+#include <map>
+#include <memory>
+#include <set>
+#include <sstream>
+#include <utility>
+#include <vector>
+
+#include <cctk.h>
+
+#include "bbox.hh"
+#include "defs.hh"
+#include "vect.hh"
+
+using namespace std;
+
+
+
+#ifndef CARPET_NO_BBOXSET2
+
+namespace bboxset2 {
+
+
+
+template<typename T, int D>
+class bboxset {
+ template<typename, int> friend class bboxset;
+ typedef ::vect<T,D> vect;
+ typedef ::bbox<T,D> bbox;
+ typedef ::vect<T,D-1> vect1;
+ typedef ::bbox<T,D-1> bbox1;
+ typedef bboxset<T,D-1> bboxset1;
+
+#if 0 && !defined __ICC
+ template<typename U>
+ using ptr = unique_ptr<U>;
+#else
+ template<typename U>
+ using ptr = shared_ptr<U>;
+#endif
+
+ typedef map<T, ptr<bboxset1>> subsets_t;
+ subsets_t subsets;
+
+ vect stride, offset;
+
+ bool is_poison_;
+
+ template<typename F>
+ void traverse_subsets(const F& f) const
+ {
+ bboxset1 decoded_subset;
+ assert(decoded_subset.empty());
+ for (const auto& pos_subset: subsets) {
+ const T& pos = pos_subset.first;
+ const bboxset1* const subsetp = pos_subset.second.get();
+ decoded_subset ^= *subsetp;
+ f(pos, decoded_subset);
+ }
+ assert(decoded_subset.empty());
+ }
+
+ template<typename F>
+ void traverse_subsets(const F& f, const bboxset& other) const
+ {
+ bboxset1 decoded_subset0;
+ bboxset1 decoded_subset1;
+ assert(decoded_subset0.empty());
+ assert(decoded_subset1.empty());
+
+ typedef typename subsets_t::const_iterator subsets_iter_t;
+ subsets_iter_t iter0 = subsets.begin();
+ subsets_iter_t iter1 = other.subsets.begin();
+ subsets_iter_t const end0 = subsets.end();
+ subsets_iter_t const end1 = other.subsets.end();
+ while (iter0!=end0 or iter1!=end1) {
+ const T next_pos0 = iter0!=end0 ? iter0->first : numeric_limits<T>::max();
+ const T next_pos1 = iter1!=end1 ? iter1->first : numeric_limits<T>::max();
+ const T pos = min(next_pos0, next_pos1);
+ const bool active0 = next_pos0==pos;
+ const bool active1 = next_pos1==pos;
+ const bboxset1* const subset0p = active0 ? iter0->second.get() : 0;
+ const bboxset1* const subset1p = active1 ? iter1->second.get() : 0;
+ if (active0) decoded_subset0 ^= *subset0p;
+ if (active1) decoded_subset1 ^= *subset1p;
+
+ f(pos, decoded_subset0, decoded_subset1);
+
+ if (active0) ++iter0;
+ if (active1) ++iter1;
+ }
+ assert(decoded_subset0.empty());
+ assert(decoded_subset1.empty());
+ }
+
+ template<typename F>
+ bboxset binary_operator(const F& op, const bboxset& other) const
+ {
+ assert(not is_poison() and not other.is_poison());
+
+#if 0
+ // TODO: This assumes that the empty set is a neutral element for
+ // the operator
+ if (other.empty()) return *this;
+ if (empty()) return other;
+ assert(all(other.stride == stride));
+ assert(all(offset == other.offset));
+#endif
+
+ if (not empty() and not other.empty()) {
+ assert(all(stride == other.stride));
+ // TODO
+ if (not (all(offset == other.offset))) {
+ cout << "*this=" << *this << "\n"
+ << "other=" << other << "\n";
+ }
+ assert(all(offset == other.offset));
+ }
+
+ bboxset res;
+ if (empty()) {
+ res.stride = other.stride;
+ res.offset = other.offset;
+ } else {
+ res.stride = stride;
+ res.offset = offset;
+ }
+
+ bboxset1 old_decoded_subsetr;
+ traverse_subsets
+ ([&](const T& pos,
+ const bboxset1& decoded_subset0,
+ const bboxset1& decoded_subset1)
+ {
+ bboxset1 decoded_subsetr = op(decoded_subset0, decoded_subset1);
+ const auto subsetrp =
+ make_shared<bboxset1>(decoded_subsetr ^ old_decoded_subsetr);
+ if (not subsetrp->empty()) {
+ res.subsets.insert(res.subsets.end(), make_pair(pos, subsetrp));
+ }
+ swap(old_decoded_subsetr, decoded_subsetr);
+ },
+ other);
+ assert(old_decoded_subsetr.empty());
+
+ return res;
+ }
+
+public:
+
+ /** Invariant */
+ bool invariant() const
+ {
+ if (is_poison() and not empty()) return false;
+ if (any(stride <= vect(0))) return false;
+ if (any(offset < vect(0) or offset >= stride)) return false;
+ for (const auto& pos_subset: subsets) {
+ if (pos_subset.second.get()->empty()) return false;
+ if (any(pos_subset.second.get()->stride != init(stride))) return false;
+ if (any(pos_subset.second.get()->offset != init(offset))) return false;
+ }
+ if (chi_size() % 2 != 0) return false;
+ return true;
+ }
+
+
+
+ /** Create empty set */
+ bboxset(): stride(vect(1)), offset(vect(0)), is_poison_(false)
+ {
+ }
+
+ /** Return a is_poisoned set */
+ static bboxset poison()
+ {
+ bboxset bs;
+ bs.is_poison_ = true;
+ return bs;
+ }
+
+ /** Copy constructor */
+ bboxset(const bboxset& other);
+
+ /** Copy constructor */
+ // bboxset(bboxset&& other);
+
+ /** Assignment */
+ bboxset& operator=(const bboxset& other);
+
+ /** Assignment */
+ // bboxset& operator=(bboxset&& other);
+
+ /** Create set from bbox */
+ bboxset(const bbox& b);
+
+ /** Create set from container of bboxes, bboxsets, or container
+ thereof */
+ template<typename C>
+ bboxset(const C& elts);
+
+ /** Create set from container of structs containing a bbox, bboxset,
+ or container thereof */
+ template<typename C, typename S, typename B>
+ bboxset(const C& elts, const B S::* const mptr);
+
+
+
+ /** Test for is_poison */
+ bool is_poison() const
+ {
+ return is_poison_;
+ }
+
+ /** Test for emptiness */
+ bool empty() const
+ {
+ assert(not is_poison());
+ return subsets.empty();
+ }
+
+ /** Find a point contained in this set */
+ vect front() const;
+
+ /** Number of characteristic points */
+ T chi_size() const;
+
+ /** Number of elements */
+ T size() const;
+
+ /** Container (min and max) */
+ bbox container() const;
+
+
+
+ /** Test for equality */
+ bool operator==(const bboxset& other) const;
+
+ /** Test for is-subset-of */
+ bool operator<=(const bboxset& other) const;
+
+ /** Test for is-strict-subset-of */
+ bool operator<(const bboxset& other) const;
+
+ bool operator!=(const bboxset& other) const
+ {
+ return not (*this == other);
+ }
+
+ bool operator>=(const bboxset& other) const
+ {
+ return other <= *this;
+ }
+
+ bool operator>(const bboxset& other) const
+ {
+ return other < *this;
+ }
+
+ bool contains(const vect& v) const
+ {
+ return bbox(v, v, stride) <= *this;
+ }
+
+ bool intersects(const bbox& b) const
+ {
+ return not (b & *this).empty();
+ }
+
+
+
+ /** Symmetric set difference */
+ bboxset operator^(const bboxset& other) const;
+
+ /** Symmetric set difference */
+ bboxset& operator^=(const bboxset& other);
+
+ /** Set intersection */
+ bboxset operator&(const bboxset& other) const;
+
+ /** Set intersection */
+ bboxset& operator&=(const bboxset& other);
+
+ /** Set Union */
+ bboxset operator|(const bboxset& other) const;
+
+ /** Set union */
+ bboxset& operator|=(const bboxset& other);
+
+ /** Symmetric set union */
+ bboxset operator+(const bboxset& other) const;
+
+ /** Symmetric set union */
+ bboxset& operator+=(const bboxset& other);
+
+ /** Set difference */
+ bboxset operator-(const bboxset& other) const;
+
+ /** Set difference */
+ bboxset& operator-=(const bboxset& other);
+
+
+
+ /** Operators taking bboxes as arguments */
+ bboxset& operator&=(const bbox& b)
+ {
+ return *this &= bboxset(b);
+ }
+ bboxset& operator^=(const bbox& b)
+ {
+ return *this ^= bboxset(b);
+ }
+ bboxset& operator|=(const bbox& b)
+ {
+ return *this |= bboxset(b);
+ }
+ bboxset& operator+=(const bbox& b)
+ {
+ return *this += bboxset(b);
+ }
+ bboxset& operator-=(const bbox& b)
+ {
+ return *this -= bboxset(b);
+ }
+
+
+
+ /** Shift all points */
+ bboxset shift(const vect& dist, const vect& dist_denom = vect(1)) const;
+
+ /** Expand the set (convolute with a bbox) */
+ bboxset expand(const vect& lo, const vect& hi) const;
+
+ /** Expand the set (convolute with a bbox) */
+ bboxset expand(const ::vect<vect,2>& lohi) const
+ {
+ return expand(lohi[0], lohi[1]);
+ }
+
+
+
+private:
+ static bool strides_are_compatible(const vect& str1, const vect& str2);
+public:
+
+ /** Expande the set (changing the stride): find the smallest
+ b-compatible bboxset around this bboxset */
+ // Possible implementation:
+ // - convert bboxset into set of bboxes
+ // - expand each bbox
+ // - create new bboxset as union of these bboxes
+ // Alternative, possibly more efficient implementation:
+ // - for each dimension:
+ // - BS1 := shifted left to align with bbox
+ // - BS2 := shifted right to align with bbox
+ // - create union of BS1 and BS2
+ bboxset expanded_for(const bbox& target) const;
+
+ /** Shrink the set (changing the stride): find the largest
+ b-compatible bboxset inside this bboxset */
+ bboxset contracted_for(const bbox& target) const;
+
+ /** Expande the set (changing the stride): find the smallest open
+ b-compatible bboxset around this bboxset */
+ bboxset anti_contracted_for(const bbox& target) const;
+
+
+
+ /** Serialise the set */
+ template<typename C>
+ void serialise(C& out) const;
+
+ /** Iterate over a serialised set */
+private:
+ typedef vector<bbox> iter_memo_t;
+public:
+ typedef typename iter_memo_t::const_iterator const_iterator;
+private:
+ mutable iter_memo_t iter_memo;
+public:
+ int setsize() const
+ {
+ iter_memo_t im;
+ serialise(im);
+ return im.size();
+ }
+ const_iterator begin() const
+ {
+ iter_memo.clear();
+ serialise(iter_memo);
+ return iter_memo.begin();
+ }
+ const_iterator end() const
+ {
+ return iter_memo.end();
+ }
+
+
+ /** Memory usage */
+ size_t memory() const;
+
+ /** Input */
+ istream& input(istream& is);
+
+ /** Output */
+ ostream& output(ostream& os) const;
+};
+
+
+
+template<typename T>
+class bboxset<T,0> {
+ template<typename, int> friend class bboxset;
+ typedef ::vect<T,0> vect;
+ typedef ::bbox<T,0> bbox;
+
+ bool state;
+ vect stride, offset;
+ bool is_poison_;
+public:
+ bool invariant() const
+ {
+ if (is_poison() and state) return false;
+ return true;
+ }
+
+ bboxset():
+ state(false), is_poison_(false)
+ {
+ }
+ static bboxset poison()
+ {
+ bboxset bs; bs.is_poison_ = true; return bs;
+ }
+ bboxset(const bbox& b):
+ state(true), is_poison_(false)
+ {
+ assert(not b.is_poison());
+ }
+ bboxset(const bboxset& other):
+ state(other.state), is_poison_(other.is_poison_)
+ {
+ }
+ bboxset& operator=(const bboxset& other)
+ {
+ state = other.state;
+ is_poison_ = other.is_poison_;
+ return *this;
+ }
+
+ bool is_poison() const
+ {
+ return is_poison_;
+ }
+ bool empty() const
+ {
+ assert(not is_poison());
+ return not state;
+ }
+ vect front() const
+ {
+ assert(not is_poison());
+ assert(not empty());
+ return vect();
+ }
+ T chi_size() const
+ {
+ assert(not is_poison());
+ return 1;
+ }
+ T size() const
+ {
+ assert(not is_poison());
+ return state;
+ }
+ bbox container() const
+ {
+ assert(not is_poison());
+ return bbox();
+ }
+
+ bool operator==(const bboxset& other) const
+ {
+ assert(not is_poison() and not other.is_poison());
+ return state == other.state;
+ }
+ bool operator<=(const bboxset& other) const
+ {
+ assert(not is_poison() and not other.is_poison());
+ return state <= other.state;
+ }
+ bool operator!=(const bboxset& other) const
+ {
+ return not (*this == other);
+ }
+ bool operator>=(const bboxset& other) const
+ {
+ return other <= *this;
+ }
+ bool operator<(const bboxset& other) const
+ {
+ return not (*this >= other);
+ }
+ bool operator>(const bboxset& other) const
+ {
+ return not (*this <= other);
+ }
+
+ bboxset& operator&=(const bboxset& other)
+ {
+ assert(not is_poison() and not other.is_poison());
+ state &= other.state;
+ return *this;
+ }
+ bboxset operator&(const bboxset& other) const
+ {
+ bboxset res = *this;
+ return res &= other;
+ }
+ bboxset& operator^=(const bboxset& other)
+ {
+ assert(not is_poison() and not other.is_poison());
+ state ^= other.state;
+ return *this;
+ }
+ bboxset operator^(const bboxset& other) const
+ {
+ bboxset res = *this;
+ return res ^= other;
+ }
+ bboxset& operator|=(const bboxset& other)
+ {
+ assert(not is_poison() and not other.is_poison());
+ state |= other.state;
+ return *this;
+ }
+ bboxset operator|(const bboxset& other) const
+ {
+ bboxset res = *this;
+ return res |= other;
+ }
+ bboxset& operator+=(const bboxset& other)
+ {
+ assert(not is_poison() and not other.is_poison());
+ assert(not (state and other.state));
+ return *this |= other;
+ }
+ bboxset operator+(const bboxset& other) const
+ {
+ bboxset res = *this;
+ return res += other;
+ }
+ bboxset& operator-=(const bboxset& other)
+ {
+ assert(not is_poison() and not other.is_poison());
+ state &= not other.state;
+ return *this;
+ }
+ bboxset operator-(const bboxset& other) const
+ {
+ bboxset res = *this;
+ return res -= other;
+ }
+
+ bboxset& operator&=(const bbox& b)
+ {
+ return *this &= bboxset(b);
+ }
+ bboxset& operator^=(const bbox& b)
+ {
+ return *this ^= bboxset(b);
+ }
+ bboxset& operator|=(const bbox& b)
+ {
+ return *this |= bboxset(b);
+ }
+ bboxset& operator+=(const bbox& b)
+ {
+ return *this += bboxset(b);
+ }
+ bboxset& operator-=(const bbox& b)
+ {
+ return *this -= bboxset(b);
+ }
+
+ bboxset shift(const vect& dist, const vect& dist_denom = vect(1)) const
+ {
+ assert(not is_poison());
+ return *this;
+ }
+ bboxset expand(const vect& lo, const vect& hi) const
+ {
+ assert(not is_poison());
+ return *this;
+ }
+ bboxset expand(const ::vect<vect,2>& lohi) const
+ {
+ return expand(lohi[0], lohi[1]);
+ }
+
+ template<typename C>
+ void serialise(C& out) const
+ {
+ assert(not is_poison());
+ if (state) out.insert(out.end(), bbox());
+ }
+
+ size_t memory() const
+ {
+ return sizeof *this;
+ }
+
+ istream& input(istream& is);
+
+ ostream& output(ostream& os) const;
+};
+
+
+
+/** Operators involving bboxes */
+template<typename T, int D>
+inline bboxset<T,D> operator&(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs & bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator&(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) & bs;
+}
+// Note: bbox & bbox -> bbox
+
+template<typename T, int D>
+inline bboxset<T,D> operator^(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs ^ bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator^(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) ^ bs;
+}
+template<typename T, int D>
+inline bboxset<T,D> operator^(const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+ return bboxset<T,D>(b1) ^ bboxset<T,D>(b2);
+}
+
+template<typename T, int D>
+inline bboxset<T,D> operator|(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs | bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator|(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) | bs;
+}
+template<typename T, int D>
+inline bboxset<T,D> operator|(const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+ return bboxset<T,D>(b1) | bboxset<T,D>(b2);
+}
+
+template<typename T, int D>
+inline bboxset<T,D> operator+(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs + bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator+(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) + bs;
+}
+template<typename T, int D>
+inline bboxset<T,D> operator+(const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+ return bboxset<T,D>(b1) + bboxset<T,D>(b2);
+}
+
+template<typename T, int D>
+inline bboxset<T,D> operator-(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs - bboxset<T,D>(b);
+}
+template<typename T, int D>
+inline bboxset<T,D> operator-(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) - bs;
+}
+template<typename T, int D>
+inline bboxset<T,D> operator-(const bbox<T,D>& b1, const bbox<T,D>& b2)
+{
+ return bboxset<T,D>(b1) - bboxset<T,D>(b2);
+}
+
+template<typename T, int D>
+inline bool operator==(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) == bs;
+}
+template<typename T, int D>
+inline bool operator==(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs == bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator!=(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) != bs;
+}
+template<typename T, int D>
+inline bool operator!=(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs != bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator<=(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) <= bs;
+}
+template<typename T, int D>
+inline bool operator<=(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs <= bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator<(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) < bs;
+}
+template<typename T, int D>
+inline bool operator<(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs < bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator>=(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) >= bs;
+}
+template<typename T, int D>
+inline bool operator>=(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs >= bboxset<T,D>(b);
+}
+
+template<typename T, int D>
+inline bool operator>(const bbox<T,D>& b, const bboxset<T,D>& bs)
+{
+ return bboxset<T,D>(b) > bs;
+}
+template<typename T, int D>
+inline bool operator>(const bboxset<T,D>& bs, const bbox<T,D>& b)
+{
+ return bs > bboxset<T,D>(b);
+}
+
+
+
+/** Memory usage */
+template<typename T, int D>
+inline size_t memoryof(bboxset<T,D> const& bs)
+{
+ return bs.memory();
+}
+
+/** Input */
+template<typename T, int D>
+inline istream& operator>>(istream& is, bboxset<T,D>& bs)
+{
+ return bs.input(is);
+}
+
+/** Output */
+template<typename T, int D>
+inline ostream& operator<<(ostream& os, const bboxset<T,D>& bs)
+{
+ return bs.output(os);
+}
+
+
+
+ //////////////////////////////////////////////////////////////////////////////
+ // bboxset<T,D>
+ //////////////////////////////////////////////////////////////////////////////
+
+ /** Copy constructor */
+ template<typename T, int D>
+ bboxset<T,D>::bboxset(const bboxset& other):
+ stride(other.stride), offset(other.offset), is_poison_(other.is_poison_)
+ {
+ for (const auto& pos_subset: other.subsets) {
+ const T& pos = pos_subset.first;
+ const auto new_subsetp = make_shared<bboxset1>(*pos_subset.second.get());
+ subsets.insert(subsets.end(), make_pair(pos, new_subsetp));
+ }
+ }
+
+ /** Copy constructor */
+ // bboxset(bboxset&& other): stride(other.stride), offset(other.offset)
+ // {
+ // swap(subsets, other.subsets);
+ // }
+
+ /** Assignment */
+ template<typename T, int D>
+ bboxset<T,D>& bboxset<T,D>::operator=(const bboxset& other)
+ {
+ if (&other == this) return *this;
+ subsets.clear();
+ for (const auto& pos_subset: other.subsets) {
+ const T& pos = pos_subset.first;
+ const auto new_subsetp = make_shared<bboxset1>(*pos_subset.second.get());
+ subsets.insert(subsets.end(), make_pair(pos, new_subsetp));
+ }
+ stride = other.stride;
+ offset = other.offset;
+ is_poison_ = other.is_poison_;
+ return *this;
+ }
+
+ /** Assignment */
+ // bboxset& operator=(bboxset&& other)
+ // {
+ // if (&other == this) return *this;
+ // swap(subsets, other.subsets);
+ // stride = other.stride;
+ // offset = other.offset;
+ // return *this;
+ // }
+
+ /** Create set from bbox */
+ template<typename T, int D>
+ bboxset<T,D>::bboxset(const bbox& b):
+ stride(b.stride()), offset(imod(b.lower(), b.stride())), is_poison_(false)
+ {
+ assert(not b.is_poison());
+ if (b.empty()) return;
+ const T lo = last(b.lower());
+ const T hi = last(b.upper()) + last(b.stride());
+ const bbox1 b1(init(b.lower()), init(b.upper()), init(b.stride()));
+ const auto lo_subsetp = make_shared<bboxset1>(b1);
+ const auto hi_subsetp = make_shared<bboxset1>(b1);
+ // subsets.emplace_hint(subsets.end(), lo, lo_subsetp);
+ // subsets.emplace_hint(subsets.end(), hi, hi_subsetp);
+ subsets.insert(subsets.end(), make_pair(lo, lo_subsetp));
+ subsets.insert(subsets.end(), make_pair(hi, hi_subsetp));
+ }
+
+ /** Create set from container of bboxes, bboxsets, or container
+ thereof */
+ template<typename T, int D>
+ template<typename C>
+ bboxset<T,D>::bboxset(const C& elts)
+ {
+ *this = bboxset();
+ for (const auto& elt: elts) {
+ *this += bboxset(elt);
+ }
+ }
+
+ /** Create set from container of structs containing a bbox, bboxset,
+ or container thereof */
+ template<typename T, int D>
+ template<typename C, typename S, typename B>
+ bboxset<T,D>::bboxset(const C& elts, const B S::* const mptr)
+ {
+ *this = bboxset();
+ for (const auto& elt: elts) {
+ *this += bboxset(elt.*mptr);
+ }
+ }
+
+
+
+ /** Find a point contained in this set */
+ template<typename T, int D>
+ vect<T,D> bboxset<T,D>::front() const
+ {
+ assert(not is_poison());
+ assert(not empty());
+ const auto& pos_subset = *subsets.begin();
+ const auto& pos = pos_subset.first;
+ const auto& subset = *pos_subset.second.get();
+ const vect1 f1 = subset.front();
+ return vect(f1, pos);
+ }
+
+ /** Number of characteristic points */
+ template<typename T, int D>
+ T bboxset<T,D>::chi_size() const
+ {
+ assert(not is_poison());
+ T sum = 0;
+ for (const auto& pos_subset: subsets) sum += pos_subset.second->chi_size();
+ return sum;
+ }
+
+ /** Number of elements */
+ template<typename T, int D>
+ T bboxset<T,D>::size() const
+ {
+ assert(not is_poison());
+ T total_size = 0; // accumulated total number of points
+ T old_pos = numeric_limits<T>::min(); // location of last subset
+ T old_subset_size = 0; // number of points in the last subset
+ traverse_subsets
+ ([&](const T& pos, const bboxset1& subset)
+ {
+ const T subset_size = subset.size();
+ total_size += (pos - old_pos) * old_subset_size;
+ old_pos = pos;
+ old_subset_size = subset_size;
+ });
+ assert(old_subset_size == 0);
+ return total_size;
+ }
+
+ /** Container (min and max) */
+ template<typename T, int D>
+ bbox<T,D> bboxset<T,D>::container() const
+ {
+ assert(not is_poison());
+ if (empty()) return bbox();
+ const T lo = subsets.begin()->first;
+ const T hi = subsets.rbegin()->first;
+ bbox1 container1;
+ traverse_subsets
+ ([&](const T& pos, const bboxset1& subset)
+ {
+ container1 = container1.expanded_containing(subset.container());
+ });
+ return bbox(vect(container1.lower(), lo),
+ vect(container1.upper(), hi - last(stride)),
+ vect(container1.stride(), last(stride)));
+ }
+
+
+
+ /** Test for equality */
+ template<typename T, int D>
+ bool bboxset<T,D>::operator==(const bboxset& other) const
+ {
+ return (*this ^ other).empty();
+ }
+
+ /** Test for is-subset-of */
+ template<typename T, int D>
+ bool bboxset<T,D>::operator<=(const bboxset& other) const
+ {
+ return (*this | other) == other;
+ }
+
+ /** Test for is-strict-subset-of */
+ template<typename T, int D>
+ bool bboxset<T,D>::operator<(const bboxset& other) const
+ {
+ return *this != other and *this <= other;
+ }
+
+
+
+ /** Symmetric set difference */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::operator^(const bboxset& other) const
+ {
+ // TODO: If other is much smaller than this, direct insertion may
+ // be faster
+ return binary_operator
+ ([](const bboxset1& set0, const bboxset1& set1) { return set0 ^ set1; },
+ other);
+ }
+
+ /** Symmetric set difference */
+ template<typename T, int D>
+ bboxset<T,D>& bboxset<T,D>::operator^=(const bboxset& other)
+ {
+ bboxset res = *this;
+ return *this = *this ^ other;
+ }
+
+ /** Set intersection */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::operator&(const bboxset& other) const
+ {
+ return binary_operator
+ ([](const bboxset1& set0, const bboxset1& set1) { return set0 & set1; },
+ other);
+ }
+
+ /** Set intersection */
+ template<typename T, int D>
+ bboxset<T,D>& bboxset<T,D>::operator&=(const bboxset& other)
+ {
+ return *this = *this & other;
+ }
+
+ /** Set Union */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::operator|(const bboxset& other) const
+ {
+ return binary_operator
+ ([](const bboxset1& set0, const bboxset1& set1) { return set0 | set1; },
+ other);
+ }
+
+ /** Set union */
+ template<typename T, int D>
+ bboxset<T,D>& bboxset<T,D>::operator|=(const bboxset& other)
+ {
+ return *this = *this | other;
+ }
+
+ /** Symmetric set union */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::operator+(const bboxset& other) const
+ {
+ // return binary_operator
+ // ([](const bboxset1& set0, const bboxset1& set1) { return set0 + set1; },
+ // other);
+#ifdef CARPET_DEBUG
+ assert((*this & other).empty());
+#endif
+ return *this ^ other;
+ }
+
+ /** Symmetric set union */
+ template<typename T, int D>
+ bboxset<T,D>& bboxset<T,D>::operator+=(const bboxset& other)
+ {
+ return *this = *this + other;
+ }
+
+ /** Set difference */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::operator-(const bboxset& other) const
+ {
+ return binary_operator
+ ([](const bboxset1& set0, const bboxset1& set1) { return set0 - set1; },
+ other);
+ }
+
+ /** Set difference */
+ template<typename T, int D>
+ bboxset<T,D>& bboxset<T,D>::operator-=(const bboxset& other)
+ {
+ return *this = *this - other;
+ }
+
+
+
+ /** Shift all points */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::shift(const vect& dist, const vect& dist_denom)
+ const
+ {
+ assert(not is_poison());
+ assert(all(stride % dist_denom == 0));
+ if (all(dist == 0)) return *this;
+ bboxset res;
+ res.stride = stride;
+ res.offset = imod(offset + dist * stride / dist_denom, res.stride);
+ for (const auto& pos_subset: subsets) {
+ const T& pos = pos_subset.first;
+ const bboxset1& subset = *pos_subset.second.get();
+ const T new_pos = pos + last(dist * stride / dist_denom);
+ const auto new_subsetp =
+ make_shared<bboxset1>(subset.shift(init(dist), init(dist_denom)));
+ res.subsets.insert(res.subsets.end(), make_pair(new_pos, new_subsetp));
+ }
+ return res;
+ }
+
+ /** Expand the set (convolute with a bbox) */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>:: expand(const vect& lo, const vect& hi) const
+ {
+ assert(not is_poison());
+ assert(all(lo >= 0 and hi >= 0));
+ bboxset res = shift(-lo);
+ for (int d=0; d<D; ++d) {
+ T to_expand = (hi + lo)[d];
+ T current_size = 1;
+ while (to_expand > 0) {
+ const T this_expand = min(to_expand, current_size);
+ res |= res.shift(vect::dir(d) * this_expand);
+ current_size += this_expand;
+ to_expand -= this_expand;
+ }
+ assert(to_expand == 0);
+ }
+ return res;
+ }
+
+
+
+ template<typename T, int D>
+ bool bboxset<T,D>::strides_are_compatible(const vect& str1, const vect& str2)
+ {
+ if (all(str1 >= str2)) {
+ return all(imod(str1, str2) == 0);
+ } else if (all(str1 <= str2)) {
+ return all(imod(str2, str1) == 0);
+ }
+ return false;
+ }
+
+ /** Expande the set (changing the stride): find the smallest
+ b-compatible bboxset around this bboxset */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::expanded_for(const bbox& target) const
+ {
+ // Check preconditions
+ assert(not is_poison() and not target.is_poison());
+ assert(not target.empty());
+ assert(strides_are_compatible(stride, target.stride()));
+ vector<bbox> bs;
+ serialise(bs);
+ bboxset res;
+ for (const auto& b: bs) {
+ res |= b.expanded_for(target);
+ }
+ return res;
+ }
+
+ /** Shrink the set (changing the stride): find the largest
+ b-compatible bboxset inside this bboxset */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::contracted_for(const bbox& target) const
+ {
+ // Check preconditions
+ assert(not is_poison() and not target.is_poison());
+ assert(not target.empty());
+ assert(strides_are_compatible(stride, target.stride()));
+ const bbox cont = container();
+ const vect safety = 10;
+ const bbox good_world = cont.anti_contracted_for(target).expand(safety);
+ const bbox world1 = good_world.anti_contracted_for(cont).expand(safety);
+ const bbox world2 = world1.anti_contracted_for(target).expand(safety);
+ return (world2 ^ (world1 ^ *this).anti_contracted_for(target)) & good_world;
+ }
+
+ /** Expande the set (changing the stride): find the smallest open
+ b-compatible bboxset around this bboxset */
+ template<typename T, int D>
+ bboxset<T,D> bboxset<T,D>::anti_contracted_for(const bbox& target) const
+ {
+ // Check preconditions
+ assert(not is_poison() and not target.is_poison());
+ assert(not target.empty());
+ assert(strides_are_compatible(stride, target.stride()));
+ vector<bbox> bs;
+ serialise(bs);
+ bboxset res;
+ for (const auto& b: bs) {
+ res |= b.anti_contracted_for(target);
+ }
+ return res;
+ }
+
+
+
+ /** Serialise the set */
+ template<typename T, int D>
+ template<typename C>
+ void bboxset<T,D>::serialise(C& out) const
+ {
+ assert(not is_poison());
+ typedef map<bbox1,T> subboxes_t;
+ typedef set<bbox1> subboxes1_t;
+ // TODO: Instead of copying from old_subboxes to subboxes,
+ // maintain subboxes via a non-const iterator
+ subboxes_t old_subboxes;
+ traverse_subsets
+ ([&](const T& pos, const bboxset1& subset)
+ {
+ // Convert subset to bboxes
+ subboxes1_t subboxes1;
+ subset.serialise(subboxes1);
+
+ const bbox1 dummy1;
+ subboxes_t subboxes;
+ // subboxes.reserve(old_subboxes.size() + subboxes1.size());
+ typedef typename subboxes_t::const_iterator subboxes_iter_t;
+ typedef typename subboxes1_t::const_iterator subboxes1_iter_t;
+ subboxes_iter_t iter0 = old_subboxes.begin();
+ subboxes1_iter_t iter1 = subboxes1.begin();
+ subboxes_iter_t const end0 = old_subboxes.end();
+ subboxes1_iter_t const end1 = subboxes1.end();
+ while (iter0!=end0 or iter1!=end1) {
+ bool active0 = iter0!=end0;
+ bool active1 = iter1!=end1;
+ const bbox1& subbox0 = active0 ? iter0->first : dummy1;
+ const bbox1& subbox1 = active1 ? *iter1 : dummy1;
+ // When both subboxes are active, keep only the first (as
+ // determined by less<>)
+ if (active0 and active1) {
+ if (subbox0 != subbox1) {
+ /*const*/ less<bbox1> bbox1_less;
+ if (bbox1_less(subbox0, subbox1)) {
+ active1 = false;
+ } else {
+ active0 = false;
+ }
+ }
+ }
+ const T old_pos = active0 ? iter0->second : T();
+
+ if ((active0 and active1) and (subbox0 == subbox1)) {
+ // The current bbox continues unchanged -- keep it
+ subboxes.insert(subboxes.end(), *iter0);
+ } else {
+ if (active0) {
+ // The current box changed; finalize it
+ const bbox new_box(vect(subbox0.lower(), old_pos),
+ vect(subbox0.upper(), pos - last(stride)),
+ stride);
+ out.insert(out.end(), new_box);
+ }
+ if (active1) {
+ // There is a new box; add it
+ subboxes.insert(subboxes.end(), make_pair(subbox1, pos));
+ }
+ }
+
+ if (active0) ++iter0;
+ if (active1) ++iter1;
+ }
+ swap(old_subboxes, subboxes);
+ });
+ assert(old_subboxes.empty());
+ }
+
+
+
+ /** Memory usage */
+ template<typename T, int D>
+ size_t bboxset<T,D>::memory() const
+ {
+ size_t s = sizeof *this;
+ for (const auto& pos_subset: subsets) {
+ s += sizeof pos_subset;
+ const auto* const subsetp = pos_subset.second.get();
+ s += memoryof(*subsetp);
+ }
+ return s;
+ }
+
+ /** Input */
+ template<typename T, int D>
+ istream& bboxset<T,D>::input(istream& is)
+ {
+ T Tdummy;
+ try {
+ skipws(is);
+ consume(is, "bboxset<");
+ consume(is, typestring(Tdummy));
+ consume(is, ",");
+ int D_;
+ is >> D_;
+ if (D_ != D) {
+ ostringstream msg;
+ msg << "Input error: Wrong bboxset dimension " << D_ << ", expected " << D;
+ CCTK_WARN(CCTK_WARN_ALERT, msg.str().c_str());
+ throw input_error();
+ }
+ consume(is, ">(");
+ consume(is, "set<bbox>:");
+ set<bbox> bs;
+ is >> bs;
+ *this = bboxset(bs);
+ consume(is, ",");
+ consume(is, "stride:");
+ is >> stride;
+ consume(is, ",");
+ consume(is, "offset:");
+ is >> offset;
+ consume(is, ")");
+ } catch (input_error& err) {
+ ostringstream msg;
+ msg << "Input error while reading a bboxset<" << typestring(Tdummy) << ",0>";
+ CCTK_WARN(CCTK_WARN_ALERT, msg.str().c_str());
+ throw err;
+ }
+ is_poison_ = false;
+ return is;
+ }
+
+ /** Output */
+ template<typename T, int D>
+ ostream& bboxset<T,D>::output(ostream& os) const
+ {
+ assert(not is_poison());
+ T Tdummy;
+ set<bbox> bs;
+ serialise(bs);
+ return os << "bboxset<" << typestring(Tdummy) << "," << D << ">("
+ << "set<bbox>:" << bs << ","
+ << "stride:" << stride << ","
+ << "offset:" << offset << ")";
+ }
+
+
+
+ //////////////////////////////////////////////////////////////////////////////
+ // bboxset<T,0>
+ //////////////////////////////////////////////////////////////////////////////
+
+ template<typename T>
+ istream& bboxset<T,0>::input(istream& is)
+ {
+ T Tdummy;
+ try {
+ skipws(is);
+ consume(is, "bboxset<");
+ consume(is, typestring(Tdummy));
+ consume(is, ",0>(");
+ consume(is, "state:");
+ is >> state;
+ consume(is, ",");
+ consume(is, "stride:");
+ is >> stride;
+ consume(is, ",");
+ consume(is, "offset:");
+ is >> offset;
+ consume(is, ")");
+ } catch (input_error& err) {
+ ostringstream msg;
+ msg << "Input error while reading a bboxset<" << typestring(Tdummy) << ",0>";
+ CCTK_WARN(CCTK_WARN_ALERT, msg.str().c_str());
+ throw err;
+ }
+ is_poison_ = false;
+ return is;
+ }
+
+
+
+ template<typename T>
+ ostream& bboxset<T,0>::output(ostream& os) const
+ {
+ assert(not is_poison());
+ T Tdummy;
+ return os << "bboxset<" << typestring(Tdummy) << ",0>("
+ << "state:" << state << ","
+ << "stride:" << stride << ","
+ << "offset:" << offset << ")";
+ }
+
+
+
+} // namespace bboxset2
+
+#endif // #ifndef CARPET_NO_BBOXSET2
+
+
+
+#endif // #ifndef BBOXSET2_HH
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-06 13:20:36 +0000