important numeric fixes

Author: derohde

README.md

@@ -98,7 +98,7 @@ All simplifications are vertex-restricted!
 - signature: `fred.Cluster_Assignment`
 - methods: 
     - `len(fred.Cluster_Assignment)`: number of centers
-    -`fred.Cluster_Assignment.count(i)`: number of curves assigned to center `i`
+    - `fred.Cluster_Assignment.count(i)`: number of curves assigned to center `i`
     - `fred.Cluster_Assignment.get(i,j)`: get index of `j`th curve assigned to center `i`
 
 ### Dimension Reduction via Gaussian Random Projection 

include/point.hpp

@@ -165,8 +165,8 @@ class Point : public Coordinates {
         const Vector u = line_end-line_start, v = *this - line_start;
         const parameter_t ulen_sqr = u.length_sqr(), vlen_sqr = v.length_sqr();
 
-        if (ulen_sqr == 0) {
-            if (vlen_sqr <= distance_sqr) return Interval(0, 1);
+        if (near_eq(ulen_sqr, parameter_t(0))) {
+            if (vlen_sqr <= distance_sqr) return Interval(parameter_t(0), parameter_t(1));
             else return Interval();
         }
 

include/types.hpp

@@ -11,6 +11,8 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI
 #pragma once
 
 #include <vector>
+#include <cmath>
+#include <limits>
 
 typedef double distance_t;
 typedef double coordinate_t;
@@ -29,3 +31,9 @@ using Vector = Point;
 using Intervals = std::vector<Interval>;
 using Coordinates = std::vector<coordinate_t>;
 
+template<typename T>
+inline bool near_eq(T x, T y) {
+  return std::abs(x - y) <= std::min(std::abs(x), std::abs(y)) * std::numeric_limits<T>::epsilon();
+}
+
+

setup.py

@@ -74,7 +74,7 @@ def build_extension(self, ext):
 
 setup(
     name='Fred-Frechet',
-    version='1.9.3',
+    version='1.9.4',
     author='Dennis Rohde',
     author_email='[email protected]',
     description='A fast, scalable and light-weight C++ Fréchet distance library, exposed to python and focused on (k,l)-clustering of polygonal curves.',

src/frechet.cpp

@@ -80,7 +80,8 @@ Distance _distance(const Curve &curve1, const Curve &curve2, distance_t ub, dist
         //Binary search over the feasible distances
         while (ub - lb > p_error) {
             ++number_searches;
-            split = (ub + lb)/2;
+            split = (ub + lb)/distance_t(2);
+            if (split == lb or split == ub) break;
             auto isLessThan = _less_than_or_equal(split, curve1, curve2, reachable1, reachable2, free_intervals1, free_intervals2);
             if (isLessThan) {
                 ub = split;
@@ -93,7 +94,7 @@ Distance _distance(const Curve &curve1, const Curve &curve2, distance_t ub, dist
     }
 
     const auto end = std::chrono::high_resolution_clock::now();
-    result.value = (ub + lb)/2.;
+    result.value = lb;
     result.time_searches = std::chrono::duration_cast<std::chrono::seconds>(end - start).count();
     result.number_searches = number_searches;
     return result;

src/fred_python_wrapper.cpp

@@ -50,8 +50,6 @@ Curve approximate_minimum_error_simplification(const Curve &curve, const curve_s
 
 PYBIND11_MODULE(backend, m) {
 
-    m.def("set_maximum_number_threads", set_number_threads);
-    
     py::class_<Config::Config>(m, "Config")
         .def(py::init<>())
         .def_property("continuous_frechet_error", [&](Config::Config&) { return fc::error; }, [&](Config::Config&, const bool error) { fc::error = error; })

src/simplification.cpp

@@ -88,9 +88,10 @@ Curve Simplification::approximate_minimum_error_simplification(const Curve &curv
     }
 
     if (Config::verbose) std::cout << "ASIMPL: binary search using upper bound" << std::endl;
-    while (max_distance - min_distance > min_distance * Frechet::Continuous::error / 100) {
-        mid_distance = (min_distance + max_distance) / 2.;
-        
+    const distance_t epsilon = std::max(min_distance * Frechet::Continuous::error / 100, std::numeric_limits<distance_t>::epsilon());
+    while (max_distance - min_distance > epsilon) {
+        mid_distance = (min_distance + max_distance) / distance_t(2);
+        if (mid_distance == max_distance or mid_distance == min_distance) break;
         new_simplification = Simplification::approximate_minimum_link_simplification(curve, mid_distance);
 
         if (new_simplification.complexity() > ell) min_distance = mid_distance;