diff options
Diffstat (limited to 'geom_bottleneck/include/basic_defs_bt.h')
| -rw-r--r-- | geom_bottleneck/include/basic_defs_bt.h | 749 |
1 files changed, 399 insertions, 350 deletions
diff --git a/geom_bottleneck/include/basic_defs_bt.h b/geom_bottleneck/include/basic_defs_bt.h index 69dc709..954696e 100644 --- a/geom_bottleneck/include/basic_defs_bt.h +++ b/geom_bottleneck/include/basic_defs_bt.h @@ -33,6 +33,7 @@ derivative works thereof, in binary and source code form. #include <ciso646> #endif +#include <utility> #include <vector> #include <stdexcept> #include <math.h> @@ -45,432 +46,480 @@ derivative works thereof, in binary and source code form. #include "def_debug_bt.h" #ifndef FOR_R_TDA -#include <iostream> -#endif -namespace hera { -namespace bt { +#include <iostream> -typedef int IdType; -constexpr IdType MinValidId = 10; +#endif -template<class Real = double> -struct Point { - Real x, y; +namespace hera { - bool operator==(const Point<Real> &other) const + template<class Real = double> + Real get_infinity() { - return ((this->x == other.x) and (this->y == other.y)); + return Real(-1.0); } - bool operator!=(const Point<Real> &other) const - { - return !(*this == other); - } + namespace bt { + - Point(Real ax, Real ay) : x(ax), y(ay) {} + typedef int IdType; + constexpr IdType MinValidId = 10; - Point() : x(0.0), y(0.0) {} + template<class Real = double> + struct Point + { + Real x, y; + + bool operator==(const Point<Real>& other) const + { + return ((this->x == other.x) and (this->y == other.y)); + } + + bool operator!=(const Point<Real>& other) const + { + return !(*this == other); + } + + Point(Real ax, Real ay) : + x(ax), y(ay) + {} + + Point() : + x(0.0), y(0.0) + {} #ifndef FOR_R_TDA - template<class R> - friend std::ostream& operator<<(std::ostream& output, const Point<R>& p) - { - output << "(" << p.x << ", " << p.y << ")"; - return output; - } + template<class R> + friend std::ostream& operator<<(std::ostream& output, const Point<R>& p) + { + output << "(" << p.x << ", " << p.y << ")"; + return output; + } #endif -}; - -template<class Real = double> -struct DiagramPoint { - // Points above the diagonal have type NORMAL - // Projections onto the diagonal have type DIAG - // for DIAG points only x-coordinate is relevant - // to-do: add getters/setters, checks in constructors, etc - enum Type { - NORMAL, DIAG - }; - // data members -private: - Real x, y; -public: - Type type; - IdType id; - - // operators, constructors - bool operator==(const DiagramPoint<Real> &other) const - { - // compare by id only - assert(this->id >= MinValidId); - assert(other.id >= MinValidId); - bool areEqual{ this->id == other.id }; - assert(!areEqual or ((this->x == other.x) and (this->y == other.y) and (this->type == other.type))); - return areEqual; - } + }; + + template<class Real = double> + struct DiagramPoint + { + // Points above the diagonal have type NORMAL + // Projections onto the diagonal have type DIAG + // for DIAG points only x-coordinate is relevant + // to-do: add getters/setters, checks in constructors, etc + enum Type + { + NORMAL, DIAG + }; + // data members + private: + Real x, y; + public: + Type type; + IdType id; + + // operators, constructors + bool operator==(const DiagramPoint<Real>& other) const + { + // compare by id only + assert(this->id >= MinValidId); + assert(other.id >= MinValidId); + bool areEqual { this->id == other.id }; + assert(!areEqual or ((this->x == other.x) and (this->y == other.y) and (this->type == other.type))); + return areEqual; + } - bool operator!=(const DiagramPoint &other) const - { - return !(*this == other); - } + bool operator!=(const DiagramPoint& other) const + { + return !(*this == other); + } - DiagramPoint(Real _x, Real _y, Type _type, IdType _id) : - x(_x), - y(_y), - type(_type), - id(_id) - { - if ( _y == _x and _type != DIAG) - throw std::runtime_error("Point on the main diagonal must have DIAG type"); + DiagramPoint() : + x(0.0), + y(0.0), + type(DiagramPoint::DIAG), + id(MinValidId - 1) + { + } - } + DiagramPoint(Real _x, Real _y, Type _type, IdType _id) : + x(_x), + y(_y), + type(_type), + id(_id) + { + if (_y == _x and _type != DIAG) { + throw std::runtime_error("Point on the main diagonal must have DIAG type"); + } + } - bool isDiagonal() const { return type == DIAG; } - bool isNormal() const { return type == NORMAL; } + bool isDiagonal() const + { return type == DIAG; } - bool isInfinity() const - { - return x == std::numeric_limits<Real>::infinity() or - x == -std::numeric_limits<Real>::infinity() or - y == std::numeric_limits<Real>::infinity() or - y == -std::numeric_limits<Real>::infinity(); - } + bool isNormal() const + { return type == NORMAL; } - Real inline getRealX() const // return the x-coord - { - return x; - } + bool isInfinity() const + { + return x == std::numeric_limits<Real>::infinity() or + x == -std::numeric_limits<Real>::infinity() or + y == std::numeric_limits<Real>::infinity() or + y == -std::numeric_limits<Real>::infinity(); + } - Real inline getRealY() const // return the y-coord - { - return y; - } + Real inline getRealX() const // return the x-coord + { + return x; + } + + Real inline getRealY() const // return the y-coord + { + return y; + } #ifndef FOR_R_TDA - template<class R> - friend std::ostream& operator<<(std::ostream& output, const DiagramPoint<R>& p) - { - if ( p.isDiagonal() ) { - output << "(" << p.x << ", " << p.y << ", " << 0.5 * (p.x + p.y) << ", " << p.id << " DIAG )"; - } else { - output << "(" << p.x << ", " << p.y << ", " << p.id << " NORMAL)"; - } - return output; - } -#endif -}; + template<class R> + friend std::ostream& operator<<(std::ostream& output, const DiagramPoint<R>& p) + { + if (p.isDiagonal()) { + output << "(" << p.x << ", " << p.y << ", " << 0.5 * (p.x + p.y) << ", " << p.id << " DIAG )"; + } else { + output << "(" << p.x << ", " << p.y << ", " << p.id << " NORMAL)"; + } + return output; + } -// compute l-inf distance between two diagram points -template<class Real> -Real distLInf(const DiagramPoint<Real>& a, const DiagramPoint<Real>& b) -{ - if ( a.isDiagonal() and b.isDiagonal() ) { - // distance between points on the diagonal is 0 - return 0.0; - } - // otherwise distance is a usual l-inf distance - return std::max(fabs(a.getRealX() - b.getRealX()), fabs(a.getRealY() - b.getRealY())); -} +#endif + }; -template <class T> -inline void hash_combine(std::size_t & seed, const T & v) -{ - std::hash<T> hasher; - seed ^= hasher(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2); -} + template<class Real> + using MatchingEdge = std::pair<DiagramPoint<Real>, DiagramPoint<Real>>; -template<class Real = double> -struct DiagramPointHash { - size_t operator()(const DiagramPoint<Real> &p) const - { - assert(p.id >= MinValidId); - return std::hash<int>()(p.id); - } -}; + // compute l-inf distance between two diagram points + template<class Real> + Real distLInf(const DiagramPoint<Real>& a, const DiagramPoint<Real>& b) + { + if (a.isDiagonal() and b.isDiagonal()) { + // distance between points on the diagonal is 0 + return 0.0; + } + // otherwise distance is a usual l-inf distance + return std::max(fabs(a.getRealX() - b.getRealX()), fabs(a.getRealY() - b.getRealY())); + } -template<class Real = double> -Real distLInf(const DiagramPoint<Real> &a, const DiagramPoint<Real> &b); -//template<class Real = double> -//typedef std::unordered_set<Point, PointHash> PointSet; -template<class Real_ = double> -class DiagramPointSet; + template<class Real = double> + Real get_infinity() + { + return Real(-1.0); + } -template<class Real> -void addProjections(DiagramPointSet<Real> &A, DiagramPointSet<Real> &B); + template<class T> + inline void hash_combine(std::size_t& seed, const T& v) + { + std::hash<T> hasher; + seed ^= hasher(v) + 0x9e3779b9 + (seed << 6) + (seed >> 2); + } -template<class Real_> -class DiagramPointSet { -public: + template<class Real = double> + struct DiagramPointHash + { + size_t operator()(const DiagramPoint<Real>& p) const + { + assert(p.id >= MinValidId); + return std::hash<int>()(p.id); + } + }; - using Real = Real_; - using DgmPoint = DiagramPoint<Real>; - using DgmPointHash = DiagramPointHash<Real>; - using const_iterator = typename std::unordered_set<DgmPoint, DgmPointHash>::const_iterator; - using iterator = typename std::unordered_set<DgmPoint, DgmPointHash>::iterator; + template<class Real = double> + Real distLInf(const DiagramPoint<Real>& a, const DiagramPoint<Real>& b); -private: + //template<class Real = double> + //typedef std::unordered_set<Point, PointHash> PointSet; + template<class Real_ = double> + class DiagramPointSet; - bool isLinked { false }; - IdType maxId { MinValidId + 1 }; - std::unordered_set<DgmPoint, DgmPointHash> points; + template<class Real> + void addProjections(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B); -public: + template<class Real_> + class DiagramPointSet + { + public: - void insert(const DgmPoint& p) - { - points.insert(p); - if (p.id > maxId) { - maxId = p.id + 1; - } - } + using Real = Real_; + using DgmPoint = DiagramPoint<Real>; + using DgmPointHash = DiagramPointHash<Real>; + using const_iterator = typename std::unordered_set<DgmPoint, DgmPointHash>::const_iterator; + using iterator = typename std::unordered_set<DgmPoint, DgmPointHash>::iterator; - void erase(const DgmPoint& p, bool doCheck = true) - { - // if doCheck, erasing non-existing elements causes assert - auto it = points.find(p); - if (it != points.end()) { - points.erase(it); - } else { - assert(!doCheck); - } - } + private: + bool isLinked { false }; + IdType maxId { MinValidId + 1 }; + std::unordered_set<DgmPoint, DgmPointHash> points; - void erase(const const_iterator it) - { - points.erase(it); - } + public: - void removeDiagonalPoints() - { - if (isLinked) { - auto ptIter = points.begin(); - while(ptIter != points.end()) { - if (ptIter->isDiagonal()) { - ptIter = points.erase(ptIter); + void insert(const DgmPoint& p) + { + points.insert(p); + if (p.id > maxId) { + maxId = p.id + 1; + } + } + + void erase(const DgmPoint& p, bool doCheck = true) + { + // if doCheck, erasing non-existing elements causes assert + auto it = points.find(p); + if (it != points.end()) { + points.erase(it); } else { - ptIter++; + assert(!doCheck); } } - isLinked = false; - } - } - size_t size() const - { - return points.size(); - } - void reserve(const size_t newSize) - { - points.reserve(newSize); - } + void erase(const const_iterator it) + { + points.erase(it); + } - void clear() - { - points.clear(); - } + void removeDiagonalPoints() + { + if (isLinked) { + auto ptIter = points.begin(); + while (ptIter != points.end()) { + if (ptIter->isDiagonal()) { + ptIter = points.erase(ptIter); + } else { + ptIter++; + } + } + isLinked = false; + } + } - bool empty() const - { - return points.empty(); - } + size_t size() const + { + return points.size(); + } - bool hasElement(const DgmPoint &p) const - { - return points.find(p) != points.end(); - } + void reserve(const size_t newSize) + { + points.reserve(newSize); + } - iterator find(const DgmPoint &p) - { - return points.find(p); - } + void clear() + { + points.clear(); + } - iterator begin() - { - return points.begin(); - } + bool empty() const + { + return points.empty(); + } - iterator end() - { - return points.end(); - } + bool hasElement(const DgmPoint& p) const + { + return points.find(p) != points.end(); + } - const_iterator cbegin() const - { - return points.cbegin(); - } + iterator find(const DgmPoint& p) + { + return points.find(p); + } - const_iterator cend() const - { - return points.cend(); - } + iterator begin() + { + return points.begin(); + } + iterator end() + { + return points.end(); + } - const_iterator find(const DgmPoint &p) const - { - return points.find(p); - } + const_iterator cbegin() const + { + return points.cbegin(); + } + + const_iterator cend() const + { + return points.cend(); + } + + + const_iterator find(const DgmPoint& p) const + { + return points.find(p); + } #ifndef FOR_R_TDA - template<class R> - friend std::ostream& operator<<(std::ostream& output, const DiagramPointSet<R>& ps) - { - output << "{ "; - for(auto pit = ps.cbegin(); pit != ps.cend(); ++pit) { - output << *pit << ", "; - } - output << "\b\b }"; - return output; - } + + template<class R> + friend std::ostream& operator<<(std::ostream& output, const DiagramPointSet<R>& ps) + { + output << "{ "; + for (auto pit = ps.cbegin(); pit != ps.cend(); ++pit) { + output << *pit << ", "; + } + output << "\b\b }"; + return output; + } + #endif - friend void addProjections<Real_>(DiagramPointSet<Real_>& A, DiagramPointSet<Real_>& B); + friend void addProjections<Real_>(DiagramPointSet<Real_>& A, DiagramPointSet<Real_>& B); - template<class PairIterator> - void fillIn(PairIterator begin_iter, PairIterator end_iter) - { - isLinked = false; - clear(); - IdType uniqueId = MinValidId + 1; - for (auto iter = begin_iter; iter != end_iter; ++iter) { - insert(DgmPoint(iter->first, iter->second, DgmPoint::NORMAL, uniqueId++)); - } - } + template<class PairIterator> + void fillIn(PairIterator begin_iter, PairIterator end_iter) + { + isLinked = false; + clear(); + IdType uniqueId = MinValidId + 1; + for (auto iter = begin_iter; iter != end_iter; ++iter) { + insert(DgmPoint(iter->first, iter->second, DgmPoint::NORMAL, uniqueId++)); + } + } - template<class PointContainer> - void fillIn(const PointContainer& dgm_cont) - { - using Traits = DiagramTraits<PointContainer>; - isLinked = false; - clear(); - IdType uniqueId = MinValidId + 1; - for (const auto& pt : dgm_cont) { - Real x = Traits::get_x(pt); - Real y = Traits::get_y(pt); - insert(DgmPoint(x, y, DgmPoint::NORMAL, uniqueId++)); - } - } + template<class PointContainer> + void fillIn(const PointContainer& dgm_cont) + { + using Traits = DiagramTraits<PointContainer>; + isLinked = false; + clear(); + IdType uniqueId = MinValidId + 1; + for (const auto& pt : dgm_cont) { + Real x = Traits::get_x(pt); + Real y = Traits::get_y(pt); + insert(DgmPoint(x, y, DgmPoint::NORMAL, uniqueId++)); + } + } - // ctor from range - template<class PairIterator> - DiagramPointSet(PairIterator begin_iter, PairIterator end_iter) - { - fillIn(begin_iter, end_iter); - } + // ctor from range + template<class PairIterator> + DiagramPointSet(PairIterator begin_iter, PairIterator end_iter) + { + fillIn(begin_iter, end_iter); + } - // ctor from container, uses DiagramTraits - template<class PointContainer> - DiagramPointSet(const PointContainer& dgm) - { - fillIn(dgm); - } + // ctor from container, uses DiagramTraits + template<class PointContainer> + DiagramPointSet(const PointContainer& dgm) + { + fillIn(dgm); + } - // default ctor, empty diagram - DiagramPointSet(IdType minId = MinValidId + 1) : maxId(minId + 1) {}; + // default ctor, empty diagram + DiagramPointSet(IdType minId = MinValidId + 1) : + maxId(minId + 1) + {}; - IdType nextId() { return maxId + 1; } + IdType nextId() + { return maxId + 1; } -}; // DiagramPointSet + }; // DiagramPointSet -template<class Real, class DiagPointContainer> -Real getFurthestDistance3Approx(DiagPointContainer& A, DiagPointContainer& B) { - Real result{0.0}; - DiagramPoint<Real> begA = *(A.begin()); - DiagramPoint<Real> optB = *(B.begin()); - for (const auto &pointB : B) { - if (distLInf(begA, pointB) > result) { - result = distLInf(begA, pointB); - optB = pointB; - } - } - for (const auto &pointA : A) { - if (distLInf(pointA, optB) > result) { - result = distLInf(pointA, optB); - } - } - return result; -} - -// preprocess diagrams A and B by adding projections onto diagonal of points of -// A to B and vice versa. Also removes points at infinity! -// NB: ids of points will be changed! -template<class Real_> -void addProjections(DiagramPointSet<Real_>& A, DiagramPointSet<Real_>& B) -{ - - using Real = Real_; - using DgmPoint = DiagramPoint<Real>; - using DgmPointSet = DiagramPointSet<Real>; - - IdType uniqueId {MinValidId + 1}; - DgmPointSet newA, newB; - - // copy normal points from A to newA - // add projections to newB - for(auto& pA : A) { - if (pA.isNormal() and not pA.isInfinity()) { - // add pA's projection to B - DgmPoint dpA {pA.getRealX(), pA.getRealY(), DgmPoint::NORMAL, uniqueId++}; - DgmPoint dpB {(pA.getRealX() +pA.getRealY())/2, (pA.getRealX() +pA.getRealY())/2, DgmPoint::DIAG, uniqueId++}; - newA.insert(dpA); - newB.insert(dpB); + template<class Real, class DiagPointContainer> + Real getFurthestDistance3Approx(DiagPointContainer& A, DiagPointContainer& B) + { + Real result { 0.0 }; + DiagramPoint<Real> begA = *(A.begin()); + DiagramPoint<Real> optB = *(B.begin()); + for (const auto& pointB : B) { + if (distLInf(begA, pointB) > result) { + result = distLInf(begA, pointB); + optB = pointB; + } + } + for (const auto& pointA : A) { + if (distLInf(pointA, optB) > result) { + result = distLInf(pointA, optB); + } + } + return result; } - } - for(auto& pB : B) { - if (pB.isNormal() and not pB.isInfinity()) { - // add pB's projection to A - DgmPoint dpB {pB.getRealX(), pB.getRealY(), DgmPoint::NORMAL, uniqueId++}; - DgmPoint dpA {(pB.getRealX() +pB.getRealY())/2, (pB.getRealX() +pB.getRealY())/2, DgmPoint::DIAG, uniqueId++}; - newB.insert(dpB); - newA.insert(dpA); + // preprocess diagrams A and B by adding projections onto diagonal of points of + // A to B and vice versa. Also removes points at infinity! + // NB: ids of points will be changed! + template<class Real_> + void addProjections(DiagramPointSet<Real_>& A, DiagramPointSet<Real_>& B) + { + + using Real = Real_; + using DgmPoint = DiagramPoint<Real>; + using DgmPointSet = DiagramPointSet<Real>; + + IdType uniqueId { MinValidId + 1 }; + DgmPointSet newA, newB; + + // copy normal points from A to newA + // add projections to newB + for (auto& pA : A) { + if (pA.isNormal() and not pA.isInfinity()) { + // add pA's projection to B + DgmPoint dpA { pA.getRealX(), pA.getRealY(), DgmPoint::NORMAL, uniqueId++ }; + DgmPoint dpB { (pA.getRealX() + pA.getRealY()) / 2, (pA.getRealX() + pA.getRealY()) / 2, + DgmPoint::DIAG, uniqueId++ }; + newA.insert(dpA); + newB.insert(dpB); + } + } + + for (auto& pB : B) { + if (pB.isNormal() and not pB.isInfinity()) { + // add pB's projection to A + DgmPoint dpB { pB.getRealX(), pB.getRealY(), DgmPoint::NORMAL, uniqueId++ }; + DgmPoint dpA { (pB.getRealX() + pB.getRealY()) / 2, (pB.getRealX() + pB.getRealY()) / 2, + DgmPoint::DIAG, uniqueId++ }; + newB.insert(dpB); + newA.insert(dpA); + } + } + + A = newA; + B = newB; + A.isLinked = true; + B.isLinked = true; } - } - A = newA; - B = newB; - A.isLinked = true; - B.isLinked = true; -} - - -//#ifndef FOR_R_TDA - -//template<class Real> -//std::ostream& operator<<(std::ostream& output, const DiagramPoint<Real>& p) -//{ -// if ( p.isDiagonal() ) { -// output << "(" << p.x << ", " << p.y << ", " << 0.5 * (p.x + p.y) << ", " << p.id << " DIAG )"; -// } else { -// output << "(" << p.x << ", " << p.y << ", " << p.id << " NORMAL)"; -// } -// return output; -//} - -//template<class Real> -//std::ostream& operator<<(std::ostream& output, const DiagramPointSet<Real>& ps) -//{ -// output << "{ "; -// for(auto pit = ps.cbegin(); pit != ps.cend(); ++pit) { -// output << *pit << ", "; -// } -// output << "\b\b }"; -// return output; -//} -//#endif // FOR_R_TDA - - -} // end namespace bt + //#ifndef FOR_R_TDA + + //template<class Real> + //std::ostream& operator<<(std::ostream& output, const DiagramPoint<Real>& p) + //{ + // if ( p.isDiagonal() ) { + // output << "(" << p.x << ", " << p.y << ", " << 0.5 * (p.x + p.y) << ", " << p.id << " DIAG )"; + // } else { + // output << "(" << p.x << ", " << p.y << ", " << p.id << " NORMAL)"; + // } + // return output; + //} + + //template<class Real> + //std::ostream& operator<<(std::ostream& output, const DiagramPointSet<Real>& ps) + //{ + // output << "{ "; + // for(auto pit = ps.cbegin(); pit != ps.cend(); ++pit) { + // output << *pit << ", "; + // } + // output << "\b\b }"; + // return output; + //} + //#endif // FOR_R_TDA + + + } // end namespace bt } // end namespace hera #endif // HERA_BASIC_DEFS_BT_H |