summaryrefslogtreecommitdiff
path: root/geom_bottleneck/include
diff options
context:
space:
mode:
Diffstat (limited to 'geom_bottleneck/include')
-rw-r--r--geom_bottleneck/include/basic_defs_bt.h475
-rw-r--r--geom_bottleneck/include/bottleneck.h116
-rw-r--r--geom_bottleneck/include/bottleneck_detail.h85
-rw-r--r--geom_bottleneck/include/bottleneck_detail.hpp785
-rw-r--r--geom_bottleneck/include/bound_match.h107
-rw-r--r--geom_bottleneck/include/bound_match.hpp473
-rw-r--r--geom_bottleneck/include/def_debug_bt.h42
-rw-r--r--geom_bottleneck/include/diagram_reader.h139
-rw-r--r--geom_bottleneck/include/diagram_traits.h45
-rw-r--r--geom_bottleneck/include/dnn/geometry/euclidean-fixed.h162
-rw-r--r--geom_bottleneck/include/dnn/local/kd-tree.h106
-rw-r--r--geom_bottleneck/include/dnn/local/kd-tree.hpp296
-rw-r--r--geom_bottleneck/include/dnn/local/search-functors.h119
-rw-r--r--geom_bottleneck/include/dnn/parallel/tbb.h235
-rw-r--r--geom_bottleneck/include/dnn/parallel/utils.h100
-rw-r--r--geom_bottleneck/include/dnn/utils.h47
-rw-r--r--geom_bottleneck/include/neighb_oracle.h295
17 files changed, 3627 insertions, 0 deletions
diff --git a/geom_bottleneck/include/basic_defs_bt.h b/geom_bottleneck/include/basic_defs_bt.h
new file mode 100644
index 0000000..ad09986
--- /dev/null
+++ b/geom_bottleneck/include/basic_defs_bt.h
@@ -0,0 +1,475 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef HERA_BASIC_DEFS_BT_H
+#define HERA_BASIC_DEFS_BT_H
+
+#ifdef _WIN32
+#include <ciso646>
+#endif
+
+#include <vector>
+#include <stdexcept>
+#include <math.h>
+#include <cstddef>
+#include <unordered_map>
+#include <unordered_set>
+#include <string>
+#include <assert.h>
+
+#include "def_debug_bt.h"
+
+#ifndef FOR_R_TDA
+#include <iostream>
+#endif
+
+namespace hera {
+namespace bt {
+
+typedef int IdType;
+constexpr IdType MinValidId = 10;
+
+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;
+ }
+
+#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;
+ }
+
+ 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");
+
+ }
+
+
+ bool isDiagonal(void) const { return type == DIAG; }
+
+ bool isNormal(void) const { return type == NORMAL; }
+
+ 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
+
+};
+
+// 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 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 = double>
+//struct PointHash {
+// size_t operator()(const Point<Real> &p) const
+// {
+// size_t seed = 0;
+// hash_combine(seed, p.x);
+// hash_combine(seed, p.y);
+// return seed;
+// }
+//};
+
+template<class Real = double>
+struct DiagramPointHash {
+ size_t operator()(const DiagramPoint<Real> &p) const
+ {
+ assert(p.id >= MinValidId);
+ return std::hash<int>()(p.id);
+ }
+};
+
+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>
+void addProjections(DiagramPointSet<Real> &A, DiagramPointSet<Real> &B);
+
+template<class Real_>
+class DiagramPointSet {
+public:
+
+ 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;
+
+private:
+
+ bool isLinked { false };
+ IdType maxId { MinValidId + 1 };
+ std::unordered_set<DgmPoint, DgmPointHash> points;
+
+public:
+
+ 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 {
+ assert(!doCheck);
+ }
+ }
+
+
+ void erase(const const_iterator it)
+ {
+ points.erase(it);
+ }
+
+ void removeDiagonalPoints()
+ {
+ if (isLinked) {
+ auto ptIter = points.begin();
+ while(ptIter != points.end()) {
+ if (ptIter->isDiagonal()) {
+ ptIter = points.erase(ptIter);
+ } else {
+ ptIter++;
+ }
+ }
+ isLinked = false;
+ }
+ }
+
+ size_t size() const
+ {
+ return points.size();
+ }
+
+ void reserve(const size_t newSize)
+ {
+ points.reserve(newSize);
+ }
+
+ void clear()
+ {
+ points.clear();
+ }
+
+ bool empty() const
+ {
+ return points.empty();
+ }
+
+ bool hasElement(const DgmPoint &p) const
+ {
+ return points.find(p) != points.end();
+ }
+
+ iterator find(const DgmPoint &p)
+ {
+ return points.find(p);
+ }
+
+ 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();
+ }
+
+#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;
+ }
+#endif
+
+ 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 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 container, uses DiagramTraits
+ template<class PointContainer>
+ DiagramPointSet(const PointContainer& dgm)
+ {
+ fillIn(dgm);
+ }
+
+
+ // default ctor, empty diagram
+ DiagramPointSet(IdType minId = MinValidId + 1) : maxId(minId + 1) {};
+
+ IdType nextId() { return maxId + 1; }
+
+}; // 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. 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()) {
+ 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()) {
+ 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;
+}
+
+
+//#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
diff --git a/geom_bottleneck/include/bottleneck.h b/geom_bottleneck/include/bottleneck.h
new file mode 100644
index 0000000..d0d82b6
--- /dev/null
+++ b/geom_bottleneck/include/bottleneck.h
@@ -0,0 +1,116 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef HERA_BOTTLENECK_H
+#define HERA_BOTTLENECK_H
+
+
+#include <fstream>
+#include <vector>
+#include <algorithm>
+#include <limits>
+#include <random>
+
+#include "diagram_traits.h"
+#include "diagram_reader.h"
+#include "bottleneck_detail.h"
+#include "basic_defs_bt.h"
+#include "bound_match.h"
+
+namespace hera {
+
+// functions taking containers as input
+// template parameter PairContainer must be a container of pairs of real
+// numbers (pair.first = x-coordinate, pair.second = y-coordinate)
+// PairContainer class must support iteration of the form
+// for(it = pairContainer.begin(); it != pairContainer.end(); ++it)
+
+
+// get exact bottleneck distance,
+template<class PairContainer>
+typename DiagramTraits<PairContainer>::RealType
+bottleneckDistExact(PairContainer& dgm_A, PairContainer& dgm_B)
+{
+ using Real = typename DiagramTraits<PairContainer>::RealType;
+ hera::bt::DiagramPointSet<Real> a(dgm_A);
+ hera::bt::DiagramPointSet<Real> b(dgm_B);
+ return hera::bt::bottleneckDistExact<Real>(a, b, 14);
+}
+
+// get exact bottleneck distance,
+template<class PairContainer>
+typename DiagramTraits<PairContainer>::RealType
+bottleneckDistExact(PairContainer& dgm_A, PairContainer& dgm_B, const int decPrecision)
+{
+ using Real = typename DiagramTraits<PairContainer>::RealType;
+ hera::bt::DiagramPointSet<Real> a(dgm_A);
+ hera::bt::DiagramPointSet<Real> b(dgm_B);
+ return hera::bt::bottleneckDistExact(a, b, decPrecision);
+}
+
+// return the interval (distMin, distMax) such that:
+// a) actual bottleneck distance between A and B is contained in the interval
+// b) if the interval is not (0,0), then (distMax - distMin) / distMin < delta
+template<class PairContainer>
+std::pair<typename DiagramTraits<PairContainer>::RealType, typename DiagramTraits<PairContainer>::RealType>
+bottleneckDistApproxInterval(PairContainer& dgm_A, PairContainer& dgm_B, const typename DiagramTraits<PairContainer>::RealType delta)
+{
+ using Real = typename DiagramTraits<PairContainer>::RealType;
+ hera::bt::DiagramPointSet<Real> a(dgm_A);
+ hera::bt::DiagramPointSet<Real> b(dgm_B);
+ return hera::bt::bottleneckDistApproxInterval(a, b, delta);
+}
+
+
+template<class PairContainer>
+typename DiagramTraits<PairContainer>::RealType
+bottleneckDistApproxHeur(PairContainer& dgm_A, PairContainer& dgm_B, const typename DiagramTraits<PairContainer>::RealType delta)
+{
+ using Real = typename DiagramTraits<PairContainer>::RealType;
+ hera::bt::DiagramPointSet<Real> a(dgm_A);
+ hera::bt::DiagramPointSet<Real> b(dgm_B);
+ std::pair<Real, Real> resPair = hera::bt::bottleneckDistApproxIntervalHeur(a, b, delta);
+ return resPair.second;
+}
+
+
+// get approximate distance,
+// see bottleneckDistApproxInterval
+template<class PairContainer>
+typename DiagramTraits<PairContainer>::RealType
+bottleneckDistApprox(PairContainer& A, PairContainer& B, const typename DiagramTraits<PairContainer>::RealType delta)
+{
+ using Real = typename DiagramTraits<PairContainer>::RealType;
+ hera::bt::DiagramPointSet<Real> a(A.begin(), A.end());
+ hera::bt::DiagramPointSet<Real> b(B.begin(), B.end());
+ return hera::bt::bottleneckDistApprox(a, b, delta);
+}
+
+} // end namespace hera
+
+#endif
diff --git a/geom_bottleneck/include/bottleneck_detail.h b/geom_bottleneck/include/bottleneck_detail.h
new file mode 100644
index 0000000..27c3c5d
--- /dev/null
+++ b/geom_bottleneck/include/bottleneck_detail.h
@@ -0,0 +1,85 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef HERA_BOTTLENECK_DETAIL_H
+#define HERA_BOTTLENECK_DETAIL_H
+
+
+#include <fstream>
+#include <vector>
+#include <algorithm>
+#include <limits>
+#include <random>
+
+#include "diagram_traits.h"
+#include "basic_defs_bt.h"
+#include "bound_match.h"
+
+namespace hera {
+
+
+namespace bt {
+
+
+
+// functions taking DiagramPointSet as input.
+// ATTENTION: parameters A and B (diagrams) will be changed after the call
+// (projections added).
+
+// return the interval (distMin, distMax) such that:
+// a) actual bottleneck distance between A and B is contained in the interval
+// b) if the interval is not (0,0), then (distMax - distMin) / distMin < epsilon
+template<class Real>
+std::pair<Real, Real> bottleneckDistApproxInterval(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const Real epsilon);
+
+
+// heuristic (sample diagram to estimate the distance)
+template<class Real>
+std::pair<Real, Real> bottleneckDistApproxIntervalHeur(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const Real epsilon);
+
+// get approximate distance,
+// see bottleneckDistApproxInterval
+template<class Real>
+Real bottleneckDistApprox(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const Real epsilon);
+
+// get exact bottleneck distance,
+template<class Real>
+Real bottleneckDistExact(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const int decPrecision);
+
+// get exact bottleneck distance,
+template<class Real>
+Real bottleneckDistExact(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B);
+
+} // end namespace bt
+
+
+} // end namespace hera
+
+#include "bottleneck_detail.hpp"
+
+#endif // HERA_BOTTLENECK_DETAIL_H
diff --git a/geom_bottleneck/include/bottleneck_detail.hpp b/geom_bottleneck/include/bottleneck_detail.hpp
new file mode 100644
index 0000000..64c6696
--- /dev/null
+++ b/geom_bottleneck/include/bottleneck_detail.hpp
@@ -0,0 +1,785 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef HERA_BOTTLENECK_HPP
+#define HERA_BOTTLENECK_HPP
+
+#ifdef FOR_R_TDA
+#undef DEBUG_BOUND_MATCH
+#undef DEBUG_MATCHING
+#undef VERBOSE_BOTTLENECK
+#endif
+
+
+#include <iomanip>
+#include <sstream>
+#include <string>
+#include <cctype>
+
+#include "bottleneck_detail.h"
+
+namespace hera {
+namespace bt {
+
+// return the interval (distMin, distMax) such that:
+// a) actual bottleneck distance between A and B is contained in the interval
+// b) if the interval is not (0,0), then (distMax - distMin) / distMin < epsilon
+template<class Real>
+std::pair<Real, Real> bottleneckDistApproxInterval(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const Real epsilon)
+{
+ // empty diagrams are not considered as error
+ if (A.empty() and B.empty())
+ return std::make_pair(0.0, 0.0);
+
+ // link diagrams A and B by adding projections
+ addProjections(A, B);
+
+ // TODO: think about that!
+ // we need one threshold for checking if the distance is 0,
+ // another one for the oracle!
+ constexpr Real epsThreshold { 1.0e-10 };
+ std::pair<Real, Real> result { 0.0, 0.0 };
+ bool useRangeSearch { true };
+ // construct an oracle
+ BoundMatchOracle<Real> oracle(A, B, epsThreshold, useRangeSearch);
+ // check for distance = 0
+ if (oracle.isMatchLess(2*epsThreshold)) {
+ return result;
+ }
+ // get a 3-approximation of maximal distance between A and B
+ // as a starting value for probe distance
+ Real distProbe { getFurthestDistance3Approx<Real, DiagramPointSet<Real>>(A, B) };
+ // aliases for result components
+ Real& distMin {result.first};
+ Real& distMax {result.second};
+
+ if ( oracle.isMatchLess(distProbe) ) {
+ // distProbe is an upper bound,
+ // find lower bound with binary search
+ do {
+ distMax = distProbe;
+ distProbe /= 2.0;
+ } while (oracle.isMatchLess(distProbe));
+ distMin = distProbe;
+ } else {
+ // distProbe is a lower bound,
+ // find upper bound with exponential search
+ do {
+ distMin = distProbe;
+ distProbe *= 2.0;
+ } while (!oracle.isMatchLess(distProbe));
+ distMax = distProbe;
+ }
+ // bounds are found, perform binary search
+ //std::cout << "Bounds found, distMin = " << distMin << ", distMax = " << distMax << ", ratio = " << ( distMax - distMin ) / distMin << std::endl ;
+ distProbe = ( distMin + distMax ) / 2.0;
+ while ( ( distMax - distMin ) / distMin >= epsilon ) {
+ if (oracle.isMatchLess(distProbe)) {
+ distMax = distProbe;
+ } else {
+ distMin = distProbe;
+ }
+ distProbe = ( distMin + distMax ) / 2.0;
+ }
+ return result;
+}
+
+template<class Real>
+void sampleDiagramForHeur(const DiagramPointSet<Real>& dgmIn, DiagramPointSet<Real>& dgmOut)
+{
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "Entered sampleDiagramForHeur, dgmIn.size = " << dgmIn.size() << std::endl;
+#endif
+ struct pair_hash {
+ std::size_t operator()(const std::pair<Real, Real> p) const
+ {
+ return std::hash<Real>()(p.first) ^ std::hash<Real>()(p.second);
+ }
+ };
+ std::unordered_map<std::pair<Real, Real>, int, pair_hash> m;
+ for(auto ptIter = dgmIn.cbegin(); ptIter != dgmIn.cend(); ++ptIter) {
+ if (ptIter->isNormal()) {
+ m[std::make_pair(ptIter->getRealX(), ptIter->getRealY())]++;
+ }
+ }
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "map filled in, m.size = " << m.size() << std::endl;
+#endif
+ if (m.size() < 2) {
+ dgmOut = dgmIn;
+ return;
+ }
+ std::vector<int> v;
+ for(const auto& ptQtyPair : m) {
+ v.push_back(ptQtyPair.second);
+ }
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "v filled in, v.size = " << v.size() << std::endl;
+#endif
+ std::sort(v.begin(), v.end());
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "v sorted" << std::endl;
+#endif
+ int maxLeap = v[1] - v[0];
+ int cutVal = v[0];
+ for(int i = 1; i < static_cast<int>(v.size())- 1; ++i) {
+ int currLeap = v[i+1] - v[i];
+ if (currLeap > maxLeap) {
+ maxLeap = currLeap;
+ cutVal = v[i];
+ }
+ }
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "cutVal found, cutVal = " << cutVal << std::endl;
+#endif
+ std::vector<std::pair<Real, Real>> vv;
+ // keep points whose multiplicites are at most cutVal
+ // quick-and-dirty: fill in vv with copies of each point
+ // to construct DiagramPointSet from it later
+ for(const auto& ptQty : m) {
+ if (ptQty.second < cutVal) {
+ for(int i = 0; i < ptQty.second; ++i) {
+ vv.push_back(std::make_pair(ptQty.first.first, ptQty.first.second));
+ }
+ }
+ }
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "vv filled in, vv.size = " << v.size() << std::endl;
+#endif
+ dgmOut.clear();
+ dgmOut = DiagramPointSet<Real>(vv.begin(), vv.end());
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "dgmOut filled in, dgmOut.size = " << dgmOut.size() << std::endl;
+#endif
+}
+
+
+// return the interval (distMin, distMax) such that:
+// a) actual bottleneck distance between A and B is contained in the interval
+// b) if the interval is not (0,0), then (distMax - distMin) / distMin < epsilon
+template<class Real>
+std::pair<Real, Real> bottleneckDistApproxIntervalWithInitial(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const Real epsilon, const std::pair<Real, Real> initialGuess)
+{
+ // empty diagrams are not considered as error
+ if (A.empty() and B.empty())
+ return std::make_pair(0.0, 0.0);
+
+ // link diagrams A and B by adding projections
+ addProjections(A, B);
+
+ constexpr Real epsThreshold { 1.0e-10 };
+ std::pair<Real, Real> result { 0.0, 0.0 };
+ bool useRangeSearch { true };
+ // construct an oracle
+ BoundMatchOracle<Real> oracle(A, B, epsThreshold, useRangeSearch);
+ Real& distMin {result.first};
+ Real& distMax {result.second};
+
+ // initialize search interval from initialGuess
+ distMin = initialGuess.first;
+ distMax = initialGuess.second;
+
+ assert(distMin <= distMax);
+
+ // make sure that distMin is a lower bound
+ while(oracle.isMatchLess(distMin)) {
+ // distMin is in fact an upper bound, so assign it to distMax
+ distMax = distMin;
+ // and decrease distMin by 5 %
+ distMin = 0.95 * distMin;
+ }
+
+ // make sure that distMax is an upper bound
+ while(not oracle.isMatchLess(distMax)) {
+ // distMax is in fact a lower bound, so assign it to distMin
+ distMin = distMax;
+ // and increase distMax by 5 %
+ distMax = 1.05 * distMax;
+ }
+
+ // bounds are found, perform binary search
+ //std::cout << "Bounds found, distMin = " << distMin << ", distMax = " << distMax << ", ratio = " << ( distMax - distMin ) / distMin << std::endl ;
+ Real distProbe = ( distMin + distMax ) / 2.0;
+ while ( ( distMax - distMin ) / distMin >= epsilon ) {
+ if (oracle.isMatchLess(distProbe)) {
+ distMax = distProbe;
+ } else {
+ distMin = distProbe;
+ }
+ distProbe = ( distMin + distMax ) / 2.0;
+ }
+ return result;
+}
+
+// return the interval (distMin, distMax) such that:
+// a) actual bottleneck distance between A and B is contained in the interval
+// b) if the interval is not (0,0), then (distMax - distMin) / distMin < epsilon
+// use heuristic: initial estimate on sampled diagrams
+template<class Real>
+std::pair<Real, Real> bottleneckDistApproxIntervalHeur(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const Real epsilon)
+{
+ // empty diagrams are not considered as error
+ if (A.empty() and B.empty())
+ return std::make_pair(0.0, 0.0);
+
+ DiagramPointSet<Real> sampledA, sampledB;
+ sampleDiagramForHeur(A, sampledA);
+ sampleDiagramForHeur(B, sampledB);
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "A : " << A.size() << ", sampled: " << sampledA.size() << std::endl;
+ std::cout << "B : " << B.size() << ", sampled: " << sampledB.size() << std::endl;
+#endif
+ std::pair<Real, Real> initGuess = bottleneckDistApproxInterval(sampledA, sampledB, epsilon);
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "initial guess with sampling: " << initGuess.first << ", " << initGuess.second << std::endl;
+ std::cout << "running on the original diagrams" << std::endl;
+#endif
+ return bottleneckDistApproxIntervalWithInitial<Real>(A, B, epsilon, initGuess);
+}
+
+
+
+// get approximate distance,
+// see bottleneckDistApproxInterval
+template<class Real>
+Real bottleneckDistApprox(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const Real epsilon)
+{
+ auto interval = bottleneckDistApproxInterval<Real>(A, B, epsilon);
+ return interval.second;
+}
+
+
+template<class Real>
+Real bottleneckDistExactFromSortedPwDist(DiagramPointSet<Real>&A, DiagramPointSet<Real>& B, std::vector<Real>& pairwiseDist, const int decPrecision)
+{
+ //for(size_t k = 0; k < pairwiseDist.size(); ++k) {
+ //std::cout << "pairwiseDist[" << k << "] = " << std::setprecision(15) << pairwiseDist[k] << std::endl;
+ //}
+ // trivial case: we have only one candidate
+ if (pairwiseDist.size() == 1)
+ return pairwiseDist[0];
+
+ bool useRangeSearch = true;
+ Real distEpsilon = std::numeric_limits<Real>::max();
+ Real diffThreshold = 0.1;
+ for(int k = 0; k < decPrecision; ++k) {
+ diffThreshold /= 10.0;
+ }
+ for(size_t k = 0; k < pairwiseDist.size() - 2; ++k) {
+ auto diff = pairwiseDist[k+1]- pairwiseDist[k];
+ //std::cout << "diff = " << diff << ", pairwiseDist[k] = " << pairwiseDist[k] << std::endl;
+ if ( diff > diffThreshold and diff < distEpsilon ) {
+ distEpsilon = diff;
+ }
+ }
+ distEpsilon /= 3.0;
+ //std::cout << "decPrecision = " << decPrecision << ", distEpsilon = " << distEpsilon << std::endl;
+
+ BoundMatchOracle<Real> oracle(A, B, distEpsilon, useRangeSearch);
+ // binary search
+ size_t iterNum {0};
+ size_t idxMin {0}, idxMax {pairwiseDist.size() - 1};
+ size_t idxMid;
+ while(idxMax > idxMin) {
+ idxMid = static_cast<size_t>(floor(idxMin + idxMax) / 2.0);
+ //std::cout << "while begin: min = " << idxMin << ", idxMax = " << idxMax << ", idxMid = " << idxMid << ", testing d = " << std::setprecision(15) << pairwiseDist[idxMid] << std::endl;
+ iterNum++;
+ // not A[imid] < dist <=> A[imid] >= dist <=> A[imid[ >= dist + eps
+ if (oracle.isMatchLess(pairwiseDist[idxMid] + distEpsilon / 2.0)) {
+ //std::cout << "isMatchLess = true" << std::endl;
+ idxMax = idxMid;
+ } else {
+ //std::cout << "isMatchLess = false " << std::endl;
+ idxMin = idxMid + 1;
+ }
+ //std::cout << "while end: idxMin = " << idxMin << ", idxMax = " << idxMax << ", idxMid = " << idxMid << std::endl;
+ }
+ idxMid = static_cast<size_t>(floor(idxMin + idxMax) / 2.0);
+ return pairwiseDist[idxMid];
+}
+
+
+template<class Real>
+Real bottleneckDistExact(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B)
+{
+ return bottleneckDistExact(A, B, 14);
+}
+
+template<class Real>
+Real bottleneckDistExact(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B, const int decPrecision)
+{
+ using DgmPoint = DiagramPoint<Real>;
+
+ constexpr Real epsilon = 0.001;
+ auto interval = bottleneckDistApproxInterval(A, B, epsilon);
+ const Real delta = 0.50001 * (interval.second - interval.first);
+ const Real approxDist = 0.5 * ( interval.first + interval.second);
+ const Real minDist = interval.first;
+ const Real maxDist = interval.second;
+ //std::cout << std::setprecision(15) << "minDist = " << minDist << ", maxDist = " << maxDist << std::endl;
+ if ( delta == 0 ) {
+ return interval.first;
+ }
+ // copy points from A to a vector
+ // todo: get rid of this?
+ std::vector<DgmPoint> pointsA;
+ pointsA.reserve(A.size());
+ for(const auto& ptA : A) {
+ pointsA.push_back(ptA);
+ }
+
+ //std::vector<Real> killdist;
+ //for(auto pta : a) {
+ //for(auto ptb : b) {
+ //if ( distlinf(pta, ptb) > mindist and distlinf(pta, ptb) < maxdist) {
+ //killdist.push_back(distlinf(pta, ptb));
+ //std::cout << pta << ", " << ptb << std::endl;
+ //}
+ //}
+ //}
+ //std::sort(killdist.begin(), killdist.end());
+ //for(auto d : killdist) {
+ //std::cout << d << std::endl;
+ //}
+ //std::cout << "*************" << std::endl;
+
+ // in this vector we store the distances between the points
+ // that are candidates to realize
+ std::vector<Real> pairwiseDist;
+ {
+ // vector to store centers of vertical stripes
+ // two for each point in A and the id of the corresponding point
+ std::vector<std::pair<Real, DgmPoint>> xCentersVec;
+ xCentersVec.reserve(2 * pointsA.size());
+ for(auto ptA : pointsA) {
+ xCentersVec.push_back(std::make_pair(ptA.getRealX() - approxDist, ptA));
+ xCentersVec.push_back(std::make_pair(ptA.getRealX() + approxDist, ptA));
+ }
+ // lambda to compare pairs <coordinate, id> w.r.t coordinate
+ auto compLambda = [](std::pair<Real, DgmPoint> a, std::pair<Real, DgmPoint> b)
+ { return a.first < b.first; };
+
+ std::sort(xCentersVec.begin(), xCentersVec.end(), compLambda);
+ //std::cout << "xCentersVec.size = " << xCentersVec.size() << std::endl;
+ //for(auto p = xCentersVec.begin(); p!= xCentersVec.end(); ++p) {
+ //if (p->second.id == 200) {
+ //std::cout << "index of 200: " << p - xCentersVec.begin() << std::endl;
+ //}
+ //}
+ //std::vector<DgmPoint>
+ // todo: sort points in B, reduce search range in lower and upper bounds
+ for(auto ptB : B) {
+ // iterator to the first stripe such that ptB lies to the left
+ // from its right boundary (x_B <= x_j + \delta iff x_j >= x_B - \delta
+ auto itStart = std::lower_bound(xCentersVec.begin(),
+ xCentersVec.end(),
+ std::make_pair(ptB.getRealX() - delta, ptB),
+ compLambda);
+ //if (ptB.id == 236) {
+ //std::cout << itStart - xCentersVec.begin() << std::endl;
+ //}
+
+ for(auto iterA = itStart; iterA < xCentersVec.end(); ++iterA) {
+ //if (ptB.id == 236) {
+ //std::cout << "consider " << iterA->second << std::endl;
+ //}
+ if ( ptB.getRealX() < iterA->first - delta) {
+ // from that moment x_B >= x_j - delta
+ // is violated: x_B no longer lies to right from the left
+ // boundary of current stripe
+ //if (ptB.id == 236) {
+ //std::cout << "break" << std::endl;
+ //}
+ break;
+ }
+ // we're here => ptB lies in vertical stripe,
+ // check if distance fits into the interval we've found
+ Real pwDist = distLInf(iterA->second, ptB);
+ //if (ptB.id == 236) {
+ //std::cout << pwDist << std::endl;
+ //}
+ //std::cout << 1000*minDist << " <= " << 1000*pwDist << " <= " << 1000*maxDist << std::endl;
+ if (pwDist >= minDist and pwDist <= maxDist) {
+ pairwiseDist.push_back(pwDist);
+ }
+ }
+ }
+ }
+
+ {
+ // for y
+ // vector to store centers of vertical stripes
+ // two for each point in A and the id of the corresponding point
+ std::vector<std::pair<Real, DgmPoint>> yCentersVec;
+ yCentersVec.reserve(2 * pointsA.size());
+ for(auto ptA : pointsA) {
+ yCentersVec.push_back(std::make_pair(ptA.getRealY() - approxDist, ptA));
+ yCentersVec.push_back(std::make_pair(ptA.getRealY() + approxDist, ptA));
+ }
+ // lambda to compare pairs <coordinate, id> w.r.t coordinate
+ auto compLambda = [](std::pair<Real, DgmPoint> a, std::pair<Real, DgmPoint> b)
+ { return a.first < b.first; };
+
+ std::sort(yCentersVec.begin(), yCentersVec.end(), compLambda);
+
+ // std::cout << "Sorted vector of y-centers:" << std::endl;
+ //for(auto coordPtPair : yCentersVec) {
+ //std::cout << coordPtPair.first << ", id = " << coordPtPair.second.id << std::endl;
+ //}
+ /*std::cout << "End of sorted vector of y-centers:" << std::endl;*/
+
+ //std::vector<DgmPoint>
+ // todo: sort points in B, reduce search range in lower and upper bounds
+ for(auto ptB : B) {
+ auto itStart = std::lower_bound(yCentersVec.begin(),
+ yCentersVec.end(),
+ std::make_pair(ptB.getRealY() - delta, ptB),
+ compLambda);
+
+
+ for(auto iterA = itStart; iterA < yCentersVec.end(); ++iterA) {
+ if ( ptB.getRealY() < iterA->first - delta) {
+ break;
+ }
+ Real pwDist = distLInf(iterA->second, ptB);
+ //std::cout << 1000*minDist << " <= " << 1000*pwDist << " <= " << 1000*maxDist << std::endl;
+ if (pwDist >= minDist and pwDist <= maxDist) {
+ pairwiseDist.push_back(pwDist);
+ }
+ }
+ }
+ }
+
+ //std::cout << "pairwiseDist.size = " << pairwiseDist.size() << " out of " << A.size() * A.size() << std::endl;
+ std::sort(pairwiseDist.begin(), pairwiseDist.end());
+ //for(auto ddd : pairwiseDist) {
+ //std::cout << std::setprecision(15) << ddd << std::endl;
+ //}
+
+ return bottleneckDistExactFromSortedPwDist(A, B, pairwiseDist, decPrecision);
+}
+
+template<class Real>
+Real bottleneckDistSlow(DiagramPointSet<Real>& A, DiagramPointSet<Real>& B)
+{
+ using DistVerticesPair = std::pair<Real, std::pair<size_t, size_t>>;
+
+ // use range search when building the layer graph
+ bool useRangeSearch { true };
+ // find maximum of min. distances for each point,
+ // use this value as lower bound for bottleneck distance
+ bool useHeurMinIdx { true };
+
+ // find matching in a greedy manner to
+ // get an upper bound for a bottleneck distance
+ bool useHeurGreedyMatching { false };
+
+ // use successive multiplication of idxMin with 2 to get idxMax
+ bool goUpToFindIdxMax { false };
+ //
+ goUpToFindIdxMax = goUpToFindIdxMax and !useHeurGreedyMatching;
+
+ if (!useHeurGreedyMatching) {
+ long int N = 3 * (A.size() / 2 ) * (B.size() / 2);
+ std::vector<Real> pairwiseDist;
+ pairwiseDist.reserve(N);
+ Real maxMinDist {0.0};
+ for(auto& p_A : A) {
+ Real minDist { std::numeric_limits<Real>::max() };
+ for(auto& p_B : B) {
+ if (p_A.isNormal() or p_B.isNormal()) {
+ Real d = distLInf(p_A, p_B);
+ pairwiseDist.push_back(d);
+ if (useHeurMinIdx and p_A.isNormal()) {
+ if (d < minDist)
+ minDist = d;
+ }
+ }
+ }
+ if (useHeurMinIdx and p_A.isNormal() and minDist > maxMinDist) {
+ maxMinDist = minDist;
+ }
+ }
+
+ std::sort(pairwiseDist.begin(), pairwiseDist.end());
+
+ Real distEpsilon = std::numeric_limits<Real>::max();
+ for(size_t k = 0; k < pairwiseDist.size() - 2; ++k) {
+ auto diff = pairwiseDist[k+1]- pairwiseDist[k];
+ if ( diff > 1.0e-10 and diff < distEpsilon ) {
+ distEpsilon = diff;
+ }
+ }
+ distEpsilon /= 3.0;
+
+ BoundMatchOracle<Real> oracle(A, B, distEpsilon, useRangeSearch);
+ // binary search
+ size_t iterNum {0};
+ size_t idxMin {0}, idxMax {pairwiseDist.size() - 1};
+ if (useHeurMinIdx) {
+ auto maxMinIter = std::equal_range(pairwiseDist.begin(), pairwiseDist.end(), maxMinDist);
+ assert(maxMinIter.first != pairwiseDist.end());
+ idxMin = maxMinIter.first - pairwiseDist.begin();
+ //std::cout << "maxMinDist = " << maxMinDist << ", idxMin = " << idxMin << ", d = " << pairwiseDist[idxMin] << std::endl;
+ }
+
+ if (goUpToFindIdxMax) {
+ if ( pairwiseDist.size() == 1) {
+ return pairwiseDist[0];
+ }
+
+ idxMax = std::max<size_t>(idxMin, 1);
+ while (!oracle.isMatchLess(pairwiseDist[idxMax])) {
+ //std::cout << "entered while" << std::endl;
+ idxMin = idxMax;
+ if (2*idxMax > pairwiseDist.size() -1) {
+ idxMax = pairwiseDist.size() - 1;
+ break;
+ } else {
+ idxMax *= 2;
+ }
+ }
+ //std::cout << "size = " << pairwiseDist.size() << ", idxMax = " << idxMax << ", pw[max] = " << pairwiseDist[idxMax] << std::endl;
+ }
+
+ size_t idxMid { (idxMin + idxMax) / 2 };
+ while(idxMax > idxMin) {
+ iterNum++;
+ if (oracle.isMatchLess(pairwiseDist[idxMid])) {
+ idxMax = idxMid;
+ } else {
+ if (idxMax - idxMin == 1)
+ idxMin++;
+ else
+ idxMin = idxMid;
+ }
+ idxMid = (idxMin + idxMax) / 2;
+ }
+ return pairwiseDist[idxMid];
+ } else {
+ // with greeedy matching
+ long int N = A.size() * B.size();
+ std::vector<DistVerticesPair> pairwiseDist;
+ pairwiseDist.reserve(N);
+ Real maxMinDist {0.0};
+ size_t idxA{0}, idxB{0};
+ for(auto p_A : A) {
+ Real minDist { std::numeric_limits<Real>::max() };
+ idxB = 0;
+ for(auto p_B : B) {
+ Real d = distLInf(p_A, p_B);
+ pairwiseDist.push_back( std::make_pair(d, std::make_pair(idxA, idxB) ) );
+ if (useHeurMinIdx and p_A.isNormal()) {
+ if (d < minDist)
+ minDist = d;
+ }
+ idxB++;
+ }
+ if (useHeurMinIdx and p_A.isNormal() and minDist > maxMinDist) {
+ maxMinDist = minDist;
+ }
+ idxA++;
+ }
+
+ auto compLambda = [](DistVerticesPair a, DistVerticesPair b)
+ { return a.first < b.first;};
+
+ std::sort(pairwiseDist.begin(),
+ pairwiseDist.end(),
+ compLambda);
+
+ Real distEpsilon = std::numeric_limits<Real>::max();
+ for(size_t k = 0; k < pairwiseDist.size() - 2; ++k) {
+ auto diff = pairwiseDist[k+1].first - pairwiseDist[k].first;
+ if ( diff > 1.0e-10 and diff < distEpsilon ) {
+ distEpsilon = diff;
+ }
+ }
+ distEpsilon /= 3.0;
+
+ BoundMatchOracle<Real> oracle(A, B, distEpsilon, useRangeSearch);
+
+ // construct greedy matching
+ size_t numVert { A.size() };
+ size_t numMatched { 0 };
+ std::unordered_set<size_t> aTobMatched, bToaMatched;
+ aTobMatched.reserve(numVert);
+ bToaMatched.reserve(numVert);
+ size_t distVecIdx {0};
+ while( numMatched < numVert) {
+ auto vertPair = pairwiseDist[distVecIdx++].second;
+ //std::cout << "distVecIdx = " << distVecIdx << ", matched: " << numMatched << " out of " << numVert << std::endl;
+ //std::cout << "vertex A idx = " << vertPair.first << ", B idx: " << vertPair.second << " out of " << numVert << std::endl;
+ if ( aTobMatched.count(vertPair.first) == 0 and
+ bToaMatched.count(vertPair.second) == 0 ) {
+ aTobMatched.insert(vertPair.first);
+ bToaMatched.insert(vertPair.second);
+ numMatched++;
+ }
+ }
+ size_t idxMax = distVecIdx-1;
+ //std::cout << "idxMax = " << idxMax << ", size = " << pairwiseDist.size() << std::endl;
+ // binary search
+ size_t iterNum {0};
+ size_t idxMin {0};
+ if (useHeurMinIdx) {
+ auto maxMinIter = std::equal_range(pairwiseDist.begin(),
+ pairwiseDist.end(),
+ std::make_pair(maxMinDist, std::make_pair(0,0)),
+ compLambda);
+ assert(maxMinIter.first != pairwiseDist.end());
+ idxMin = maxMinIter.first - pairwiseDist.begin();
+ //std::cout << "maxMinDist = " << maxMinDist << ", idxMin = " << idxMin << ", d = " << pairwiseDist[idxMin].first << std::endl;
+ }
+ size_t idxMid { (idxMin + idxMax) / 2 };
+ while(idxMax > idxMin) {
+ iterNum++;
+ if (oracle.isMatchLess(pairwiseDist[idxMid].first)) {
+ idxMax = idxMid;
+ } else {
+ if (idxMax - idxMin == 1)
+ idxMin++;
+ else
+ idxMin = idxMid;
+ }
+ idxMid = (idxMin + idxMax) / 2;
+ }
+ return pairwiseDist[idxMid].first;
+ }
+ // stats
+ /*
+ // count number of edges
+ // pairwiseDist is sorted, add edges of the same length
+ int edgeNumber {idxMid};
+ while(pairwiseDist[edgeNumber + 1] == pairwiseDist[edgeNumber])
+ edgeNumber++;
+ // add edges between diagonal points
+ edgeNumber += N / 3;
+ // output stats
+ std::cout << idxMid << "\t" << N;
+ std::cout << "\t" << iterNum;
+ std::cout << "\t" << A.size() + B.size();
+ std::cout << "\t" << edgeNumber << "\t";
+ std::cout << (Real)(edgeNumber) / (Real)(A.size() + B.size()) << std::endl;
+ */
+}
+
+// wrappers
+template<class Real>
+bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<Real, Real>>& result)
+{
+ int decPrecision;
+ return readDiagramPointSet(fname.c_str(), result, decPrecision);
+}
+
+template<class Real>
+bool readDiagramPointSet(const char* fname, std::vector<std::pair<Real, Real>>& result)
+{
+ int decPrecision;
+ return readDiagramPointSet(fname, result, decPrecision);
+}
+
+template<class Real>
+bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<Real, Real>>& result, int& decPrecision)
+{
+ return readDiagramPointSet(fname.c_str(), result, decPrecision);
+}
+
+// reading function
+template<class Real>
+bool readDiagramPointSet(const char* fname, std::vector<std::pair<Real, Real>>& result, int& decPrecision)
+{
+ size_t lineNumber { 0 };
+ result.clear();
+ std::ifstream f(fname);
+ if (!f.good()) {
+#ifndef FOR_R_TDA
+ std::cerr << "Cannot open file " << fname << std::endl;
+#endif
+ return false;
+ }
+ std::string line;
+ while(std::getline(f, line)) {
+ lineNumber++;
+ // process comments: remove everything after hash
+ auto hashPos = line.find_first_of("#", 0);
+ if( std::string::npos != hashPos) {
+ line = std::string(line.begin(), line.begin() + hashPos);
+ }
+ if (line.empty()) {
+ continue;
+ }
+ // trim whitespaces
+ auto whiteSpaceFront = std::find_if_not(line.begin(),line.end(),isspace);
+ auto whiteSpaceBack = std::find_if_not(line.rbegin(),line.rend(),isspace).base();
+ if (whiteSpaceBack <= whiteSpaceFront) {
+ // line consists of spaces only - move to the next line
+ continue;
+ }
+ line = std::string(whiteSpaceFront,whiteSpaceBack);
+ bool fracPart = false;
+ int currDecPrecision = 0;
+ for(auto c : line) {
+ if (c == '.') {
+ fracPart = true;
+ } else if (fracPart) {
+ if (isdigit(c)) {
+ currDecPrecision++;
+ } else {
+ fracPart = false;
+ if (currDecPrecision > decPrecision)
+ decPrecision = currDecPrecision;
+ currDecPrecision = 0;
+ }
+ }
+ }
+ Real x, y;
+ std::istringstream iss(line);
+ if (not(iss >> x >> y)) {
+#ifndef FOR_R_TDA
+ std::cerr << "Error in file " << fname << ", line number " << lineNumber << ": cannot parse \"" << line << "\"" << std::endl;
+#endif
+ return false;
+ }
+ if ( x != y ) {
+ result.push_back(std::make_pair(x,y));
+ } else {
+#ifndef FOR_R_TDA
+#ifndef VERBOSE_BOTTLENECK
+ std::cerr << "Warning: in file " << fname << ", line number " << lineNumber << ", zero persistence point ignored: \"" << line << "\"" << std::endl;
+#endif
+#endif
+ }
+ }
+ f.close();
+ return true;
+}
+
+} // end namespace bt
+} // end namespace hera
+#endif // HERA_BOTTLENECK_HPP
diff --git a/geom_bottleneck/include/bound_match.h b/geom_bottleneck/include/bound_match.h
new file mode 100644
index 0000000..770c7df
--- /dev/null
+++ b/geom_bottleneck/include/bound_match.h
@@ -0,0 +1,107 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef HERA_BOUND_MATCH_H
+#define HERA_BOUND_MATCH_H
+
+#include <unordered_map>
+#include <memory>
+
+#include "basic_defs_bt.h"
+#include "neighb_oracle.h"
+
+
+namespace hera {
+namespace bt {
+
+template<class Real = double>
+class Matching {
+public:
+ using DgmPoint = DiagramPoint<Real>;
+ using DgmPointSet = DiagramPointSet<Real>;
+ using DgmPointHash = DiagramPointHash<Real>;
+ using Path = std::vector<DgmPoint>;
+
+ Matching(const DgmPointSet& AA, const DgmPointSet& BB) : A(AA), B(BB) {};
+ DgmPointSet getExposedVertices(bool forA = true) const ;
+ bool isExposed(const DgmPoint& p) const;
+ void getAllAdjacentVertices(const DgmPointSet& setIn, DgmPointSet& setOut, bool forA = true) const;
+ void increase(const Path& augmentingPath);
+ void checkAugPath(const Path& augPath) const;
+ bool getMatchedVertex(const DgmPoint& p, DgmPoint& result) const;
+ bool isPerfect() const;
+ void trimMatching(const Real newThreshold);
+#ifndef FOR_R_TDA
+ template<class R>
+ friend std::ostream& operator<<(std::ostream& output, const Matching<R>& m);
+#endif
+private:
+ DgmPointSet A;
+ DgmPointSet B;
+ std::unordered_map<DgmPoint, DgmPoint, DgmPointHash> AToB, BToA;
+ void matchVertices(const DgmPoint& pA, const DgmPoint& pB);
+ void sanityCheck() const;
+};
+
+
+
+template<class Real_ = double, class NeighbOracle_ = NeighbOracleDnn<Real_>>
+class BoundMatchOracle {
+public:
+ using Real = Real_;
+ using NeighbOracle = NeighbOracle_;
+ using DgmPoint = DiagramPoint<Real>;
+ using DgmPointSet = DiagramPointSet<Real>;
+ using Path = std::vector<DgmPoint>;
+
+ BoundMatchOracle(DgmPointSet psA, DgmPointSet psB, Real dEps, bool useRS = true);
+ bool isMatchLess(Real r);
+ bool buildMatchingForThreshold(const Real r);
+private:
+ DgmPointSet A, B;
+ Matching<Real> M;
+ void printLayerGraph();
+ void buildLayerGraph(Real r);
+ void buildLayerOracles(Real r);
+ bool buildAugmentingPath(const DgmPoint startVertex, Path& result);
+ void removeFromLayer(const DgmPoint& p, const int layerIdx);
+ std::unique_ptr<NeighbOracle> neighbOracle;
+ bool augPathExist;
+ std::vector<DgmPointSet> layerGraph;
+ std::vector<std::unique_ptr<NeighbOracle>> layerOracles;
+ Real distEpsilon;
+ bool useRangeSearch;
+ Real prevQueryValue;
+};
+
+} // end namespace bt
+} // end namespace hera
+
+#include "bound_match.hpp"
+
+#endif // HERA_BOUND_MATCH_H
diff --git a/geom_bottleneck/include/bound_match.hpp b/geom_bottleneck/include/bound_match.hpp
new file mode 100644
index 0000000..221bd0f
--- /dev/null
+++ b/geom_bottleneck/include/bound_match.hpp
@@ -0,0 +1,473 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef HERA_BOUND_MATCH_HPP
+#define HERA_BOUND_MATCH_HPP
+
+
+#ifdef FOR_R_TDA
+#undef DEBUG_BOUND_MATCH
+#undef DEBUG_MATCHING
+#undef VERBOSE_BOTTLENECK
+#endif
+
+#include <assert.h>
+#include "def_debug_bt.h"
+#include "bound_match.h"
+
+#ifdef VERBOSE_BOTTLENECK
+#include <chrono>
+#endif
+
+#ifndef FOR_R_TDA
+#include <iostream>
+#endif
+
+namespace hera{
+namespace bt {
+
+#ifndef FOR_R_TDA
+template<class Real>
+std::ostream& operator<<(std::ostream& output, const Matching<Real>& m)
+{
+ output << "Matching: " << m.AToB.size() << " pairs (";
+ if (!m.isPerfect()) {
+ output << "not";
+ }
+ output << " perfect)" << std::endl;
+ for(auto atob : m.AToB) {
+ output << atob.first << " <-> " << atob.second << " distance: " << distLInf(atob.first, atob.second) << std::endl;
+ }
+ return output;
+}
+#endif
+
+template<class R>
+void Matching<R>::sanityCheck() const
+{
+#ifdef DEBUG_MATCHING
+ assert( AToB.size() == BToA.size() );
+ for(auto aToBPair : AToB) {
+ auto bToAPair = BToA.find(aToBPair.second);
+ assert(bToAPair != BToA.end());
+ assert( aToBPair.first == bToAPair->second);
+ }
+#endif
+}
+
+template<class R>
+bool Matching<R>::isPerfect() const
+{
+ return AToB.size() == A.size();
+}
+
+template<class R>
+void Matching<R>::matchVertices(const DgmPoint& pA, const DgmPoint& pB)
+{
+ assert(A.hasElement(pA));
+ assert(B.hasElement(pB));
+ AToB.erase(pA);
+ AToB.insert( {{ pA, pB }} );
+ BToA.erase(pB);
+ BToA.insert( {{ pB, pA }} );
+}
+
+template<class R>
+bool Matching<R>::getMatchedVertex(const DgmPoint& p, DgmPoint& result) const
+{
+ sanityCheck();
+ auto inA = AToB.find(p);
+ if (inA != AToB.end()) {
+ result = inA->second;
+ return true;
+ } else {
+ auto inB = BToA.find(p);
+ if (inB != BToA.end()) {
+ result = inB->second;
+ return true;
+ }
+ }
+ return false;
+}
+
+
+template<class R>
+void Matching<R>::checkAugPath(const Path& augPath) const
+{
+ assert(augPath.size() % 2 == 0);
+ for(size_t idx = 0; idx < augPath.size(); ++idx) {
+ bool mustBeExposed { idx == 0 or idx == augPath.size() - 1 };
+ if (isExposed(augPath[idx]) != mustBeExposed) {
+#ifndef FOR_R_TDA
+ std::cerr << "mustBeExposed = " << mustBeExposed << ", idx = " << idx << ", point " << augPath[idx] << std::endl;
+#endif
+ }
+ assert( isExposed(augPath[idx]) == mustBeExposed );
+ DgmPoint matchedVertex {0.0, 0.0, DgmPoint::DIAG, 1};
+ if ( idx % 2 == 0 ) {
+ assert( A.hasElement(augPath[idx]));
+ if (not mustBeExposed) {
+ getMatchedVertex(augPath[idx], matchedVertex);
+ assert(matchedVertex == augPath[idx - 1]);
+ }
+ } else {
+ assert( B.hasElement(augPath[idx]));
+ if (not mustBeExposed) {
+ getMatchedVertex(augPath[idx], matchedVertex);
+ assert(matchedVertex == augPath[idx + 1]);
+ }
+ }
+ }
+}
+
+// use augmenting path to increase
+// the size of the matching
+template<class R>
+void Matching<R>::increase(const Path& augPath)
+{
+ sanityCheck();
+ // check that augPath is an augmenting path
+ checkAugPath(augPath);
+ for(size_t idx = 0; idx < augPath.size() - 1; idx += 2) {
+ matchVertices( augPath[idx], augPath[idx + 1]);
+ }
+ sanityCheck();
+}
+
+template<class R>
+DiagramPointSet<R> Matching<R>::getExposedVertices(bool forA) const
+{
+ sanityCheck();
+ DgmPointSet result;
+ const DgmPointSet* setToSearch { forA ? &A : &B };
+ const std::unordered_map<DgmPoint, DgmPoint, DgmPointHash>* mapToSearch { forA ? &AToB : &BToA };
+ for(auto it = setToSearch->cbegin(); it != setToSearch->cend(); ++it) {
+ if (mapToSearch->find((*it)) == mapToSearch->cend()) {
+ result.insert((*it));
+ }
+ }
+ return result;
+}
+
+template<class R>
+void Matching<R>::getAllAdjacentVertices(const DgmPointSet& setIn,
+ DgmPointSet& setOut,
+ bool forA) const
+{
+ sanityCheck();
+ //bool isDebug {false};
+ setOut.clear();
+ const std::unordered_map<DgmPoint, DgmPoint, DgmPointHash>* m;
+ m = ( forA ) ? &BToA : &AToB;
+ for(auto pit = setIn.cbegin(); pit != setIn.cend(); ++pit) {
+ auto findRes = m->find(*pit);
+ if (findRes != m->cend()) {
+ setOut.insert((*findRes).second);
+ }
+ }
+ sanityCheck();
+}
+
+template<class R>
+bool Matching<R>::isExposed(const DgmPoint& p) const
+{
+ return ( AToB.find(p) == AToB.end() ) && ( BToA.find(p) == BToA.end() );
+}
+
+// remove all edges whose length is > newThreshold
+template<class R>
+void Matching<R>::trimMatching(const R newThreshold)
+{
+ //bool isDebug { false };
+ sanityCheck();
+ for(auto aToBIter = AToB.begin(); aToBIter != AToB.end(); ) {
+ if ( distLInf(aToBIter->first, aToBIter->second) > newThreshold ) {
+ // remove edge from AToB and BToA
+ BToA.erase(aToBIter->second);
+ aToBIter = AToB.erase(aToBIter);
+ } else {
+ aToBIter++;
+ }
+ }
+ sanityCheck();
+}
+
+// ------- BoundMatchOracle --------------
+
+template<class R, class NO>
+BoundMatchOracle<R, NO>::BoundMatchOracle(DgmPointSet psA, DgmPointSet psB,
+ Real dEps, bool useRS) :
+ A(psA), B(psB), M(A, B), distEpsilon(dEps), useRangeSearch(useRS), prevQueryValue(0.0)
+{
+ neighbOracle = std::unique_ptr<NeighbOracle>(new NeighbOracle(psB, 0, distEpsilon));
+}
+
+template<class R, class NO>
+bool BoundMatchOracle<R, NO>::isMatchLess(Real r)
+{
+#ifdef VERBOSE_BOTTLENECK
+ std::chrono::high_resolution_clock hrClock;
+ std::chrono::time_point<std::chrono::high_resolution_clock> startMoment;
+ startMoment = hrClock.now();
+#endif
+ bool result = buildMatchingForThreshold(r);
+#ifdef VERBOSE_BOTTLENECK
+ auto endMoment = hrClock.now();
+ std::chrono::duration<double, std::milli> iterTime = endMoment - startMoment;
+ std::cout << "isMatchLess for r = " << r << " finished in " << std::chrono::duration<double, std::milli>(iterTime).count() << " ms." << std::endl;
+#endif
+ return result;
+
+}
+
+
+template<class R, class NO>
+void BoundMatchOracle<R, NO>::removeFromLayer(const DgmPoint& p, const int layerIdx) {
+ //bool isDebug {false};
+ layerGraph[layerIdx].erase(p);
+ if (layerOracles[layerIdx]) {
+ layerOracles[layerIdx]->deletePoint(p);
+ }
+}
+
+// return true, if there exists an augmenting path from startVertex
+// in this case the path is returned in result.
+// startVertex must be an exposed vertex from L_1 (layer[0])
+template<class R, class NO>
+bool BoundMatchOracle<R, NO>::buildAugmentingPath(const DgmPoint startVertex, Path& result)
+{
+ //bool isDebug {false};
+ DgmPoint prevVertexA = startVertex;
+ result.clear();
+ result.push_back(startVertex);
+ size_t evenLayerIdx {1};
+ while ( evenLayerIdx < layerGraph.size() ) {
+ //for(size_t evenLayerIdx = 1; evenLayerIdx < layerGraph.size(); evenLayerIdx += 2) {
+ DgmPoint nextVertexB{0.0, 0.0, DgmPoint::DIAG, 1}; // next vertex from even layer
+ bool neighbFound = layerOracles[evenLayerIdx]->getNeighbour(prevVertexA, nextVertexB);
+ if (neighbFound) {
+ result.push_back(nextVertexB);
+ if ( layerGraph.size() == evenLayerIdx + 1) {
+ break;
+ } else {
+ // nextVertexB must be matched with some vertex from the next odd
+ // layer
+ DgmPoint nextVertexA {0.0, 0.0, DgmPoint::DIAG, 1};
+ if (!M.getMatchedVertex(nextVertexB, nextVertexA)) {
+#ifndef FOR_R_TDA
+ std::cerr << "Vertices in even layers must be matched! Unmatched: ";
+ std::cerr << nextVertexB << std::endl;
+ std::cerr << evenLayerIdx << "; " << layerGraph.size() << std::endl;
+#endif
+ throw std::runtime_error("Unmatched vertex in even layer");
+ } else {
+ assert( ! (nextVertexA.getRealX() == 0 and nextVertexA.getRealY() == 0) );
+ result.push_back(nextVertexA);
+ prevVertexA = nextVertexA;
+ evenLayerIdx += 2;
+ continue;
+ }
+ }
+ } else {
+ // prevVertexA has no neighbours in the next layer,
+ // backtrack
+ if (evenLayerIdx == 1) {
+ // startVertex is not connected to any vertices
+ // in the next layer, augm. path doesn't exist
+ removeFromLayer(startVertex, 0);
+ return false;
+ } else {
+ assert(evenLayerIdx >= 3);
+ assert(result.size() % 2 == 1);
+ result.pop_back();
+ DgmPoint prevVertexB = result.back();
+ result.pop_back();
+ removeFromLayer(prevVertexA, evenLayerIdx-1);
+ removeFromLayer(prevVertexB, evenLayerIdx-2);
+ // we should proceed from the previous odd layer
+ assert(result.size() >= 1);
+ prevVertexA = result.back();
+ evenLayerIdx -= 2;
+ continue;
+ }
+ }
+ } // finished iterating over all layers
+ // remove all vertices in the augmenting paths
+ // the corresponding layers
+ for(size_t layerIdx = 0; layerIdx < result.size(); ++layerIdx) {
+ removeFromLayer(result[layerIdx], layerIdx);
+ }
+ return true;
+}
+
+
+
+
+template<class R, class NO>
+bool BoundMatchOracle<R, NO>::buildMatchingForThreshold(const Real r)
+{
+ //bool isDebug {false};
+ if (prevQueryValue > r) {
+ M.trimMatching(r);
+ }
+ prevQueryValue = r;
+ while(true) {
+ buildLayerGraph(r);
+ if (augPathExist) {
+ std::vector<Path> augmentingPaths;
+ DgmPointSet copyLG0;
+ for(DgmPoint p : layerGraph[0]) {
+ copyLG0.insert(p);
+ }
+ for(DgmPoint exposedVertex : copyLG0) {
+ Path augPath;
+ if (buildAugmentingPath(exposedVertex, augPath)) {
+ augmentingPaths.push_back(augPath);
+ }
+ }
+ if (augmentingPaths.empty()) {
+#ifndef FOR_R_TDA
+ std::cerr << "augmenting paths must exist, but were not found!" << std::endl;
+#endif
+ throw std::runtime_error("bad epsilon?");
+ }
+ // swap all augmenting paths with matching to increase it
+ for(auto& augPath : augmentingPaths ) {
+ M.increase(augPath);
+ }
+ } else {
+ return M.isPerfect();
+ }
+ }
+}
+
+template<class R, class NO>
+void BoundMatchOracle<R, NO>::printLayerGraph(void)
+{
+#ifdef DEBUG_BOUND_MATCH
+ for(auto& layer : layerGraph) {
+ std::cout << "{ ";
+ for(auto& p : layer) {
+ std::cout << p << "; ";
+ }
+ std::cout << "\b\b }" << std::endl;
+ }
+#endif
+}
+
+template<class R, class NO>
+void BoundMatchOracle<R, NO>::buildLayerGraph(Real r)
+{
+#ifdef VERBOSE_BOTTLENECK
+ std::cout << "Entered buildLayerGraph, r = " << r << std::endl;
+#endif
+ layerGraph.clear();
+ DgmPointSet L1 = M.getExposedVertices();
+ layerGraph.push_back(L1);
+ neighbOracle->rebuild(B, r);
+ size_t k = 0;
+ DgmPointSet layerNextEven;
+ DgmPointSet layerNextOdd;
+ bool exposedVerticesFound {false};
+ while(true) {
+ layerNextEven.clear();
+ for( auto p : layerGraph[k]) {
+ bool neighbFound;
+ DgmPoint neighbour {0.0, 0.0, DgmPoint::DIAG, 1};
+ if (useRangeSearch) {
+ std::vector<DgmPoint> neighbVec;
+ neighbOracle->getAllNeighbours(p, neighbVec);
+ neighbFound = !neighbVec.empty();
+ for(auto& neighbPt : neighbVec) {
+ layerNextEven.insert(neighbPt);
+ if (!exposedVerticesFound and M.isExposed(neighbPt))
+ exposedVerticesFound = true;
+ }
+ } else {
+ while(true) {
+ neighbFound = neighbOracle->getNeighbour(p, neighbour);
+ if (neighbFound) {
+ layerNextEven.insert(neighbour);
+ neighbOracle->deletePoint(neighbour);
+ if ((!exposedVerticesFound) && M.isExposed(neighbour)) {
+ exposedVerticesFound = true;
+ }
+ } else {
+ break;
+ }
+ }
+ } // without range search
+ } // all vertices from previous odd layer processed
+ if (layerNextEven.empty()) {
+ augPathExist = false;
+ break;
+ }
+ if (exposedVerticesFound) {
+ for(auto it = layerNextEven.cbegin(); it != layerNextEven.cend(); ) {
+ if ( ! M.isExposed(*it) ) {
+ layerNextEven.erase(it++);
+ } else {
+ ++it;
+ }
+
+ }
+ layerGraph.push_back(layerNextEven);
+ augPathExist = true;
+ break;
+ }
+ layerGraph.push_back(layerNextEven);
+ M.getAllAdjacentVertices(layerNextEven, layerNextOdd);
+ layerGraph.push_back(layerNextOdd);
+ k += 2;
+ }
+ buildLayerOracles(r);
+ printLayerGraph();
+ }
+
+// create geometric oracles for each even layer
+// odd layers have NULL in layerOracles
+template<class R, class NO>
+void BoundMatchOracle<R, NO>::buildLayerOracles(Real r)
+{
+ //bool isDebug {false};
+ // free previously constructed oracles
+ layerOracles.clear();
+ for(size_t layerIdx = 0; layerIdx < layerGraph.size(); ++layerIdx) {
+ if (layerIdx % 2 == 1) {
+ // even layer, build actual oracle
+ layerOracles.emplace_back(new NeighbOracle(layerGraph[layerIdx], r, distEpsilon));
+ } else {
+ // odd layer
+ layerOracles.emplace_back(nullptr);
+ }
+ }
+}
+
+} // end namespace bt
+} // end namespace hera
+#endif // HERA_BOUND_MATCH_HPP
diff --git a/geom_bottleneck/include/def_debug_bt.h b/geom_bottleneck/include/def_debug_bt.h
new file mode 100644
index 0000000..21557e7
--- /dev/null
+++ b/geom_bottleneck/include/def_debug_bt.h
@@ -0,0 +1,42 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef DEF_DEBUG_BT_H
+#define DEF_DEBUG_BT_H
+
+//#define DEBUG_BOUND_MATCH
+//#define DEBUG_NEIGHBOUR_ORACLE
+//#define DEBUG_MATCHING
+//#define DEBUG_AUCTION
+// This symbol should be defined only in the version
+// for R package TDA, to comply with some CRAN rules
+// like no usage of cout, cerr, cin, exit, etc.
+//#define FOR_R_TDA
+//#define VERBOSE_BOTTLENECK
+
+#endif
diff --git a/geom_bottleneck/include/diagram_reader.h b/geom_bottleneck/include/diagram_reader.h
new file mode 100644
index 0000000..5bf106d
--- /dev/null
+++ b/geom_bottleneck/include/diagram_reader.h
@@ -0,0 +1,139 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef HERA_DIAGRAM_READER_H
+#define HERA_DIAGRAM_READER_H
+
+#ifndef FOR_R_TDA
+#include <iostream>
+#endif
+
+#include <iomanip>
+#include <sstream>
+#include <string>
+#include <cctype>
+
+namespace hera {
+// fill in result with points from file fname
+// return false if file can't be opened
+// or error occurred while reading
+// decPrecision is the maximal decimal precision in the input,
+// it is zero if all coordinates in the input are integers
+template<class RealType_ = double, class ContType_ = std::vector<std::pair<RealType_, RealType_>>>
+bool readDiagramPointSet(const char* fname, ContType_& result, int& decPrecision)
+{
+ size_t lineNumber { 0 };
+ result.clear();
+ std::ifstream f(fname);
+ if (!f.good()) {
+#ifndef FOR_R_TDA
+ std::cerr << "Cannot open file " << fname << std::endl;
+#endif
+ return false;
+ }
+ std::string line;
+ while(std::getline(f, line)) {
+ lineNumber++;
+ // process comments: remove everything after hash
+ auto hashPos = line.find_first_of("#", 0);
+ if( std::string::npos != hashPos) {
+ line = std::string(line.begin(), line.begin() + hashPos);
+ }
+ if (line.empty()) {
+ continue;
+ }
+ // trim whitespaces
+ auto whiteSpaceFront = std::find_if_not(line.begin(),line.end(),isspace);
+ auto whiteSpaceBack = std::find_if_not(line.rbegin(),line.rend(),isspace).base();
+ if (whiteSpaceBack <= whiteSpaceFront) {
+ // line consists of spaces only - move to the next line
+ continue;
+ }
+ line = std::string(whiteSpaceFront,whiteSpaceBack);
+ bool fracPart = false;
+ int currDecPrecision = 0;
+ for(auto c : line) {
+ if (c == '.') {
+ fracPart = true;
+ } else if (fracPart) {
+ if (isdigit(c)) {
+ currDecPrecision++;
+ } else {
+ fracPart = false;
+ if (currDecPrecision > decPrecision)
+ decPrecision = currDecPrecision;
+ currDecPrecision = 0;
+ }
+ }
+ }
+ RealType_ x, y;
+ std::istringstream iss(line);
+ if (not(iss >> x >> y)) {
+#ifndef FOR_R_TDA
+ std::cerr << "Error in file " << fname << ", line number " << lineNumber << ": cannot parse \"" << line << "\"" << std::endl;
+#endif
+ return false;
+ }
+ if ( x != y ) {
+ result.push_back(std::make_pair(x,y));
+ } else {
+#ifdef VERBOSE_BOTTLENECK
+ std::cerr << "Warning: in file " << fname << ", line number " << lineNumber << ", zero persistence point ignored: \"" << line << "\"" << std::endl;
+#endif
+ }
+ }
+ f.close();
+ return true;
+}
+
+
+// wrappers
+template<class RealType_ = double, class ContType_ = std::vector<std::pair<RealType_, RealType_>>>
+bool readDiagramPointSet(const std::string& fname, ContType_& result, int& decPrecision)
+{
+ return readDiagramPointSet<RealType_, ContType_>(fname.c_str(), result, decPrecision);
+}
+
+// these two functions are now just wrappers for the previous ones,
+// in case someone needs them; decPrecision is ignored
+template<class RealType_ = double, class ContType_ = std::vector<std::pair<RealType_, RealType_>>>
+bool readDiagramPointSet(const char* fname, ContType_& result)
+{
+ int decPrecision;
+ return readDiagramPointSet<RealType_, ContType_>(fname, result, decPrecision);
+}
+
+template<class RealType_ = double, class ContType_ = std::vector<std::pair<RealType_, RealType_>>>
+bool readDiagramPointSet(const std::string& fname, ContType_& result)
+{
+ int decPrecision;
+ return readDiagramPointSet<RealType_, ContType_>(fname.c_str(), result, decPrecision);
+}
+
+} // end namespace hera
+#endif // HERA_DIAGRAM_READER_H
diff --git a/geom_bottleneck/include/diagram_traits.h b/geom_bottleneck/include/diagram_traits.h
new file mode 100644
index 0000000..c8d4862
--- /dev/null
+++ b/geom_bottleneck/include/diagram_traits.h
@@ -0,0 +1,45 @@
+#ifndef HERA_DIAGRAM_TRAITS_H
+#define HERA_DIAGRAM_TRAITS_H
+
+namespace hera {
+
+template<class PairContainer_, class PointType_ = typename std::remove_reference< decltype(*std::declval<PairContainer_>().begin())>::type >
+struct DiagramTraits
+{
+ using Container = PairContainer_;
+ using PointType = PointType_;
+ using RealType = typename std::remove_reference< decltype(std::declval<PointType>()[0]) >::type;
+
+ static RealType get_x(const PointType& p) { return p[0]; }
+ static RealType get_y(const PointType& p) { return p[1]; }
+};
+
+
+template<class PairContainer_>
+struct DiagramTraits<PairContainer_, std::pair<double, double>>
+{
+ using PointType = std::pair<double, double>;
+ using RealType = double;
+ using Container = std::vector<PointType>;
+
+ static RealType get_x(const PointType& p) { return p.first; }
+ static RealType get_y(const PointType& p) { return p.second; }
+};
+
+
+template<class PairContainer_>
+struct DiagramTraits<PairContainer_, std::pair<float, float>>
+{
+ using PointType = std::pair<float, float>;
+ using RealType = float;
+ using Container = std::vector<PointType>;
+
+ static RealType get_x(const PointType& p) { return p.first; }
+ static RealType get_y(const PointType& p) { return p.second; }
+};
+
+
+} // end namespace hera
+
+
+#endif // HERA_DIAGRAM_TRAITS_H
diff --git a/geom_bottleneck/include/dnn/geometry/euclidean-fixed.h b/geom_bottleneck/include/dnn/geometry/euclidean-fixed.h
new file mode 100644
index 0000000..f45b980
--- /dev/null
+++ b/geom_bottleneck/include/dnn/geometry/euclidean-fixed.h
@@ -0,0 +1,162 @@
+#ifndef HERA_BT_DNN_GEOMETRY_EUCLIDEAN_FIXED_H
+#define HERA_BT_DNN_GEOMETRY_EUCLIDEAN_FIXED_H
+
+#include <boost/operators.hpp>
+#include <boost/array.hpp>
+#include <boost/range/value_type.hpp>
+#include <boost/serialization/access.hpp>
+#include <boost/serialization/base_object.hpp>
+
+#include <fstream>
+#include <string>
+#include <sstream>
+#include <cmath>
+
+#include "../parallel/tbb.h" // for dnn::vector<...>
+
+namespace hera {
+namespace bt {
+namespace dnn
+{
+ // TODO: wrap in another namespace (e.g., euclidean)
+
+ template<size_t D, typename Real = double>
+ struct Point:
+ boost::addable< Point<D,Real>,
+ boost::subtractable< Point<D,Real>,
+ boost::dividable2< Point<D, Real>, Real,
+ boost::multipliable2< Point<D, Real>, Real > > > >,
+ public boost::array<Real, D>
+ {
+ public:
+ typedef Real Coordinate;
+ typedef Real DistanceType;
+
+
+ public:
+ Point(size_t id = 0): id_(id) {}
+ template<size_t DD>
+ Point(const Point<DD,Real>& p, size_t id = 0):
+ id_(id) { *this = p; }
+
+ static size_t dimension() { return D; }
+
+ // Assign a point of different dimension
+ template<size_t DD>
+ Point& operator=(const Point<DD,Real>& p) { for (size_t i = 0; i < (D < DD ? D : DD); ++i) (*this)[i] = p[i]; if (DD < D) for (size_t i = DD; i < D; ++i) (*this)[i] = 0; return *this; }
+
+ Point& operator+=(const Point& p) { for (size_t i = 0; i < D; ++i) (*this)[i] += p[i]; return *this; }
+ Point& operator-=(const Point& p) { for (size_t i = 0; i < D; ++i) (*this)[i] -= p[i]; return *this; }
+ Point& operator/=(Real r) { for (size_t i = 0; i < D; ++i) (*this)[i] /= r; return *this; }
+ Point& operator*=(Real r) { for (size_t i = 0; i < D; ++i) (*this)[i] *= r; return *this; }
+
+ Real norm2() const { Real n = 0; for (size_t i = 0; i < D; ++i) n += (*this)[i] * (*this)[i]; return n; }
+ Real max_norm() const
+ {
+ Real m = std::fabs((*this)[0]);
+ for (size_t i = 1; i < D; ++i)
+ if (std::fabs((*this)[i]) > m)
+ m = std::fabs((*this)[i]);
+ return m; }
+
+ // quick and dirty for now; make generic later
+ //DistanceType distance(const Point& other) const { return sqrt(sq_distance(other)); }
+ //DistanceType sq_distance(const Point& other) const { return (other - *this).norm2(); }
+
+ DistanceType distance(const Point& other) const { return (other - *this).max_norm(); }
+ DistanceType sq_distance(const Point& other) const { DistanceType d = distance(other); return d*d; }
+
+ size_t id() const { return id_; }
+ size_t& id() { return id_; }
+
+ private:
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int version) { ar & boost::serialization::base_object< boost::array<Real,D> >(*this) & id_; }
+
+ private:
+ size_t id_;
+ };
+
+ template<size_t D, typename Real>
+ std::ostream&
+ operator<<(std::ostream& out, const Point<D,Real>& p)
+ { out << p[0]; for (size_t i = 1; i < D; ++i) out << " " << p[i]; return out; }
+
+
+ template<class Point>
+ struct PointTraits; // intentionally undefined; should be specialized for each type
+
+ template<size_t D, typename Real>
+ struct PointTraits< Point<D, Real> > // specialization for dnn::Point
+ {
+ typedef Point<D,Real> PointType;
+ typedef const PointType* PointHandle;
+ typedef std::vector<PointType> PointContainer;
+
+ typedef typename PointType::Coordinate Coordinate;
+ typedef typename PointType::DistanceType DistanceType;
+
+ static DistanceType
+ distance(const PointType& p1, const PointType& p2) { return p1.distance(p2); }
+ static DistanceType
+ distance(PointHandle p1, PointHandle p2) { return distance(*p1,*p2); }
+ static DistanceType
+ sq_distance(const PointType& p1,
+ const PointType& p2) { return p1.sq_distance(p2); }
+ static DistanceType
+ sq_distance(PointHandle p1, PointHandle p2) { return sq_distance(*p1,*p2); }
+ static size_t dimension() { return D; }
+ static Real coordinate(const PointType& p, size_t i) { return p[i]; }
+ static Real& coordinate(PointType& p, size_t i) { return p[i]; }
+ static Real coordinate(PointHandle p, size_t i) { return coordinate(*p,i); }
+
+ static size_t id(const PointType& p) { return p.id(); }
+ static size_t& id(PointType& p) { return p.id(); }
+ static size_t id(PointHandle p) { return id(*p); }
+
+ static PointHandle
+ handle(const PointType& p) { return &p; }
+ static const PointType&
+ point(PointHandle ph) { return *ph; }
+
+ void swap(PointType& p1, PointType& p2) const { return std::swap(p1, p2); }
+
+ static PointContainer
+ container(size_t n = 0, const PointType& p = PointType()) { return PointContainer(n, p); }
+ static typename PointContainer::iterator
+ iterator(PointContainer& c, PointHandle ph) { return c.begin() + (ph - &c[0]); }
+ static typename PointContainer::const_iterator
+ iterator(const PointContainer& c, PointHandle ph) { return c.begin() + (ph - &c[0]); }
+
+ private:
+ friend class boost::serialization::access;
+
+ template<class Archive>
+ void serialize(Archive& ar, const unsigned int version) {}
+ };
+
+ template<class PointContainer>
+ void read_points(const std::string& filename, PointContainer& points)
+ {
+ typedef typename boost::range_value<PointContainer>::type Point;
+ typedef typename PointTraits<Point>::Coordinate Coordinate;
+
+ std::ifstream in(filename.c_str());
+ std::string line;
+ while(std::getline(in, line))
+ {
+ if (line[0] == '#') continue; // comment line in the file
+ std::stringstream linestream(line);
+ Coordinate x;
+ points.push_back(Point());
+ size_t i = 0;
+ while (linestream >> x)
+ points.back()[i++] = x;
+ }
+ }
+} // dnn
+} // bt
+} // hera
+#endif
diff --git a/geom_bottleneck/include/dnn/local/kd-tree.h b/geom_bottleneck/include/dnn/local/kd-tree.h
new file mode 100644
index 0000000..c1aed2b
--- /dev/null
+++ b/geom_bottleneck/include/dnn/local/kd-tree.h
@@ -0,0 +1,106 @@
+#ifndef HERA_BT_DNN_LOCAL_KD_TREE_H
+#define HERA_BT_DNN_LOCAL_KD_TREE_H
+
+#include "../utils.h"
+#include "search-functors.h"
+
+#include <unordered_map>
+#include <stack>
+
+#include <boost/tuple/tuple.hpp>
+#include <boost/shared_ptr.hpp>
+#include <boost/range/value_type.hpp>
+
+#include <boost/static_assert.hpp>
+#include <boost/type_traits.hpp>
+
+namespace hera {
+namespace bt {
+namespace dnn
+{
+ // Weighted KDTree
+ // Traits_ provides Coordinate, DistanceType, PointType, dimension(), distance(p1,p2), coordinate(p,i)
+ template< class Traits_ >
+ class KDTree
+ {
+ public:
+ typedef Traits_ Traits;
+ typedef hera::bt::dnn::HandleDistance<KDTree> HandleDistance;
+
+ typedef typename Traits::PointType Point;
+ typedef typename Traits::PointHandle PointHandle;
+ typedef typename Traits::Coordinate Coordinate;
+ typedef typename Traits::DistanceType DistanceType;
+ typedef std::vector<PointHandle> HandleContainer;
+ typedef std::vector<HandleDistance> HDContainer; // TODO: use tbb::scalable_allocator
+ typedef HDContainer Result;
+ typedef std::vector<DistanceType> DistanceContainer;
+ typedef std::unordered_map<PointHandle, size_t> HandleMap;
+ //private:
+ typedef typename HandleContainer::iterator HCIterator;
+ typedef std::tuple<HCIterator, HCIterator, size_t, ssize_t> KDTreeNode;
+ typedef std::tuple<HCIterator, HCIterator> KDTreeNodeNoCut;
+
+ //BOOST_STATIC_ASSERT_MSG(has_coordinates<Traits, PointHandle, int>::value, "KDTree requires coordinates");
+
+ public:
+ KDTree(const Traits& traits):
+ traits_(traits) {}
+
+ KDTree(const Traits& traits, HandleContainer&& handles);
+
+ template<class Range>
+ KDTree(const Traits& traits, const Range& range);
+
+ template<class Range>
+ void init(const Range& range);
+
+ HandleDistance find(PointHandle q) const;
+ Result findR(PointHandle q, DistanceType r) const; // all neighbors within r
+ Result findFirstR(PointHandle q, DistanceType r) const; // first neighbor within r
+ Result findK(PointHandle q, size_t k) const; // k nearest neighbors
+
+ HandleDistance find(const Point& q) const { return find(traits().handle(q)); }
+ Result findR(const Point& q, DistanceType r) const { return findR(traits().handle(q), r); }
+ Result findFirstR(const Point& q, DistanceType r) const { return findFirstR(traits().handle(q), r); }
+ Result findK(const Point& q, size_t k) const { return findK(traits().handle(q), k); }
+
+
+
+ template<class ResultsFunctor>
+ void search(PointHandle q, ResultsFunctor& rf) const;
+
+ const Traits& traits() const { return traits_; }
+
+ void get_path_to_root(const size_t idx, std::stack<KDTreeNodeNoCut>& s);
+ // to support deletion
+ void init_n_elems();
+ void delete_point(const size_t idx);
+ void delete_point(PointHandle p);
+ void update_n_elems(const ssize_t idx, const int delta);
+ void increase_n_elems(const ssize_t idx);
+ void decrease_n_elems(const ssize_t idx);
+ size_t get_num_points() const { return num_points_; }
+ //private:
+ void init();
+
+
+ struct CoordinateComparison;
+ struct OrderTree;
+
+ //private:
+ Traits traits_;
+ HandleContainer tree_;
+ std::vector<char> delete_flags_;
+ std::vector<int> subtree_n_elems;
+ HandleMap indices_;
+ std::vector<ssize_t> parents_;
+
+ size_t num_points_;
+ };
+} // dnn
+} // bt
+} // hera
+#include "kd-tree.hpp"
+
+#endif
diff --git a/geom_bottleneck/include/dnn/local/kd-tree.hpp b/geom_bottleneck/include/dnn/local/kd-tree.hpp
new file mode 100644
index 0000000..249fa55
--- /dev/null
+++ b/geom_bottleneck/include/dnn/local/kd-tree.hpp
@@ -0,0 +1,296 @@
+#include <boost/range/counting_range.hpp>
+#include <boost/range/algorithm_ext/push_back.hpp>
+#include <boost/range.hpp>
+
+#include <queue>
+
+#include "../parallel/tbb.h"
+
+template<class T>
+hera::bt::dnn::KDTree<T>::KDTree(const Traits& traits, HandleContainer&& handles):
+ traits_(traits),
+ tree_(std::move(handles)),
+ delete_flags_(handles.size(), static_cast<char>(0) ),
+ subtree_n_elems(handles.size(), static_cast<size_t>(0)),
+ num_points_(handles.size())
+{
+ init();
+}
+
+template<class T>
+template<class Range>
+hera::bt::dnn::KDTree<T>::KDTree(const Traits& traits, const Range& range):
+ traits_(traits)
+{
+ init(range);
+}
+
+template<class T>
+template<class Range>
+void hera::bt::dnn::KDTree<T>::init(const Range& range)
+{
+ size_t sz = std::distance(std::begin(range), std::end(range));
+ subtree_n_elems = std::vector<int>(sz, 0);
+ delete_flags_ = std::vector<char>(sz, 0);
+ num_points_ = sz;
+ tree_.reserve(sz);
+ for (PointHandle h : range)
+ tree_.push_back(h);
+ parents_.resize(sz, -1);
+ init();
+}
+
+template<class T>
+void hera::bt::dnn::KDTree<T>::init()
+{
+ if (tree_.empty())
+ return;
+
+#if defined(TBB)
+ task_group g;
+ g.run(OrderTree(this, tree_.begin(), tree_.end(), -1, 0, traits()));
+ g.wait();
+#else
+ OrderTree(this, tree_.begin(), tree_.end(), -1, 0, traits()).serial();
+#endif
+
+ for (size_t i = 0; i < tree_.size(); ++i)
+ indices_[tree_[i]] = i;
+ init_n_elems();
+}
+
+template<class T>
+struct
+hera::bt::dnn::KDTree<T>::OrderTree
+{
+ OrderTree(KDTree* tree_, HCIterator b_, HCIterator e_, ssize_t p_, size_t i_, const Traits& traits_):
+ tree(tree_), b(b_), e(e_), p(p_), i(i_), traits(traits_) {}
+
+ void operator()() const
+ {
+ if (e - b < 1000)
+ {
+ serial();
+ return;
+ }
+
+ HCIterator m = b + (e - b)/2;
+ ssize_t im = m - tree->tree_.begin();
+ tree->parents_[im] = p;
+
+ CoordinateComparison cmp(i, traits);
+ std::nth_element(b,m,e, cmp);
+ size_t next_i = (i + 1) % traits.dimension();
+
+ task_group g;
+ if (b < m - 1) g.run(OrderTree(tree, b, m, im, next_i, traits));
+ if (e > m + 2) g.run(OrderTree(tree, m+1, e, im, next_i, traits));
+ g.wait();
+ }
+
+ void serial() const
+ {
+ std::queue<KDTreeNode> q;
+ q.push(KDTreeNode(b,e,p,i));
+ while (!q.empty())
+ {
+ HCIterator b, e; ssize_t p; size_t i;
+ std::tie(b,e,p,i) = q.front();
+ q.pop();
+ HCIterator m = b + (e - b)/2;
+ ssize_t im = m - tree->tree_.begin();
+ tree->parents_[im] = p;
+
+ CoordinateComparison cmp(i, traits);
+ std::nth_element(b,m,e, cmp);
+ size_t next_i = (i + 1) % traits.dimension();
+
+ // Replace with a size condition instead?
+ if (b < m - 1)
+ q.push(KDTreeNode(b, m, im, next_i));
+ else if (b < m)
+ tree->parents_[im - 1] = im;
+ if (e > m + 2)
+ q.push(KDTreeNode(m+1, e, im, next_i));
+ else if (e > m + 1)
+ tree->parents_[im + 1] = im;
+ }
+ }
+
+ KDTree* tree;
+ HCIterator b, e;
+ ssize_t p;
+ size_t i;
+ const Traits& traits;
+};
+
+template<class T>
+void hera::bt::dnn::KDTree<T>::update_n_elems(ssize_t idx, const int delta)
+// add delta to the number of points in node idx and update subtree_n_elems
+// for all parents of the node idx
+{
+ //std::cout << "subtree_n_elems.size = " << subtree_n_elems.size() << std::endl;
+ // update the node itself
+ while (idx != -1)
+ {
+ //std::cout << idx << std::endl;
+ subtree_n_elems[idx] += delta;
+ idx = parents_[idx];
+ }
+}
+
+template<class T>
+void hera::bt::dnn::KDTree<T>::increase_n_elems(const ssize_t idx)
+{
+ update_n_elems(idx, static_cast<ssize_t>(1));
+}
+
+template<class T>
+void hera::bt::dnn::KDTree<T>::decrease_n_elems(const ssize_t idx)
+{
+ update_n_elems(idx, static_cast<ssize_t>(-1));
+}
+
+template<class T>
+void hera::bt::dnn::KDTree<T>::init_n_elems()
+{
+ for(size_t idx = 0; idx < tree_.size(); ++idx) {
+ increase_n_elems(idx);
+ }
+}
+
+
+template<class T>
+template<class ResultsFunctor>
+void hera::bt::dnn::KDTree<T>::search(PointHandle q, ResultsFunctor& rf) const
+{
+ typedef typename HandleContainer::const_iterator HCIterator;
+ typedef std::tuple<HCIterator, HCIterator, size_t> KDTreeNode;
+
+ if (tree_.empty())
+ return;
+
+ DistanceType D = std::numeric_limits<DistanceType>::infinity();
+
+ // TODO: use tbb::scalable_allocator for the queue
+ std::queue<KDTreeNode> nodes;
+
+ nodes.push(KDTreeNode(tree_.begin(), tree_.end(), 0));
+
+ //std::cout << "started kdtree::search" << std::endl;
+
+ while (!nodes.empty())
+ {
+ HCIterator b, e; size_t i;
+ std::tie(b,e,i) = nodes.front();
+ nodes.pop();
+
+ CoordinateComparison cmp(i, traits());
+ i = (i + 1) % traits().dimension();
+
+ HCIterator m = b + (e - b)/2;
+ size_t m_idx = m - tree_.begin();
+ // ignore deleted points
+ if ( delete_flags_[m_idx] == 0 ) {
+ DistanceType dist = traits().distance(q, *m);
+ // + weights_[m - tree_.begin()];
+ //std::cout << "Supplied to functor: m : ";
+ //std::cout << "(" << (*(*m))[0] << ", " << (*(*m))[1] << ")";
+ //std::cout << " and q : ";
+ //std::cout << "(" << (*q)[0] << ", " << (*q)[1] << ")" << std::endl;
+ //std::cout << "dist^q + weight = " << dist << std::endl;
+ //std::cout << "weight = " << weights_[m - tree_.begin()] << std::endl;
+ //std::cout << "dist = " << traits().distance(q, *m) << std::endl;
+ //std::cout << "dist^q = " << pow(traits().distance(q, *m), wassersteinPower) << std::endl;
+
+ D = rf(*m, dist);
+ }
+ // we are really searching w.r.t L_\infty ball; could prune better with an L_2 ball
+ Coordinate diff = cmp.diff(q, *m); // diff returns signed distance
+ DistanceType diffToWasserPower = (diff > 0 ? 1.0 : -1.0) * fabs(diff);
+
+ size_t lm = m + 1 + (e - (m+1))/2 - tree_.begin();
+ if ( subtree_n_elems[lm] > 0 ) {
+ if (e > m + 1 && diffToWasserPower >= -D) {
+ nodes.push(KDTreeNode(m+1, e, i));
+ }
+ }
+
+ size_t rm = b + (m - b) / 2 - tree_.begin();
+ if ( subtree_n_elems[rm] > 0 ) {
+ if (b < m && diffToWasserPower <= D) {
+ nodes.push(KDTreeNode(b, m, i));
+ }
+ }
+ }
+ //std::cout << "exited kdtree::search" << std::endl;
+}
+
+template<class T>
+typename hera::bt::dnn::KDTree<T>::HandleDistance hera::bt::dnn::KDTree<T>::find(PointHandle q) const
+{
+ hera::bt::dnn::NNRecord<HandleDistance> nn;
+ search(q, nn);
+ return nn.result;
+}
+
+template<class T>
+typename hera::bt::dnn::KDTree<T>::Result hera::bt::dnn::KDTree<T>::findR(PointHandle q, DistanceType r) const
+{
+ hera::bt::dnn::rNNRecord<HandleDistance> rnn(r);
+ search(q, rnn);
+ //std::sort(rnn.result.begin(), rnn.result.end());
+ return rnn.result;
+}
+
+template<class T>
+typename hera::bt::dnn::KDTree<T>::Result hera::bt::dnn::KDTree<T>::findFirstR(PointHandle q, DistanceType r) const
+{
+ hera::bt::dnn::firstrNNRecord<HandleDistance> rnn(r);
+ search(q, rnn);
+ return rnn.result;
+}
+
+template<class T>
+typename hera::bt::dnn::KDTree<T>::Result hera::bt::dnn::KDTree<T>::findK(PointHandle q, size_t k) const
+{
+ hera::bt::dnn::kNNRecord<HandleDistance> knn(k);
+ search(q, knn);
+ // do we need this???
+ std::sort(knn.result.begin(), knn.result.end());
+ return knn.result;
+}
+
+template<class T>
+struct hera::bt::dnn::KDTree<T>::CoordinateComparison
+{
+ CoordinateComparison(size_t i, const Traits& traits):
+ i_(i), traits_(traits) {}
+
+ bool operator()(PointHandle p1, PointHandle p2) const { return coordinate(p1) < coordinate(p2); }
+ Coordinate diff(PointHandle p1, PointHandle p2) const { return coordinate(p1) - coordinate(p2); }
+
+ Coordinate coordinate(PointHandle p) const { return traits_.coordinate(p, i_); }
+ size_t axis() const { return i_; }
+
+ private:
+ size_t i_;
+ const Traits& traits_;
+};
+
+template<class T>
+void hera::bt::dnn::KDTree<T>::delete_point(const size_t idx)
+{
+ // prevent double deletion
+ assert(delete_flags_[idx] == 0);
+ delete_flags_[idx] = 1;
+ decrease_n_elems(idx);
+ --num_points_;
+}
+
+template<class T>
+void hera::bt::dnn::KDTree<T>::delete_point(PointHandle p)
+{
+ delete_point(indices_[p]);
+}
+
diff --git a/geom_bottleneck/include/dnn/local/search-functors.h b/geom_bottleneck/include/dnn/local/search-functors.h
new file mode 100644
index 0000000..63ad11d
--- /dev/null
+++ b/geom_bottleneck/include/dnn/local/search-functors.h
@@ -0,0 +1,119 @@
+#ifndef HERA_BT_DNN_LOCAL_SEARCH_FUNCTORS_H
+#define HERA_BT_DNN_LOCAL_SEARCH_FUNCTORS_H
+
+#include <boost/range/algorithm/heap_algorithm.hpp>
+
+namespace hera
+{
+namespace bt
+{
+namespace dnn
+{
+
+template<class NN>
+struct HandleDistance
+{
+ typedef typename NN::PointHandle PointHandle;
+ typedef typename NN::DistanceType DistanceType;
+ typedef typename NN::HDContainer HDContainer;
+
+ HandleDistance() {}
+ HandleDistance(PointHandle pp, DistanceType dd):
+ p(pp), d(dd) {}
+ bool operator<(const HandleDistance& other) const { return d < other.d; }
+
+ PointHandle p;
+ DistanceType d;
+};
+
+template<class HandleDistance>
+struct NNRecord
+{
+ typedef typename HandleDistance::PointHandle PointHandle;
+ typedef typename HandleDistance::DistanceType DistanceType;
+
+ NNRecord() { result.d = std::numeric_limits<DistanceType>::infinity(); }
+ DistanceType operator()(PointHandle p, DistanceType d) { if (d < result.d) { result.p = p; result.d = d; } return result.d; }
+ HandleDistance result;
+};
+
+template<class HandleDistance>
+struct rNNRecord
+{
+ typedef typename HandleDistance::PointHandle PointHandle;
+ typedef typename HandleDistance::DistanceType DistanceType;
+ typedef typename HandleDistance::HDContainer HDContainer;
+
+ rNNRecord(DistanceType r_): r(r_) {}
+ DistanceType operator()(PointHandle p, DistanceType d)
+ {
+ if (d <= r)
+ result.push_back(HandleDistance(p,d));
+ return r;
+ }
+
+ DistanceType r;
+ HDContainer result;
+};
+
+template<class HandleDistance>
+struct firstrNNRecord
+{
+ typedef typename HandleDistance::PointHandle PointHandle;
+ typedef typename HandleDistance::DistanceType DistanceType;
+ typedef typename HandleDistance::HDContainer HDContainer;
+
+ firstrNNRecord(DistanceType r_): r(r_) {}
+
+ DistanceType operator()(PointHandle p, DistanceType d)
+ {
+ if (d <= r) {
+ result.push_back(HandleDistance(p,d));
+ return -100000000.0;
+ } else {
+ return r;
+ }
+ }
+
+ DistanceType r;
+ HDContainer result;
+};
+
+
+template<class HandleDistance>
+struct kNNRecord
+{
+ typedef typename HandleDistance::PointHandle PointHandle;
+ typedef typename HandleDistance::DistanceType DistanceType;
+ typedef typename HandleDistance::HDContainer HDContainer;
+
+ kNNRecord(unsigned k_): k(k_) {}
+ DistanceType operator()(PointHandle p, DistanceType d)
+ {
+ if (result.size() < k)
+ {
+ result.push_back(HandleDistance(p,d));
+ boost::push_heap(result);
+ if (result.size() < k)
+ return std::numeric_limits<DistanceType>::infinity();
+ } else if (d < result[0].d)
+ {
+ boost::pop_heap(result);
+ result.back() = HandleDistance(p,d);
+ boost::push_heap(result);
+ }
+ if ( result.size() > 1 ) {
+ assert( result[0].d >= result[1].d );
+ }
+ return result[0].d;
+ }
+
+ unsigned k;
+ HDContainer result;
+};
+
+} // dnn
+} // bt
+} // hera
+
+#endif // HERA_BT_DNN_LOCAL_SEARCH_FUNCTORS_H
diff --git a/geom_bottleneck/include/dnn/parallel/tbb.h b/geom_bottleneck/include/dnn/parallel/tbb.h
new file mode 100644
index 0000000..14f0093
--- /dev/null
+++ b/geom_bottleneck/include/dnn/parallel/tbb.h
@@ -0,0 +1,235 @@
+#ifndef HERA_BT_PARALLEL_H
+#define HERA_BT_PARALLEL_H
+
+#ifndef FOR_R_TDA
+#include <iostream>
+#endif
+
+#include <vector>
+
+#include <boost/range.hpp>
+#include <boost/bind.hpp>
+#include <boost/foreach.hpp>
+
+#ifdef TBB
+
+#include <tbb/tbb.h>
+#include <tbb/concurrent_hash_map.h>
+#include <tbb/scalable_allocator.h>
+
+#include <boost/serialization/split_free.hpp>
+#include <boost/serialization/collections_load_imp.hpp>
+#include <boost/serialization/collections_save_imp.hpp>
+
+namespace hera {
+namespace bt {
+namespace dnn
+{
+ using tbb::mutex;
+ using tbb::task_scheduler_init;
+ using tbb::task_group;
+ using tbb::task;
+
+ template<class T>
+ struct vector
+ {
+ typedef tbb::concurrent_vector<T> type;
+ };
+
+ template<class T>
+ struct atomic
+ {
+ typedef tbb::atomic<T> type;
+ static T compare_and_swap(type& v, T n, T o) { return v.compare_and_swap(n,o); }
+ };
+
+ template<class Iterator, class F>
+ void do_foreach(Iterator begin, Iterator end, const F& f) { tbb::parallel_do(begin, end, f); }
+
+ template<class Range, class F>
+ void for_each_range_(const Range& r, const F& f)
+ {
+ for (typename Range::iterator cur = r.begin(); cur != r.end(); ++cur)
+ f(*cur);
+ }
+
+ template<class F>
+ void for_each_range(size_t from, size_t to, const F& f)
+ {
+ //static tbb::affinity_partitioner ap;
+ //tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::range_type, F>, _1, f), ap);
+ tbb::parallel_for(from, to, f);
+ }
+
+ template<class Container, class F>
+ void for_each_range(const Container& c, const F& f)
+ {
+ //static tbb::affinity_partitioner ap;
+ //tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::range_type, F>, _1, f), ap);
+ tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::const_range_type, F>, _1, f));
+ }
+
+ template<class Container, class F>
+ void for_each_range(Container& c, const F& f)
+ {
+ //static tbb::affinity_partitioner ap;
+ //tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::range_type, F>, _1, f), ap);
+ tbb::parallel_for(c.range(), boost::bind(&for_each_range_<typename Container::range_type, F>, _1, f));
+ }
+
+ template<class ID, class NodePointer, class IDTraits, class Allocator>
+ struct map_traits
+ {
+ typedef tbb::concurrent_hash_map<ID, NodePointer, IDTraits, Allocator> type;
+ typedef typename type::range_type range;
+ };
+
+ struct progress_timer
+ {
+ progress_timer(): start(tbb::tick_count::now()) {}
+ ~progress_timer()
+ {
+#ifndef FOR_R_TDA
+ std::cout << (tbb::tick_count::now() - start).seconds() << " s" << std::endl;
+#endif
+ }
+
+ tbb::tick_count start;
+ };
+}
+}
+}
+
+// Serialization for tbb::concurrent_vector<...>
+namespace boost
+{
+ namespace serialization
+ {
+ template<class Archive, class T, class A>
+ void save(Archive& ar, const tbb::concurrent_vector<T,A>& v, const unsigned int file_version)
+ { stl::save_collection(ar, v); }
+
+ template<class Archive, class T, class A>
+ void load(Archive& ar, tbb::concurrent_vector<T,A>& v, const unsigned int file_version)
+ {
+ stl::load_collection<Archive,
+ tbb::concurrent_vector<T,A>,
+ stl::archive_input_seq< Archive, tbb::concurrent_vector<T,A> >,
+ stl::reserve_imp< tbb::concurrent_vector<T,A> >
+ >(ar, v);
+ }
+
+ template<class Archive, class T, class A>
+ void serialize(Archive& ar, tbb::concurrent_vector<T,A>& v, const unsigned int file_version)
+ { split_free(ar, v, file_version); }
+
+ template<class Archive, class T>
+ void save(Archive& ar, const tbb::atomic<T>& v, const unsigned int file_version)
+ { T v_ = v; ar << v_; }
+
+ template<class Archive, class T>
+ void load(Archive& ar, tbb::atomic<T>& v, const unsigned int file_version)
+ { T v_; ar >> v_; v = v_; }
+
+ template<class Archive, class T>
+ void serialize(Archive& ar, tbb::atomic<T>& v, const unsigned int file_version)
+ { split_free(ar, v, file_version); }
+ }
+}
+
+#else
+
+#include <algorithm>
+#include <map>
+#include <boost/progress.hpp>
+
+namespace hera {
+namespace bt {
+namespace dnn
+{
+ template<class T>
+ struct vector
+ {
+ typedef ::std::vector<T> type;
+ };
+
+ template<class T>
+ struct atomic
+ {
+ typedef T type;
+ static T compare_and_swap(type& v, T n, T o) { if (v != o) return v; v = n; return o; }
+ };
+
+ template<class Iterator, class F>
+ void do_foreach(Iterator begin, Iterator end, const F& f) { std::for_each(begin, end, f); }
+
+ template<class F>
+ void for_each_range(size_t from, size_t to, const F& f)
+ {
+ for (size_t i = from; i < to; ++i)
+ f(i);
+ }
+
+ template<class Container, class F>
+ void for_each_range(Container& c, const F& f)
+ {
+ BOOST_FOREACH(const typename Container::value_type& i, c)
+ f(i);
+ }
+
+ template<class Container, class F>
+ void for_each_range(const Container& c, const F& f)
+ {
+ BOOST_FOREACH(const typename Container::value_type& i, c)
+ f(i);
+ }
+
+ struct mutex
+ {
+ struct scoped_lock
+ {
+ scoped_lock() {}
+ scoped_lock(mutex& ) {}
+ void acquire(mutex& ) const {}
+ void release() const {}
+ };
+ };
+
+ struct task_scheduler_init
+ {
+ task_scheduler_init(unsigned) {}
+ void initialize(unsigned) {}
+ static const unsigned automatic = 0;
+ static const unsigned deferred = 0;
+ };
+
+ struct task_group
+ {
+ template<class Functor>
+ void run(const Functor& f) const { f(); }
+ void wait() const {}
+ };
+
+ template<class ID, class NodePointer, class IDTraits, class Allocator>
+ struct map_traits
+ {
+ typedef std::map<ID, NodePointer,
+ typename IDTraits::Comparison,
+ Allocator> type;
+ typedef type range;
+ };
+
+ using boost::progress_timer;
+}
+}
+}
+
+#endif // TBB
+
+namespace dnn
+{
+ template<class Range, class F>
+ void do_foreach(const Range& range, const F& f) { do_foreach(boost::begin(range), boost::end(range), f); }
+}
+
+#endif
diff --git a/geom_bottleneck/include/dnn/parallel/utils.h b/geom_bottleneck/include/dnn/parallel/utils.h
new file mode 100644
index 0000000..9809e77
--- /dev/null
+++ b/geom_bottleneck/include/dnn/parallel/utils.h
@@ -0,0 +1,100 @@
+#ifndef HERA_BT_PARALLEL_UTILS_H
+#define HERA_BT_PARALLEL_UTILS_H
+
+#include "../utils.h"
+
+namespace hera
+{
+namespace bt
+{
+namespace dnn
+{
+ // Assumes rng is synchronized across ranks
+ template<class DataVector, class RNGType, class SwapFunctor>
+ void shuffle(mpi::communicator& world, DataVector& data, RNGType& rng, const SwapFunctor& swap, DataVector empty = DataVector());
+
+ template<class DataVector, class RNGType>
+ void shuffle(mpi::communicator& world, DataVector& data, RNGType& rng)
+ {
+ typedef decltype(data[0]) T;
+ shuffle(world, data, rng, [](T& x, T& y) { std::swap(x,y); });
+ }
+}
+}
+}
+
+template<class DataVector, class RNGType, class SwapFunctor>
+void
+hera::bt::dnn::shuffle(mpi::communicator& world, DataVector& data, RNGType& rng, const SwapFunctor& swap, DataVector empty)
+{
+ // This is not a perfect shuffle: it dishes out data in chunks of 1/size.
+ // (It can be interpreted as generating a bistochastic matrix by taking the
+ // sum of size random permutation matrices.) Hopefully, it works for our purposes.
+
+ typedef typename RNGType::result_type RNGResult;
+
+ int size = world.size();
+ int rank = world.rank();
+
+ // Generate local seeds
+ boost::uniform_int<RNGResult> uniform;
+ RNGResult seed;
+ for (size_t i = 0; i < size; ++i)
+ {
+ RNGResult v = uniform(rng);
+ if (i == rank)
+ seed = v;
+ }
+ RNGType local_rng(seed);
+
+ // Shuffle local data
+ hera::bt::dnn::random_shuffle(data.begin(), data.end(), local_rng, swap);
+
+ // Decide how much of our data goes to i-th processor
+ std::vector<size_t> out_counts(size);
+ std::vector<int> ranks(boost::counting_iterator<int>(0),
+ boost::counting_iterator<int>(size));
+ for (size_t i = 0; i < size; ++i)
+ {
+ hera::bt::dnn::random_shuffle(ranks.begin(), ranks.end(), rng);
+ ++out_counts[ranks[rank]];
+ }
+
+ // Fill the outgoing array
+ size_t total = 0;
+ std::vector< DataVector > outgoing(size, empty);
+ for (size_t i = 0; i < size; ++i)
+ {
+ size_t count = data.size()*out_counts[i]/size;
+ if (total + count > data.size())
+ count = data.size() - total;
+
+ outgoing[i].reserve(count);
+ for (size_t j = total; j < total + count; ++j)
+ outgoing[i].push_back(data[j]);
+
+ total += count;
+ }
+
+ boost::uniform_int<size_t> uniform_outgoing(0,size-1); // in range [0,size-1]
+ while(total < data.size()) // send leftover to random processes
+ {
+ outgoing[uniform_outgoing(local_rng)].push_back(data[total]);
+ ++total;
+ }
+ data.clear();
+
+ // Exchange the data
+ std::vector< DataVector > incoming(size, empty);
+ mpi::all_to_all(world, outgoing, incoming);
+ outgoing.clear();
+
+ // Assemble our data
+ for(const DataVector& vec : incoming)
+ for (size_t i = 0; i < vec.size(); ++i)
+ data.push_back(vec[i]);
+ hera::bt::dnn::random_shuffle(data.begin(), data.end(), local_rng, swap);
+ // XXX: the final shuffle is irrelevant for our purposes. But it's also cheap.
+}
+
+#endif
diff --git a/geom_bottleneck/include/dnn/utils.h b/geom_bottleneck/include/dnn/utils.h
new file mode 100644
index 0000000..f4ce632
--- /dev/null
+++ b/geom_bottleneck/include/dnn/utils.h
@@ -0,0 +1,47 @@
+#ifndef HERA_BT_DNN_UTILS_H
+#define HERA_BT_DNN_UTILS_H
+
+#include <boost/random/uniform_int.hpp>
+#include <boost/foreach.hpp>
+#include <boost/typeof/typeof.hpp>
+
+namespace hera
+{
+namespace bt
+{
+namespace dnn
+{
+
+template <typename T, typename... Args>
+struct has_coordinates
+{
+ template <typename C, typename = decltype( std::declval<C>().coordinate(std::declval<Args>()...) )>
+ static std::true_type test(int);
+
+ template <typename C>
+ static std::false_type test(...);
+
+ static constexpr bool value = decltype(test<T>(0))::value;
+};
+
+template<class RandomIt, class UniformRandomNumberGenerator, class SwapFunctor>
+void random_shuffle(RandomIt first, RandomIt last, UniformRandomNumberGenerator& g, const SwapFunctor& swap)
+{
+ size_t n = last - first;
+ boost::uniform_int<size_t> uniform(0,n);
+ for (size_t i = n-1; i > 0; --i)
+ swap(first[i], first[uniform(g,i+1)]); // picks a random number in [0,i] range
+}
+
+template<class RandomIt, class UniformRandomNumberGenerator>
+void random_shuffle(RandomIt first, RandomIt last, UniformRandomNumberGenerator& g)
+{
+ typedef decltype(*first) T;
+ random_shuffle(first, last, g, [](T& x, T& y) { std::swap(x,y); });
+}
+
+} // dnn
+} // bt
+} // hera
+
+#endif
diff --git a/geom_bottleneck/include/neighb_oracle.h b/geom_bottleneck/include/neighb_oracle.h
new file mode 100644
index 0000000..c3751b3
--- /dev/null
+++ b/geom_bottleneck/include/neighb_oracle.h
@@ -0,0 +1,295 @@
+/*
+
+Copyright (c) 2015, M. Kerber, D. Morozov, A. Nigmetov
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
+
+2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
+
+3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+
+You are under no obligation whatsoever to provide any bug fixes, patches, or
+upgrades to the features, functionality or performance of the source code
+(Enhancements) to anyone; however, if you choose to make your Enhancements
+available either publicly, or directly to copyright holder,
+without imposing a separate written license agreement for such Enhancements,
+then you hereby grant the following license: a non-exclusive, royalty-free
+perpetual license to install, use, modify, prepare derivative works, incorporate
+into other computer software, distribute, and sublicense such enhancements or
+derivative works thereof, in binary and source code form.
+
+*/
+
+#ifndef HERA_NEIGHB_ORACLE_H
+#define HERA_NEIGHB_ORACLE_H
+
+#include <unordered_map>
+#include <algorithm>
+#include <memory>
+
+#include "basic_defs_bt.h"
+#include "dnn/geometry/euclidean-fixed.h"
+#include "dnn/local/kd-tree.h"
+
+
+
+namespace hera {
+namespace bt {
+
+template<class Real>
+class NeighbOracleSimple
+{
+public:
+ using DgmPoint = DiagramPoint<Real>;
+ using DgmPointSet = DiagramPointSet<Real>;
+
+private:
+ Real r;
+ Real distEpsilon;
+ DgmPointSet pointSet;
+
+public:
+
+ NeighbOracleSimple() : r(0.0) {}
+
+ NeighbOracleSimple(const DgmPointSet& _pointSet, const Real _r, const Real _distEpsilon) :
+ r(_r),
+ distEpsilon(_distEpsilon),
+ pointSet(_pointSet)
+ {}
+
+ void deletePoint(const DgmPoint& p)
+ {
+ pointSet.erase(p);
+ }
+
+ void rebuild(const DgmPointSet& S, const double rr)
+ {
+ pointSet = S;
+ r = rr;
+ }
+
+ bool getNeighbour(const DgmPoint& q, DgmPoint& result) const
+ {
+ for(auto pit = pointSet.cbegin(); pit != pointSet.cend(); ++pit) {
+ if ( distLInf(*pit, q) <= r) {
+ result = *pit;
+ return true;
+ }
+ }
+ return false;
+ }
+
+ void getAllNeighbours(const DgmPoint& q, std::vector<DgmPoint>& result)
+ {
+ result.clear();
+ for(const auto& point : pointSet) {
+ if ( distLInf(point, q) <= r) {
+ result.push_back(point);
+ }
+ }
+ for(auto& pt : result) {
+ deletePoint(pt);
+ }
+ }
+
+};
+
+template<class Real_>
+class NeighbOracleDnn
+{
+public:
+
+ using Real = Real_;
+ using DnnPoint = dnn::Point<2, double>;
+ using DnnTraits = dnn::PointTraits<DnnPoint>;
+ using DgmPoint = DiagramPoint<Real>;
+ using DgmPointSet = DiagramPointSet<Real>;
+ using DgmPointHash = DiagramPointHash<Real>;
+
+ Real r;
+ Real distEpsilon;
+ std::vector<DgmPoint> allPoints;
+ DgmPointSet diagonalPoints;
+ std::unordered_map<DgmPoint, size_t, DgmPointHash> pointIdxLookup;
+ // dnn-stuff
+ std::unique_ptr<dnn::KDTree<DnnTraits>> kdtree;
+ std::vector<DnnPoint> dnnPoints;
+ std::vector<DnnPoint*> dnnPointHandles;
+ std::vector<size_t> kdtreeItems;
+
+ NeighbOracleDnn(const DgmPointSet& S, const Real rr, const Real dEps) :
+ kdtree(nullptr)
+ {
+ assert(dEps >= 0);
+ distEpsilon = dEps;
+ rebuild(S, rr);
+ }
+
+
+ void deletePoint(const DgmPoint& p)
+ {
+ auto findRes = pointIdxLookup.find(p);
+ assert(findRes != pointIdxLookup.end());
+ //std::cout << "Deleting point " << p << std::endl;
+ size_t pointIdx { (*findRes).second };
+ //std::cout << "pointIdx = " << pointIdx << std::endl;
+ diagonalPoints.erase(p, false);
+ kdtree->delete_point(dnnPointHandles[kdtreeItems[pointIdx]]);
+ }
+
+ void rebuild(const DgmPointSet& S, const Real rr)
+ {
+ //std::cout << "Entered rebuild, r = " << rr << std::endl;
+ r = rr;
+ size_t dnnNumPoints = S.size();
+ //printDebug(isDebug, "S = ", S);
+ if (dnnNumPoints > 0) {
+ pointIdxLookup.clear();
+ pointIdxLookup.reserve(S.size());
+ allPoints.clear();
+ allPoints.reserve(S.size());
+ diagonalPoints.clear();
+ diagonalPoints.reserve(S.size() / 2);
+ for(auto pit = S.cbegin(); pit != S.cend(); ++pit) {
+ allPoints.push_back(*pit);
+ if (pit->isDiagonal()) {
+ diagonalPoints.insert(*pit);
+ }
+ }
+
+ size_t pointIdx = 0;
+ for(auto& dataPoint : allPoints) {
+ pointIdxLookup.insert( { dataPoint, pointIdx++ } );
+ }
+
+ size_t dnnItemIdx { 0 };
+ size_t trueIdx { 0 };
+ dnnPoints.clear();
+ kdtreeItems.clear();
+ dnnPointHandles.clear();
+ dnnPoints.clear();
+ kdtreeItems.reserve(S.size() );
+ // store normal items in kd-tree
+ for(const auto& g : allPoints) {
+ if (true) {
+ kdtreeItems[trueIdx] = dnnItemIdx;
+ // index of items is id of dnn-point
+ DnnPoint p(trueIdx);
+ p[0] = g.getRealX();
+ p[1] = g.getRealY();
+ dnnPoints.push_back(p);
+ assert(dnnItemIdx == dnnPoints.size() - 1);
+ dnnItemIdx++;
+ }
+ trueIdx++;
+ }
+ assert(dnnPoints.size() == allPoints.size() );
+ for(size_t i = 0; i < dnnPoints.size(); ++i) {
+ dnnPointHandles.push_back(&dnnPoints[i]);
+ }
+ DnnTraits traits;
+ //std::cout << "kdtree: " << dnnPointHandles.size() << " points" << std::endl;
+ kdtree.reset(new dnn::KDTree<DnnTraits>(traits, dnnPointHandles));
+ }
+ }
+
+
+ bool getNeighbour(const DgmPoint& q, DgmPoint& result) const
+ {
+ //std::cout << "getNeighbour for q = " << q << ", r = " << r << std::endl;
+ //std::cout << *this << std::endl;
+ // distance between two diagonal points
+ // is 0
+ if (q.isDiagonal()) {
+ if (!diagonalPoints.empty()) {
+ result = *diagonalPoints.cbegin();
+ //std::cout << "Neighbour found in diagonal points, res = " << result;
+ return true;
+ }
+ }
+ // check if kdtree is not empty
+ if (0 == kdtree->get_num_points() ) {
+ //std::cout << "empty tree, no neighb." << std::endl;
+ return false;
+ }
+ // if no neighbour found among diagonal points,
+ // search in kd_tree
+ DnnPoint queryPoint;
+ queryPoint[0] = q.getRealX();
+ queryPoint[1] = q.getRealY();
+ auto kdtreeResult = kdtree->findFirstR(queryPoint, r);
+ if (kdtreeResult.empty()) {
+ //std::cout << "no neighbour within " << r << "found." << std::endl;
+ return false;
+ }
+ if (kdtreeResult[0].d <= r + distEpsilon) {
+ result = allPoints[kdtreeResult[0].p->id()];
+ //std::cout << "Neighbour found with kd-tree, index = " << kdtreeResult[0].p->id() << std::endl;
+ //std::cout << "result = " << result << std::endl;
+ return true;
+ }
+ //std::cout << "No neighbour found for r = " << r << std::endl;
+ return false;
+ }
+
+
+
+ void getAllNeighbours(const DgmPoint& q, std::vector<DgmPoint>& result)
+ {
+ //std::cout << "Entered getAllNeighbours for q = " << q << std::endl;
+ result.clear();
+ // add diagonal points, if necessary
+ if ( q.isDiagonal() ) {
+ for( auto& diagPt : diagonalPoints ) {
+ result.push_back(diagPt);
+ }
+ }
+ // delete diagonal points we found
+ // to prevent finding them again
+ for(auto& pt : result) {
+ //std::cout << "deleting DIAG point pt = " << pt << std::endl;
+ deletePoint(pt);
+ }
+ size_t diagOffset = result.size();
+ std::vector<size_t> pointIndicesOut;
+ // perorm range search on kd-tree
+ DnnPoint queryPoint;
+ queryPoint[0] = q.getRealX();
+ queryPoint[1] = q.getRealY();
+ auto kdtreeResult = kdtree->findR(queryPoint, r);
+ pointIndicesOut.reserve(kdtreeResult.size());
+ for(auto& handleDist : kdtreeResult) {
+ if (handleDist.d <= r + distEpsilon) {
+ pointIndicesOut.push_back(handleDist.p->id());
+ } else {
+ break;
+ }
+ }
+ // get actual points in result
+ for(auto& ptIdx : pointIndicesOut) {
+ result.push_back(allPoints[ptIdx]);
+ }
+ // delete all points we found
+ for(auto ptIt = result.begin() + diagOffset; ptIt != result.end(); ++ptIt) {
+ //printDebug(isDebug, "deleting point pt = ", *ptIt);
+ deletePoint(*ptIt);
+ }
+ }
+
+ //DgmPointSet originalPointSet;
+ template<class R>
+ friend std::ostream& operator<<(std::ostream& out, const NeighbOracleDnn<R>& oracle);
+
+};
+
+} // end namespace bt
+} // end namespace hera
+
+#endif // HERA_NEIGHB_ORACLE_H
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-01 16:07:25 +0000