Minor fixes

README.md

@@ -74,6 +74,7 @@ The parameter `distance_func` controls which distance function to use. Possible
 #### Consecutive call option
 
 The parameter `consecutive_call` controls whether distances and simplifications already computed in a previous clustering call are resused (set to `true` then); defaults to `false`.
+Has no effect when `fred.config.use_distance_matrix == False`.
 
 #### discrete (k,l)-center clustering (continuous Fréchet)
 - from [**Approximating (k,l)-center clustering for curves**](https://dl.acm.org/doi/10.5555/3310435.3310616)

include/dynamic_time_warping.hpp

@@ -32,6 +32,7 @@ namespace Discrete {
         std::string repr() const;
 
         std::vector<std::pair<curve_number_t, curve_number_t>> matching;
+        curve_size_t n, m;
     };
 
     Points vertices_matching_points(const Curve&, const Curve&, const Distance&);

setup.py

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

src/clustering.cpp

@@ -86,36 +86,58 @@ py::list Clustering_Result::compute_center_enclosing_balls(const Curves &in, con
         center_matching_points.push_back(std::vector<Points>(get(i).complexity(), get(i).dimensions()));
     }
 
+    curve_number_t ii, jj;    
+    Points tpoints(0);
+
     for (curve_number_t i = 0; i < size(); ++i) {
         if (Config::verbosity > 2) py::print("Clustering Result: computing points for center ", i);
+
         py::list center_list;
-        for (curve_number_t j = 0; j < (*assignment)[i].size(); ++j) {
-            Points tpoints(get(i).dimensions());
+        const Curve &center_curve = get(i);
+
+        for (curve_number_t j = 0; j < (*assignment)[i].size(); ++j) {           
+            ii = (*assignment)[i][j];
+            const Curve &input_curve = in[ii];
+
+            if (consecutive_call) jj = center_indices[i];
+            else jj = i;
+
             switch (distance_func) {
                 case 0:
-                    if (Config::use_distance_matrix) tpoints = Frechet::Continuous::vertices_matching_points(get(i), in[(*assignment)[i][j]], dynamic_cast<const Frechet::Continuous::Distance&>(*distances[(*assignment)[i][j]][i]));
-                    else {
-                        Frechet::Continuous::Distance curr_dist = Frechet::Continuous::distance(get(i), in[(*assignment)[i][j]]);
-                        tpoints = Frechet::Continuous::vertices_matching_points(get(i), in[(*assignment)[i][j]], curr_dist);
+                    if (Config::use_distance_matrix) {
+                        const auto &dist = dynamic_cast<const Frechet::Continuous::Distance&>(*distances[ii][jj]);
+                        tpoints = Frechet::Continuous::vertices_matching_points(input_curve, center_curve, dist);
+                    } else {
+                        const Frechet::Continuous::Distance dist = Frechet::Continuous::distance(input_curve, center_curve);
+                        tpoints = Frechet::Continuous::vertices_matching_points(input_curve, center_curve, dist);
                     }
                     break;
                 case 1:
                     break;
                 case 2:
-                    if (Config::use_distance_matrix) tpoints = Dynamic_Time_Warping::Discrete::vertices_matching_points(get(i), in[(*assignment)[i][j]], dynamic_cast<const Dynamic_Time_Warping::Discrete::Distance&>(*distances[(*assignment)[i][j]][i]));
-                    else {
-                        Dynamic_Time_Warping::Discrete::Distance curr_dist = Dynamic_Time_Warping::Discrete::distance(get(i), in[(*assignment)[i][j]]);
-                        tpoints = Dynamic_Time_Warping::Discrete::vertices_matching_points(get(i), in[(*assignment)[i][j]], curr_dist);
+                    if (Config::use_distance_matrix) {
+                        const auto &dist = dynamic_cast<const Dynamic_Time_Warping::Discrete::Distance&>(*distances[ii][jj]);
+                        tpoints = Dynamic_Time_Warping::Discrete::vertices_matching_points(input_curve, center_curve, dist);
+                    } else {
+                        const Dynamic_Time_Warping::Discrete::Distance dist = Dynamic_Time_Warping::Discrete::distance(input_curve, center_curve);
+                        tpoints = Dynamic_Time_Warping::Discrete::vertices_matching_points(input_curve, center_curve, dist);
                     }
                     break;
                 default:
                     py::print("not implemented!");
             }
-            for (curve_size_t k = 0; k < get(i).complexity(); ++k) { 
-                center_matching_points[i][k].push_back(tpoints[k]);
+
+            if (Config::verbosity > 2) py::print("Clustering Result: collecting matching points for center ", i);
+
+            for (curve_size_t k = 0; k < center_curve.complexity(); ++k) { 
+                if (Config::verbosity > 3) py::print("Clustering Result: collecting matching points for center ", i, " vertex ", k);
+                center_matching_points[i][k].push_back(tpoints.at(k));
             }
         }
-        for (curve_size_t k = 0; k < get(i).complexity(); ++k) {
+
+        if (Config::verbosity > 2) py::print("Clustering Result: Computing input for stabbing algorithm");
+
+        for (curve_size_t k = 0; k < center_curve.complexity(); ++k) {
             py::list b_r;
             switch (distance_func) {
                 case 0:
@@ -130,6 +152,7 @@ py::list Clustering_Result::compute_center_enclosing_balls(const Curves &in, con
                     break;
                 case 2:
                     {
+                        if (Config::verbosity > 2) py::print("Clustering Result: Computing mean of matching points for vertex ", k);
                         const auto mean = center_matching_points[i][k].centroid();
                         b_r.append(mean.as_ndarray());
                         b_r.append(0);
@@ -329,6 +352,7 @@ Clustering_Result kl_cluster(const curve_number_t num_centers, const curve_size_
                 }
             }
         }
+        curr_maxdist = cost;
     }
 
     if (median) {

src/dynamic_time_warping.cpp

@@ -24,28 +24,38 @@ std::string Distance::repr() const {
     return ss.str();
 }
 
-Points vertices_matching_points(const Curve &curve1, const Curve &curve2, const Distance &dist) {
-    if ((curve1.complexity() < 2) or (curve2.complexity() < 2)) {
+Points vertices_matching_points(const Curve &input_curve, const Curve &center_curve, const Distance &dist) {
+    if ((input_curve.complexity() < 2) or (center_curve.complexity() < 2)) {
         py::print("WARNING: curves must be of at least two points");
-        Points result(curve1.dimensions());
+        Points result(center_curve.dimensions());
         return result;
     }
+
+    if (Config::verbosity > 1) py::print("DDTW: computing matching points from center_curve of complexity ", center_curve.complexity(), " to input_curve of complexity ", input_curve.complexity());
+    if (Config::verbosity > 2) py::print("DDTW: distance between input_curve and center_curve is ", dist.value);
 
-    if (Config::verbosity > 1) py::print("DDTW: computing matching points from curve1 of complexity ", curve1.complexity(), " to curve2 of complexity ", curve2.complexity());
-    if (Config::verbosity > 2) py::print("DDTW: distance between curve1 and curve2 is ", dist.value);
+    std::vector<Points> matching_points(center_curve.size(), Points(center_curve.dimensions()));    
 
-    std::vector<Points> matching_points(curve1.size(), Points(curve1.dimensions()));    
+    curve_number_t j, k;
 
     for (curve_number_t i = 0; i < dist.matching.size(); ++i) {
-        curve_number_t j = dist.matching[i].first, k = dist.matching[i].second;
-        matching_points[j].push_back(curve2[k]);
+        j = dist.matching[i].first;
+        k = dist.matching[i].second;
+
+        if (Config::verbosity > 2) py::print("DDTW: matching point ", j, " on input_curve to point ", k, " on curve 2");
+        matching_points[k].push_back(input_curve[j]);
     }
 
-    Points result(curve1.size(), curve1.dimensions());
+    Points result(center_curve.size(), center_curve.dimensions());
 
-    for (curve_size_t i = 0; i < curve1.size(); ++i) {
+    if (Config::verbosity > 2) py::print("DDTW: computing centroids to aggregate multi-matching points");
+
+    for (curve_size_t i = 0; i < center_curve.size(); ++i) {
+        if (Config::verbosity > 2) py::print("DDTW: computing centroid ", i);
         result[i] = matching_points[i].centroid();
     }
+
+    if (Config::verbosity > 2) py::print("DDTW: matching points computed");
 
     return result;
 }
@@ -70,12 +80,12 @@ Distance distance(const Curve &curve1, const Curve &curve2) {
 
     std::vector<std::vector<std::pair<curve_number_t, curve_number_t>>> multi_warp_counter(n1 + 1, std::vector<std::pair<curve_number_t, curve_number_t>>(n2 + 1, std::make_pair(0, 0)));
     std::vector<std::vector<std::pair<curve_number_t, curve_number_t>>> b(n1 + 1, std::vector<std::pair<curve_number_t, curve_number_t>>(n2 + 1, std::make_pair(0, 0)));
-    std::vector<std::vector<distance_t>> a(curve1.complexity() + 1, std::vector<distance_t>(curve2.complexity() + 1, infty));
-    std::vector<std::vector<distance_t>> dists(curve1.complexity(), std::vector<distance_t>(curve2.complexity()));
+    std::vector<std::vector<distance_t>> a(n1 + 1, std::vector<distance_t>(n2 + 1, infty));
+    std::vector<std::vector<distance_t>> dists(n1, std::vector<distance_t>(n2));
 
     #pragma omp parallel for collapse(2)
-    for (curve_size_t i = 0; i < curve1.complexity(); ++i) {
-        for (curve_size_t j = 0; j < curve2.complexity(); ++j) {
+    for (curve_size_t i = 0; i < n1; ++i) {
+        for (curve_size_t j = 0; j < n2; ++j) {
             dists[i][j] = curve1[i].dist(curve2[j]);
         }
     }
@@ -86,8 +96,8 @@ Distance distance(const Curve &curve1, const Curve &curve2) {
     distance_t min_ele;
 
     a[0][0] = 0;
-    for (curve_size_t i = 1; i <= curve1.complexity(); ++i) {
-        for (curve_size_t j = 1; j <= curve2.complexity(); ++j) {
+    for (curve_size_t i = 1; i <= n1; ++i) {
+        for (curve_size_t j = 1; j <= n2; ++j) {
             a[i][j] = dists[i-1][j-1];
             mwd = 0;
 
@@ -127,6 +137,8 @@ Distance distance(const Curve &curve1, const Curve &curve2) {
     }
 
     result.matching.push_back(std::make_pair(0, 0));
+    result.n = n1;
+    result.m = n2;
 
     const auto end = std::clock();
     result.time = (end - start) / CLOCKS_PER_SEC;

src/frechet.cpp

@@ -27,31 +27,31 @@ std::string Distance::repr() const {
     return ss.str();
 }
 
-Points vertices_matching_points(const Curve &curve1, const Curve &curve2, const Distance &dist) {
-    if ((curve1.complexity() < 2) or (curve2.complexity() < 2)) {
+Points vertices_matching_points(const Curve &input_curve, const Curve &center_curve, const Distance &dist) {
+    if ((center_curve.complexity() < 2) or (input_curve.complexity() < 2)) {
         py::print("WARNING: curves must be of at least two points");
-        Points result(curve1.dimensions());
+        Points result(center_curve.dimensions());
         return result;
     }
 
-    if (Config::verbosity > 1) py::print("CFD: computing matching points for distance ", dist.value, " from curve1 of complexity ", curve1.complexity(), " to curve2 of complexity ", curve2.complexity());
+    if (Config::verbosity > 1) py::print("CFD: computing matching points for distance ", dist.value, " from center_curve of complexity ", center_curve.complexity(), " to input_curve of complexity ", input_curve.complexity());
 
     const distance_t dist_sqr = dist.value * dist.value;
-    const curve_size_t n1 = curve1.complexity();
-    const curve_size_t n2 = curve2.complexity();
+    const curve_size_t n1 = center_curve.complexity();
+    const curve_size_t n2 = input_curve.complexity();
 
     std::vector<Intervals> free_intervals(n1, Intervals(n2, Interval()));
 
     #pragma omp parallel for collapse(2) if (n1 * n2 > 1000)
     for (curve_size_t i = 1; i < n1; ++i) {
         for (curve_size_t j = 0; j < n2 - 1; ++j) {
-            free_intervals[i][j] = curve1[i].ball_intersection_interval(dist_sqr, curve2[j], curve2[j+1]);
+            free_intervals[i][j] = center_curve[i].ball_intersection_interval(dist_sqr, input_curve[j], input_curve[j+1]);
         }
     }
 
     if (Config::verbosity > 1) py::print("CFD: free space computed, computing matching");
 
-    Points result(n1, curve1.dimensions());
+    Points result(n1, center_curve.dimensions());
     parameter_t p;
     curve_size_t jj(0);
 
@@ -68,11 +68,11 @@ Points vertices_matching_points(const Curve &curve1, const Curve &curve2, const
                 break;
             }
         }
-        result[i] = curve2[jj].line_segment_point(curve2[jj+1], p);
-        if (Config::verbosity > 1) py::print("CFD: matching vertex ", i, "to ", p, "on segment ", jj, " to ", jj+1, " with distance ", result[i].dist(curve1[i]));
+        result[i] = input_curve[jj].line_segment_point(input_curve[jj+1], p);
+        if (Config::verbosity > 1) py::print("CFD: matching vertex ", i, "to ", p, "on segment ", jj, " to ", jj+1, " with distance ", result[i].dist(center_curve[i]));
     }
-    result[0] = curve2[0];
-    result[n1-1] = curve2[n2-1];
+    result[0] = input_curve[0];
+    result[n1-1] = input_curve[n2-1];
 
     return result;
 }