Numeric Stability and Readme

Author: derohde

README.md

@@ -29,13 +29,29 @@ By default, Fred will automatically determine the number of threads to use. If y
 - signature: `fred.discrete_frechet(curve1, curve2)`
 - returns: `fred.Discrete_Frechet_Result` with members `value` and `time`
 
+### Curve Simplification
+
+#### weak minimum error simplification
+- graph approach from [**Polygonal Approximations of a Curve — Formulations and Algorithms**](https://www.sciencedirect.com/science/article/pii/B9780444704672500114)
+- signature: `fred.weak_minimum_error_simplification(fred.Curve, int complexity)`
+- returns: `fred.Curve`that uses input curves vertices, with `complexity` number of vertices and that has minimum distance to input curve
+
+#### approximate weak minimum link simplification
+- algorithm "FS" from [**Near-Linear Time Approximation Algorithms for Curve Simplification**](https://link.springer.com/article/10.1007/s00453-005-1165-y)
+- signature: `fred.approximate_weak_minimum_error_simplification(fred.Curve, double error)`
+- returns: `fred.Curve` that uses input curves vertices, is of small complexity and with distance to input curve at most `error`
+
+#### approximate weak minimum error simplification
+- binary search on `fred.approximate_weak_minimum_link_simplification`
+- signature: `fred.approximate_weak_minimum_error_simplification(fred.Curve, int complexity)`
+- returns: `fred.Curve`that uses input curves vertices, with `complexity` number of vertices and that has small distance to input curve
+
 ### Clustering
 
 ##### Distance_Matrix
 
 A `fred.Distance_Matrix()` can be used to speed up consecutive calls of `fred.discrete_klcenter` and `fred.discrete_klmedian`. As the name suggests, it stores the distances already computed.
 
-
 #### discrete (k,l)-center clustering (continuous Fréchet) -- multiple calls
 - from [**Approximating (k,l)-center clustering for curves**](https://dl.acm.org/doi/10.5555/3310435.3310616)
 - signature: `fred.discrete_klcenter_multi(k, l, curves, distances, center_domain, random_first_center)` with parameters 

experiments/dtw_mislead_i.py

@@ -125,7 +125,7 @@ class Curve : private Points {
         np::dtype dt = np::dtype::get_builtin<coordinate_t>();
         p::list l;
         np::ndarray result = np::array(l, dt);
-        for (const auto &elem : *this) {
+        for (const Point &elem : *this) {
             l.append(elem.as_ndarray());
         }
         result = np::array(l, dt);

experiments/dtw_mislead_m.py

@@ -20,7 +20,7 @@ class Interval {
     parameter_t beg, en;
 
 public:
-    Interval() : beg{1}, en{0} {}
+    Interval() : beg{1.L}, en{0.L} {}
 
     Interval(const parameter_t begin, const parameter_t end) : beg{begin}, en{end} {}
 
@@ -29,7 +29,8 @@ class Interval {
     }
 
     inline bool is_empty() const { 
-        return beg > en; 
+        if (en - beg >= std::numeric_limits<parameter_t>::epsilon()) return beg > en;
+        else return true;
     }
 
     inline bool intersects(const Interval &other) const {
@@ -49,8 +50,8 @@ class Interval {
     }
 
     inline void reset() {
-        beg = 1;
-        en = 0;
+        beg = 1.L;
+        en = 0.L;
     }
 };
 

experiments/dtw_normal.py

@@ -72,7 +72,7 @@ class Point : public Coordinates {
     }
 
     inline Point operator+(const Point &point) const {
-        auto result = *this;
+        Point result = *this;
         #pragma omp simd
         for (dimensions_t i = 0; i < dimensions(); ++i){
             result.operator[](i) += point[i];
@@ -81,7 +81,7 @@ class Point : public Coordinates {
     }
 
     inline Point operator-(const Point &point) const {
-        auto result = *this;
+        Point result = *this;
         #pragma omp simd
         for (dimensions_t i = 0; i < dimensions(); ++i){
             result.operator[](i) -= point[i];
@@ -145,32 +145,32 @@ class Point : public Coordinates {
 
     inline Interval intersection_interval(const distance_t distance_sqr, const Point &line_start, const Point &line_end) const {
         const Vector u = line_end-line_start, v = *this - line_start;
-        const distance_t ulen_sqr = u.length_sqr(), vlen_sqr = v.length_sqr();
+        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);
             else return Interval();
         }
 
-        const distance_t p =  -2. / ulen_sqr * (u * v), q = (vlen_sqr - distance_sqr) / ulen_sqr;
+        const parameter_t p =  -2. * ((u * v) / ulen_sqr), q = vlen_sqr / ulen_sqr - distance_sqr / ulen_sqr;
 
-        const distance_t phalf_sqr = p * p / 4., discriminant = phalf_sqr - q;
+        const parameter_t phalf_sqr = p * p / 4., discriminant = phalf_sqr - q;
 
         if (discriminant < 0) return Interval();
 
-        const distance_t discriminant_sqrt = std::sqrt(discriminant);
+        const parameter_t discriminant_sqrt = std::sqrt(discriminant);
 
-        const distance_t minus_p_h = - p / 2., r1 = minus_p_h + discriminant_sqrt, r2 = minus_p_h - discriminant_sqrt;
-        const distance_t lambda1 = std::min(r1, r2), lambda2 = std::max(r1, r2);
+        const parameter_t minus_p_h = - p / 2., r1 = minus_p_h + discriminant_sqrt, r2 = minus_p_h - discriminant_sqrt;
+        const parameter_t lambda1 = std::min(r1, r2), lambda2 = std::max(r1, r2);
 
-        return Interval(std::max(0., lambda1), std::min(1., lambda2));
+        return Interval(std::max(0.L, lambda1), std::min(1.L, lambda2));
     }
 
     inline auto as_ndarray() const {
         np::dtype dt = np::dtype::get_builtin<coordinate_t>();
         p::list l;
         np::ndarray result = np::array(l, dt);
-        for (const auto &elem : *this) {
+        for (const coordinate_t &elem : *this) {
             l.append(elem);
         }
         result = np::array(l, dt);
@@ -223,7 +223,7 @@ class Points : public std::vector<Point> {
 
     inline auto as_ndarray() const {
         p::list l;
-        for (const auto &elem : *this) {
+        for (const Point &elem : *this) {
             l.append(elem.as_ndarray());
         }
         auto result = np::array(l);

experiments/dtw_old.py

@@ -75,11 +75,11 @@ class Custom_Probability_Generator {
     }
 
     inline T get() {
-        const auto n = cumulative_probabilities.size();
-        const auto r = uform_gen.get();
+        const std::size_t n = cumulative_probabilities.size();
+        const T r = uform_gen.get();
         const auto upper = std::upper_bound(cumulative_probabilities.cbegin(), cumulative_probabilities.cend(), r);
         assert(upper != cumulative_probabilities.cend());
-        const auto result = std::distance(cumulative_probabilities.cbegin(), upper);
+        const std::size_t result = std::distance(cumulative_probabilities.cbegin(), upper);
         return result;
     }