From 33bf675ba1bf103e6dbc8a50bba0d5cc8d6f8d48 Mon Sep 17 00:00:00 2001
From: Marcel Jacobse <44684927+mjacobse@users.noreply.github.com>
Date: Wed, 6 Dec 2023 18:04:13 +0100
Subject: [PATCH] Fix scaling-dependence of termination criterion

The termination criterion for bisecting alpha implicitly assumes the
optimal alpha to be in the order of magnitude of 1.0. For larger values
of alpha, lower and upper can never get close enough, meaning the
bisection never terminates successfully. For smaller values of alpha,
the bisection may terminate much too early.

This fixes that by checking the relative instead of the absolute
difference using numpy.isclose. Also, this adds tests with a large and
a tiny optimal alpha that would have failed with the previous
termination criterion (they also wouldve failed for much more moderate
scaling factors like 100).
---
 alphashape/optimizealpha.py |  2 +-
 tests/test_optimizealpha.py | 23 +++++++++++++++++++++++
 2 files changed, 24 insertions(+), 1 deletion(-)

diff --git a/alphashape/optimizealpha.py b/alphashape/optimizealpha.py
index b228c86..aa6d21e 100644
--- a/alphashape/optimizealpha.py
+++ b/alphashape/optimizealpha.py
@@ -97,7 +97,7 @@ def optimizealpha(points: Union[List[Tuple[float]], np.ndarray],
 
     # Begin the bisection loop
     counter = 0
-    while (upper - lower) > np.finfo(float).eps * 2:
+    while not np.isclose(lower, upper, atol=0.0, rtol=np.finfo(float).eps * 2):
         # Bisect the current bounds
         test_alpha = (upper + lower) * .5
 
diff --git a/tests/test_optimizealpha.py b/tests/test_optimizealpha.py
index 6609cad..43bf25c 100644
--- a/tests/test_optimizealpha.py
+++ b/tests/test_optimizealpha.py
@@ -5,6 +5,7 @@
 
 
 import unittest
+import numpy as np
 
 from alphashape import optimizealpha
 
@@ -39,3 +40,25 @@ def test_reach_max_iterations(self):
                  (0.5, 0.25), (0.5, 0.75), (0.25, 0.5), (0.75, 0.5)],
                 max_iterations=2, lower=0.0, upper=1000.0)
         self.assertEqual(alpha, 0.0)
+
+    def test_large_alpha(self):
+        """
+        Given a polygon for which a large alpha is optimal, the optimizealpha
+        function should find a large alpha
+        """
+        scale = 1e30
+        alpha = optimizealpha(np.array(
+            [(0., 0.), (0., 1.), (1., 1.), (1., 0.),
+             (0.5, 0.25), (0.5, 0.75), (0.25, 0.5), (0.75, 0.5)]) / scale)
+        assert alpha > 3. * scale and alpha < 3.5 * scale
+
+    def test_tiny_alpha(self):
+        """
+        Given a polygon for which a tiny alpha is optimal, the optimizealpha
+        function should find a tiny alpha
+        """
+        scale = 1e-30
+        alpha = optimizealpha(np.array(
+            [(0., 0.), (0., 1.), (1., 1.), (1., 0.),
+             (0.5, 0.25), (0.5, 0.75), (0.25, 0.5), (0.75, 0.5)]) / scale)
+        assert alpha > 3. * scale and alpha < 3.5 * scale