improve random projection parallelization

setup.py

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

@@ -28,7 +28,7 @@ Curve::Curve(const py::array_t<coordinate_t> &in, const std::string &name) : Poi
 
     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: " << dimensions() << std::endl;
+        std::cerr << "Current dimensions: " << array_dim << std::endl;
         std::cerr << "WARNING: constructed empty curve" << std::endl;
         return;
     }

src/jl_transform.cpp

@@ -25,27 +25,31 @@ Curves transform_naive(const Curves &in, const distance_t epsilon, const bool em
     const distance_t epsiloncu = epsilonsq * epsilon;
     const dimensions_t new_number_dimensions = empirical_k ? std::ceil(2 * std::log(number_points) * 1/epsilonsq):
         std::ceil(4 * std::log(number_points) * 1 /((epsilonsq/2) - (epsiloncu/3)));
-
-    std::vector<Coordinates> mat (new_number_dimensions);
+    const dimensions_t old_number_dimensions = in.dimensions();
+
+    std::vector<Coordinates> mat(new_number_dimensions);
 
     #if DEBUG
-    std::cout << "populating " << new_number_dimensions << "x" << in[0].dimensions() << " matrix" << std::endl;
+    std::cout << "populating " << new_number_dimensions << "x" << old_number_dimensions << " matrix" << std::endl;
     #endif
 
-    for (Coordinates &elem : mat) elem = rg.get(in[0].dimensions());
+    #pragma omp parallel for
+    for (dimensions_t i = 0; i < new_number_dimensions; ++i) {
+        mat[i] = rg.get(old_number_dimensions);
+    }
 
     Curves result(in.size(), in.get_m(), new_number_dimensions);
 
     coordinate_t sqrtk = std::sqrt(new_number_dimensions);
 
+    #pragma omp parallel for
+    for (curve_number_t l = 0; l < in.size(); ++l) result[l] = Curve(in[l].complexity(), new_number_dimensions, in[l].get_name());
 
     for (curve_number_t l = 0; l < in.size(); ++l) {
-        result[l] = Curve(in[l].complexity(), new_number_dimensions, in[l].get_name());
-
         #pragma omp parallel for collapse(2)
         for (curve_size_t i = 0; i < in[l].complexity(); ++i) {
-
             for (dimensions_t j = 0; j < new_number_dimensions; ++j) {
+
                 result[l][i][j] = mat[j][0] * in[l][i][0];
 
                 for (dimensions_t k = 1; k < in[l].dimensions(); ++k) {
@@ -55,16 +59,13 @@ Curves transform_naive(const Curves &in, const distance_t epsilon, const bool em
                 result[l][i][j] /= sqrtk;
 
             }
-
         }
 
         #if DEBUG
         std::cout << "projected curve no. " << l << " from " << in.dimensions() << " to " << new_number_dimensions << " dimensions" << std::endl;
         #endif
-
     }
 
     return result;
 }
-
 }