OpenMP and Tidy Up

LICENSE

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2020 Dennis Rohde
+Copyright (c) 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

README.md

@@ -63,9 +63,9 @@ A `fred.Distance_Matrix()` can be used to speed up consecutive calls of `fred.di
 
 #### 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)
-- signature: `fred.discrete_klcenter_multi(k, l, curves, distances, center_domain, random_first_center)` with parameters 
+- signature: `fred.discrete_klcenter(k, l, curves, distances, random_first_center, fast_simplification)` with parameters 
     - `k`: number of centers
-    - `l`: maximum complexity of the centers, only used when center_domain is default value
+    - `l`: maximum complexity of the centers
     - `distances`: `fred.Distance_Matrix`, defaults to empty `fred.Distance_Matrix`
     - `random_first_center`: determines if first center is chosen uniformly at random or first curve is used as first center, optional, defaults to true
     - `fast_simplification`: determines whether to use the weak minimum error simplification or the faster approximate weak minimum error simplification, defaults to false
@@ -76,9 +76,9 @@ A `fred.Distance_Matrix()` can be used to speed up consecutive calls of `fred.di
 
 #### discrete (k,l)-median clustering (continuous Fréchet)
 - Algorithm from section 4.3 in [**Geometric Approximation Algorithms**](http://www.ams.org/books/surv/173/) + simplification
-- signature: `fred.discrete_klmedian_multi(k, l, curves, distances, center_domain)` with parameters 
+- signature: `fred.discrete_klmedian(k, l, curves, distances, fast_simplification)` with parameters 
     - `k`: number of centers
-    - `l`: maximum complexity of the centers, only used when center_domain is default value
+    - `l`: maximum complexity of the centers
     - `distances`: `fred.Distance_Matrix`, defaults to empty `fred.Distance_Matrix`
     - `fast_simplification`: determines whether to use the weak minimum error simplification or the faster approximate weak minimum error simplification, defaults to false
 - returns: `fred.Clustering_Result` with mebers 
@@ -105,6 +105,7 @@ A `fred.Distance_Matrix()` can be used to speed up consecutive calls of `fred.di
 ### Requirements
 
 You have to have installed:
+ - cmake
  - git
  - openmp available (should be a part of your compiler)
 

include/clustering.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 

include/coreset.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 

include/curve.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
@@ -186,7 +186,7 @@ class Curves : public std::vector<Curve> {
         return py::array_t<coordinate_t>(l);
     }
 
-    Curves simplify(const curve_size_t);
+    Curves simplify(const curve_size_t, const bool);
 
     std::string str() const;
 

include/frechet.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 

include/grid.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 

include/interval.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
@@ -20,21 +20,21 @@ class Interval {
     parameter_t beg, en;
 
 public:
-    Interval() : beg{1.L}, en{0.L} {}
+    Interval() : beg{parameter_t(1)}, en{parameter_t(0)} {}
 
     Interval(const parameter_t begin, const parameter_t end) : beg{begin}, en{end} {}
 
     inline bool operator<(const Interval &other) const {
         return (beg < other.begin()) or ((beg == other.begin()) and (en < other.end()));
     }
 
-    inline bool is_empty() const { 
+    inline bool empty() const { 
         if (en - beg >= std::numeric_limits<parameter_t>::epsilon()) return beg > en;
         else return true;
     }
 
     inline bool intersects(const Interval &other) const {
-        if (is_empty() or other.is_empty()) return false;
+        if (empty() or other.empty()) return false;
 
         return ((other.beg >= beg) and (other.beg <= en)) or
             ((other.en >= beg) and (other.en <= en)) or
@@ -50,8 +50,8 @@ class Interval {
     }
 
     inline void reset() {
-        beg = 1.L;
-        en = 0.L;
+        beg = parameter_t(1);
+        en = parameter_t(0);
     }
 };
 

include/jl_transform.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 

include/point.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
@@ -24,7 +24,7 @@ THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLI
 
 namespace py = pybind11;
 
-class Point : public Coordinates {    
+class Point : public Coordinates {
 public:    
     inline Point(const dimensions_t d) : Coordinates(d) {}
 
@@ -78,7 +78,7 @@ class Point : public Coordinates {
         Point result = *this;
         #pragma omp simd
         for (dimensions_t i = 0; i < dimensions(); ++i){
-            result.operator[](i) += point[i];
+            result[i] += point[i];
         }
         return result;
     }
@@ -87,7 +87,7 @@ class Point : public Coordinates {
         Point result = *this;
         #pragma omp simd
         for (dimensions_t i = 0; i < dimensions(); ++i){
-            result.operator[](i) -= point[i];
+            result[i] -= point[i];
         }
         return result;
     }
@@ -96,7 +96,7 @@ class Point : public Coordinates {
         Point result = *this;
         #pragma omp simd
         for (dimensions_t i = 0; i < dimensions(); ++i){
-            result.operator[](i) *= mult;
+            result[i] *= mult;
         }
         return result;
     }
@@ -114,7 +114,7 @@ class Point : public Coordinates {
         Point result = *this;
         #pragma omp simd
         for (dimensions_t i = 0; i < dimensions(); ++i){
-            result.operator[](i) /= dist;
+            result[i] /= dist;
         }
         return result;
     }
@@ -166,7 +166,7 @@ class Point : public Coordinates {
         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.L, lambda1), std::min(1.L, lambda2));
+        return Interval(std::max(parameter_t(0), lambda1), std::min(parameter_t(1), lambda2));
     }
 
      inline auto as_ndarray() const {

include/random.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 

include/simplification.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
@@ -36,7 +36,6 @@ class Subcurve_Shortcut_Graph {
         auto distance = Frechet::Continuous::Distance();
 
         for (curve_size_t i = 0; i < complexity - 1; ++i) {
-
             for (curve_size_t j = i + 1; j < complexity; ++j) {
 
                 curve.set_subcurve(i, j);
@@ -75,24 +74,17 @@ class Subcurve_Shortcut_Graph {
         for (curve_size_t i = 0; i < l; ++i) {
 
             if (i == 0) {
-
+                #pragma omp parallel for
                 for (curve_size_t j = 1; j < curve.complexity(); ++j) {
-
                     distances[j][0] = edges[0][j];
                     predecessors[j][0] = 0;
-
                 }
-
             } else {
-
                 for (curve_size_t j = 1; j < curve.complexity(); ++j) {
-
                     others.resize(j);
-
+                    #pragma omp parallel for
                     for (curve_size_t k = 0; k < j; ++k) {
-
                         others[k] = std::max(distances[k][i - 1], edges[k][j]);
-
                     }
 
                     best = std::distance(others.begin(), std::min_element(others.begin(), others.end()));

include/types.hpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 

setup.py

@@ -74,7 +74,7 @@ def build_extension(self, ext):
 
 setup(
     name='Fred-Frechet',
-    version='1.7.5',
+    version='1.7.7',
     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/curve.cpp

@@ -1,5 +1,5 @@
 /*
-Copyright 2020 Dennis Rohde
+Copyright 2020-2021 Dennis Rohde
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
 
@@ -23,34 +23,30 @@ Curve::Curve(const Points &points, const std::string &name) : Points(points), vs
     #endif
 }
 
-Curve::Curve(const py::array_t<coordinate_t> &in, const std::string &name) : Points(in.request().shape[0], in.request().ndim > 1 ? in.request().shape[1] : 1), name{name}, vstart{0}, vend{Points::size() - 1} {    
-    const dimensions_t n_dimensions = in.ndim();
-    auto shape = in.shape();
-    const curve_size_t number_points = shape[0];
+Curve::Curve(const py::array_t<coordinate_t> &in, const std::string &name) : Points(in.request().shape[0], in.request().ndim > 1 ? in.request().shape[1] : 1), name{name}, vstart{0}, vend{Points::size() - 1} {
+    const auto array_dim = in.ndim();
 
-    if (n_dimensions > 2){
+    if (array_dim > 2){
         std::cerr << "A Curve requires a 1- or 2-dimensional numpy array of type " << typeid(coordinate_t).name() << "." << std::endl;
-        std::cerr << "Current dimensions: " << n_dimensions << std::endl;
+        std::cerr << "Current dimensions: " << dimensions() << std::endl;
         std::cerr << "WARNING: constructed empty curve" << std::endl;
         return;
     }
 
-    if (n_dimensions == 2) {
-        const dimensions_t point_size = shape[1];
-
+    if (array_dim == 2) {        
         #if DEBUG
-        std::cout << "constructing curve of size " << number_points << " and " << point_size << " dimensions" << std::endl;
+        std::cout << "constructing curve of size " << number() << " and " << dimensions() << " dimensions" << std::endl;
         #endif
 
-        #pragma omp parallel for simd
-        for (curve_size_t i = 0; i < number_points; ++i) {
-            for(dimensions_t j = 0; j < point_size; ++j){
+        #pragma omp parallel for simd collapse(2)
+        for (curve_size_t i = 0; i < number(); ++i) {
+            for(dimensions_t j = 0; j < dimensions(); ++j){
               Points::operator[](i)[j] = *in.data(i, j);
             }
         }
     } else {
         #pragma omp parallel for simd
-        for (curve_size_t i = 0; i < number_points; ++i) {
+        for (curve_size_t i = 0; i < number(); ++i) {
             Points::operator[](i)[0] = *in.data(i);
         }
     }
@@ -61,8 +57,7 @@ Curve::Curve(const py::array_t<coordinate_t> &in, const std::string &name) : Poi
     }
 }
 
-Curves Curves::simplify(const curve_size_t l) {
-    const bool approx = true;
+Curves Curves::simplify(const curve_size_t l, const bool approx = false) {
     Curves result(size(), l, Curves::dimensions());
     for (curve_number_t i = 0; i < size(); ++i) {
         if (approx) {

src/dynamic_time_warping.cpp

@@ -29,9 +29,11 @@ Distance distance(const Curve &curve1, const Curve &curve2) {
     auto start = std::chrono::high_resolution_clock::now();
     std::vector<std::vector<distance_t>> a(curve1.complexity() + 1, std::vector<distance_t>(curve2.complexity() + 1, std::numeric_limits<distance_t>::infinity()));
     a[0][0] = 0;
+    #pragma omp parallel for ordered collapse(2)
     for (curve_size_t i = 1; i <= curve1.complexity(); ++i) {
         for (curve_size_t j = 1; j <= curve2.complexity(); ++j) {
             a[i][j] = std::sqrt(curve1[i-1].dist_sqr(curve2[j-1]));
+            #pragma omp ordered
             a[i][j] += std::min(std::min(a[i-1][j], a[i][j-1]), a[i-1][j-1]);
         }
     }