temporary fix for exact version from files

3 files changed, 75 insertions, 10 deletions

diff --git a/geom_bottleneck/bottleneck/include/bottleneck.h b/geom_bottleneck/bottleneck/include/bottleneck.h
index 3267c5d..2078965 100644
--- a/geom_bottleneck/bottleneck/include/bottleneck.h
+++ b/geom_bottleneck/bottleneck/include/bottleneck.h
@@ -51,8 +51,12 @@ std::pair<double, double> bottleneckDistApproxInterval(DiagramPointSet& A, Diagr
 double bottleneckDistApprox(DiagramPointSet& A, DiagramPointSet& B, const double epsilon);
 
 // get exact bottleneck distance,
+double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B, const int decPrecision);
+
+// get exact bottleneck distance,
 double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B);
 
+
 // 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)
@@ -86,12 +90,31 @@ double bottleneckDistExact(PairContainer& A, PairContainer& B)
 {
     DiagramPointSet a(A.begin(), A.end());
     DiagramPointSet b(B.begin(), B.end());
-    return bottleneckDistExact(a, b);
+    return bottleneckDistExact(a, b, 14);
+}
+
+
+// get exact bottleneck distance,
+template<class PairContainer>
+double bottleneckDistExact(PairContainer& A, PairContainer& B, const int decPrecision)
+{
+    DiagramPointSet a(A.begin(), A.end());
+    DiagramPointSet b(B.begin(), B.end());
+    return bottleneckDistExact(a, b, decPrecision);
 }
 
 // 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
+bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double>>& result, int& decPrecision);
+// wrapper for standard string
+bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result, int& decPrecision);
+
+
+// these two functions are now just wrappers for the previous ones,
+// in case someone needs them; decPrecision is ignored
 bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double>>& result);
 // wrapper for standard string
 bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result);
diff --git a/geom_bottleneck/bottleneck/src/bottleneck.cpp b/geom_bottleneck/bottleneck/src/bottleneck.cpp
index 34f7ed4..92e9dbf 100644
--- a/geom_bottleneck/bottleneck/src/bottleneck.cpp
+++ b/geom_bottleneck/bottleneck/src/bottleneck.cpp
@@ -100,7 +100,7 @@ double bottleneckDistApprox(DiagramPointSet& A, DiagramPointSet& B, const double
 }
 
 
-double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B, std::vector<double>& pairwiseDist)
+double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B, std::vector<double>& pairwiseDist, const int decPrecision)
 {
     //for(size_t k = 0; k < pairwiseDist.size(); ++k) {
         //std::cout << "pairwiseDist[" << k << "] = " << std::setprecision(15) << pairwiseDist[k] << std::endl;
@@ -111,13 +111,19 @@ double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B
 
     bool useRangeSearch = true;
     double distEpsilon = std::numeric_limits<double>::max();
+    double 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];
-        if ( diff > 1.0e-14 and diff < distEpsilon ) {
+        //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 oracle(A, B, distEpsilon, useRangeSearch);
     // binary search
@@ -145,6 +151,11 @@ double bottleneckDistExactFromSortedPwDist(DiagramPointSet&A, DiagramPointSet& B
 
 double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B)
 {
+    return bottleneckDistExact(A, B, 14);
+}
+
+double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B, const int decPrecision)
+{
     constexpr double epsilon = 0.001;
     auto interval = bottleneckDistApproxInterval(A, B, epsilon);
     const double delta = 0.5 * (interval.second - interval.first);
@@ -291,7 +302,7 @@ double bottleneckDistExact(DiagramPointSet& A, DiagramPointSet& B)
         //std::cout << std::setprecision(15) << ddd << std::endl;
     //}
 
-    return bottleneckDistExactFromSortedPwDist(A, B, pairwiseDist);
+    return bottleneckDistExactFromSortedPwDist(A, B, pairwiseDist, decPrecision);
 }
 
 double bottleneckDistSlow(DiagramPointSet& A, DiagramPointSet& B)
@@ -495,13 +506,27 @@ double bottleneckDistSlow(DiagramPointSet& A, DiagramPointSet& B)
     */
 }
 
+// wrappers
 bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result)
 {
-    return readDiagramPointSet(fname.c_str(), result);
+    int decPrecision;
+    return readDiagramPointSet(fname.c_str(), result, decPrecision);
 }
 
 bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double>>& result)
 {
+    int decPrecision;
+    return readDiagramPointSet(fname, result, decPrecision);
+}
+
+bool readDiagramPointSet(const std::string& fname, std::vector<std::pair<double, double>>& result, int& decPrecision)
+{
+    return readDiagramPointSet(fname.c_str(), result, decPrecision);
+}
+
+// reading function
+bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double>>& result, int& decPrecision)
+{
     size_t lineNumber { 0 };
     result.clear();
     std::ifstream f(fname);
@@ -530,6 +555,22 @@ bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double
             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;
+                }
+            } 
+        }
         double x, y;
         std::istringstream iss(line);
         if (not(iss >> x >> y)) {
@@ -542,7 +583,7 @@ bool readDiagramPointSet(const char* fname, std::vector<std::pair<double, double
             result.push_back(std::make_pair(x,y));
         } else {
 #ifndef FOR_R_TDA
-            std::cerr << "Warning: in file " << fname << ", line number " << lineNumber << ", zero persistence point: \"" << line << "\"" << std::endl;
+            //std::cerr << "Warning: in file " << fname << ", line number " << lineNumber << ", zero persistence point ignored: \"" << line << "\"" << std::endl;
 #endif
         }
     }
diff --git a/geom_bottleneck/example/bottleneck_dist.cpp b/geom_bottleneck/example/bottleneck_dist.cpp
index 9a50ce2..655de44 100644
--- a/geom_bottleneck/example/bottleneck_dist.cpp
+++ b/geom_bottleneck/example/bottleneck_dist.cpp
@@ -14,11 +14,12 @@ int main(int argc, char* argv[])
     }
 
     PairVector diagramA, diagramB;
-    if (!geom_bt::readDiagramPointSet(argv[1], diagramA)) {
+    int decPrecision { 0 };
+    if (!geom_bt::readDiagramPointSet(argv[1], diagramA, decPrecision)) {
         std::exit(1);
     }
 
-    if (!geom_bt::readDiagramPointSet(argv[2], diagramB)) {
+    if (!geom_bt::readDiagramPointSet(argv[2], diagramB, decPrecision)) {
         std::exit(1);
     }
 
@@ -30,14 +31,14 @@ int main(int argc, char* argv[])
         if (approxEpsilon > 0.0) {
             res = geom_bt::bottleneckDistApprox(diagramA, diagramB, approxEpsilon);
         } else if (approxEpsilon == 0.0) {
-            res = geom_bt::bottleneckDistExact(diagramA, diagramB);
+            res = geom_bt::bottleneckDistExact(diagramA, diagramB, decPrecision);
         } else {
             std::cerr << "The third parameter (relative error) must be positive!" << std::endl;
             std::exit(1);
         }
     } else {
         // only filenames have been supplied, return exact distance
-        res = geom_bt::bottleneckDistExact(diagramA, diagramB);
+        res = geom_bt::bottleneckDistExact(diagramA, diagramB, decPrecision);
     }
     std::cout << std::setprecision(15) << res << std::endl;