Fix bug with exact bottleneck distance.

3 files changed, 19 insertions, 20 deletions

diff --git a/bottleneck/include/bottleneck_detail.hpp b/bottleneck/include/bottleneck_detail.hpp
index 8f51d07..f976909 100644
--- a/bottleneck/include/bottleneck_detail.hpp
+++ b/bottleneck/include/bottleneck_detail.hpp
@@ -51,7 +51,6 @@ namespace hera {
         void binarySearch(const Real epsilon,
                           std::pair<Real, Real>& result,
                           BoundMatchOracle <Real>& oracle,
-                          const Real infinityCost,
                           bool isResultInitializedCorrectly,
                           const Real distProbeInit)
         {
@@ -59,9 +58,6 @@ namespace hera {
             Real& distMin = result.first;
             Real& distMax = result.second;
 
-            distMin = std::max(distMin, infinityCost);
-            distMax = std::max(distMax, infinityCost);
-
             Real distProbe;
 
             if (not isResultInitializedCorrectly) {
@@ -87,13 +83,6 @@ namespace hera {
             // bounds are correct , perform binary search
             distProbe = (distMin + distMax) / 2.0;
             while ((distMax - distMin) / distMin >= epsilon) {
-
-                if (distMax < infinityCost) {
-                    distMin = infinityCost;
-                    distMax = infinityCost;
-                    break;
-                }
-
                 if (oracle.isMatchLess(distProbe)) {
                     distMax = distProbe;
                 } else {
@@ -102,9 +91,6 @@ namespace hera {
 
                 distProbe = (distMin + distMax) / 2.0;
             }
-
-            distMin = std::max(distMin, infinityCost);
-            distMax = std::max(distMax, infinityCost);
         }
 
         //        template<class Real>
@@ -291,7 +277,7 @@ namespace hera {
             // 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) };
-            binarySearch(epsilon, result, oracle, infinity_cost, false, distProbe);
+            binarySearch(epsilon, result, oracle, false, distProbe);
             // to compute longest edge a perfect matching is needed
             if (compute_longest_edge and result.first > infinity_cost) {
                 oracle.isMatchLess(result.second);
@@ -448,8 +434,10 @@ namespace hera {
         Real bottleneckDistApprox(DiagramPointSet <Real>& A, DiagramPointSet <Real>& B, const Real epsilon,
                                   MatchingEdge <Real>& longest_edge, bool compute_longest_edge)
         {
+            // must compute here: infinity points will be erased in bottleneckDistApproxInterval
+            Real infCost = getInfinityCost(A, B).cost;
             auto interval = bottleneckDistApproxInterval<Real>(A, B, epsilon, longest_edge, compute_longest_edge);
-            return interval.second;
+            return std::max(infCost, interval.second);
         }
 
 
@@ -518,6 +506,9 @@ namespace hera {
             using DgmPoint = DiagramPoint<Real>;
 
             constexpr Real epsilon = 0.001;
+
+            Real infCost = getInfinityCost(A, B, true).cost;
+
             auto interval = bottleneckDistApproxInterval(A, B, epsilon, longest_edge, true);
             // if the longest edge is on infinity, the answer is already exact
             // this will be detected here and all the code after if
@@ -627,8 +618,9 @@ namespace hera {
                 pw_dists.push_back(d);
             }
 
-            return bottleneckDistExactFromSortedPwDist(A, B, pw_dists, decPrecision, longest_edge,
-                                                       compute_longest_edge);
+            Real exactFinite = bottleneckDistExactFromSortedPwDist(A, B, pw_dists, decPrecision, longest_edge, compute_longest_edge);
+
+            return std::max(infCost, exactFinite);
         }
 
     } // end namespace bt
diff --git a/bottleneck/tests/data/test_list.txt b/bottleneck/tests/data/test_list.txt
index 24b462a..6fc63cc 100644
--- a/bottleneck/tests/data/test_list.txt
+++ b/bottleneck/tests/data/test_list.txt
@@ -790,3 +790,4 @@ test_876_A test_876_B 0.01 165
 test_877_A test_877_B 0.01 190
 test_878_A test_878_B 0.01 242
 test_879_A test_879_B 0.01 187
+test_880_A test_880_B 0.000001 0.0411428571428574
diff --git a/bottleneck/tests/test_hera_bottleneck.cpp b/bottleneck/tests/test_hera_bottleneck.cpp
index f22e415..b31b2b7 100644
--- a/bottleneck/tests/test_hera_bottleneck.cpp
+++ b/bottleneck/tests/test_hera_bottleneck.cpp
@@ -159,6 +159,9 @@ TEST_CASE("infinity points", "bottleneckDistApprox")
         double d = hera::bottleneckDistApprox<>(diagram_A, diagram_B, delta);
         double corr_answer = 1.0;
         REQUIRE(  fabs(d - corr_answer) <= delta * corr_answer);
+
+        double exact_d = hera::bottleneckDistExact<>(diagram_A, diagram_B);
+        REQUIRE(exact_d == corr_answer);
     }
 
     SECTION("two points at infinity") {
@@ -565,7 +568,11 @@ TEST_CASE("file cases", "bottleneck_dist")
             REQUIRE(read_file_A);
             REQUIRE(read_file_B);
 
-            double hera_answer = hera::bottleneckDistApprox(diagram_A, diagram_B, ts.delta, longest_edge, true);
+            double hera_answer;
+            if (ts.delta > 0)
+                hera_answer = hera::bottleneckDistApprox(diagram_A, diagram_B, ts.delta, longest_edge, true);
+            else
+                hera_answer = hera::bottleneckDistExact(diagram_A, diagram_B);
             std::pair<int, int> hera_le { longest_edge.first.get_user_id(), longest_edge.second.get_user_id() };
 
             REQUIRE((hera_answer == ts.answer or fabs(hera_answer - ts.answer) <= ts.delta * hera_answer));
@@ -604,7 +611,6 @@ TEST_CASE("file cases", "bottleneck_dist")
 //                          << hera_answer_exact << std::endl;
 //            }
             REQUIRE( (not check_longest_edge_cost or fabs(hera_le_cost - hera_answer_exact) < 0.0001 * hera_answer_exact) );
-            std::cout << ts << " PASSED " << std::endl;
         }
     }