draft: Cythonize

similaritymeasures/cy_similaritymeasures/benchmark.py

@@ -0,0 +1,82 @@
+'''
+just a quick benchmaerking script to showcase cythonized vs
+non cythonized performances
+'''
+
+import timeit
+import numpy as np
+from similaritymeasures import frechet_dist, dtw
+from cy_similaritymeasures import frechet_dist as cy_frechet_dist
+from cy_similaritymeasures import dtw as cy_dtw
+
+#generating data in a way similar to the tests
+x1 = np.linspace(0.0, 1.0, 100)
+y1 = np.ones(100)*2
+x2 = np.linspace(0.0, 1.0, 50)
+y2 = np.ones(50)
+
+np.random.seed(1212121)
+curve_a_rand = np.random.random((100, 2))
+curve_b_rand = np.random.random((90, 2))
+
+curve1 = np.array((x1, y1)).T
+curve2 = np.array((x2, y2)).T
+
+res = frechet_dist(curve1, curve2)
+
+cy_res = cy_frechet_dist(curve1, curve2)
+
+print(res)
+print(cy_res)
+print(res == cy_res)
+
+def run_frechet():
+ return frechet_dist(curve1, curve2)
+
+def run_cy_frechet():
+ return cy_frechet_dist(curve1, curve2)
+
+n_repeats = 5
+n_runs = 200
+
+t1 = timeit.repeat(run_frechet, repeat = n_repeats, number = n_runs)
+
+t2 = timeit.repeat(run_cy_frechet, repeat = n_repeats, number = n_runs)
+print(t1)
+print(t2)
+
+avg = sum(t1) / len(t1)
+cy_avg = sum(t2) / len(t2)
+improv_pct = (avg - cy_avg) / avg * 100
+print('average execution time for non cythonized version: %f' % avg)
+print('average execution time for cythonized version: %f' % cy_avg)
+print('improvement: %d %%' % improv_pct)
+
+#then i would benchmark dtw
+
+res, _ = dtw(curve1, curve2)
+
+cy_res, _ = cy_dtw(curve1, curve2)
+
+print(res)
+print(cy_res)
+print(res == cy_res)
+
+def run_dtw():
+ return dtw(curve1, curve2)
+
+def run_cy_dtw():
+ return cy_dtw(curve1, curve2)
+
+t1 = timeit.repeat(run_dtw, repeat = n_repeats, number = n_runs)
+
+t2 = timeit.repeat(run_cy_dtw, repeat = n_repeats, number = n_runs)
+print(t1)
+print(t2)
+
+avg = sum(t1) / len(t1)
+cy_avg = sum(t2) / len(t2)
+improv_pct = (avg - cy_avg) / avg * 100
+print('average execution time for non cythonized version: %f' % avg)
+print('average execution time for cythonized version: %f' % cy_avg)
+print('improvement: %d %%' % improv_pct)

similaritymeasures/cy_similaritymeasures/build_cython_extension.sh

@@ -0,0 +1 @@
+python3 setup.py build_ext --inplace

similaritymeasures/cy_similaritymeasures/cy_similaritymeasures.pyx

@@ -0,0 +1,212 @@
+import numpy as np
+from scipy.spatial import distance
+
+cdef double c_max(double a, double b):
+ if a > b:
+ return a
+ return b
+
+cdef double c_min(double a, double b):
+ if a < b:
+ return a
+ return b
+
+def frechet_dist(exp_data, num_data, p=2):
+ r"""
+ Compute the discrete Frechet distance
+
+ Compute the Discrete Frechet Distance between two N-D curves according to
+ [1]_. The Frechet distance has been defined as the walking dog problem.
+ From Wikipedia: "In mathematics, the Frechet distance is a measure of
+ similarity between curves that takes into account the location and
+ ordering of the points along the curves. It is named after Maurice Frechet.
+ https://en.wikipedia.org/wiki/Fr%C3%A9chet_distance
+
+ Parameters
+ ----------
+ exp_data : array_like
+ Curve from your experimental data. exp_data is of (M, N) shape, where
+ M is the number of data points, and N is the number of dimmensions
+ num_data : array_like
+ Curve from your numerical data. num_data is of (P, N) shape, where P
+ is the number of data points, and N is the number of dimmensions
+ p : float, 1 <= p <= infinity
+ Which Minkowski p-norm to use. Default is p=2 (Eculidean).
+ The manhattan distance is p=1.
+
+ Returns
+ -------
+ df : float
+ discrete Frechet distance
+
+ References
+ ----------
+ .. [1] Thomas Eiter and Heikki Mannila. Computing discrete Frechet
+ distance. Technical report, 1994.
+ http://www.kr.tuwien.ac.at/staff/eiter/et-archive/cdtr9464.pdf
+ http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.90.937&rep=rep1&type=pdf
+
+ Notes
+ -----
+ Your x locations of data points should be exp_data[:, 0], and the y
+ locations of the data points should be exp_data[:, 1]. Same for num_data.
+
+ Thanks to Arbel Amir for the issue, and Sen ZHANG for the iterative code
+ https://github.com/cjekel/similarity_measures/issues/6
+
+ Examples
+ --------
+ >>> # Generate random experimental data
+ >>> x = np.random.random(100)
+ >>> y = np.random.random(100)
+ >>> exp_data = np.zeros((100, 2))
+ >>> exp_data[:, 0] = x
+ >>> exp_data[:, 1] = y
+ >>> # Generate random numerical data
+ >>> x = np.random.random(100)
+ >>> y = np.random.random(100)
+ >>> num_data = np.zeros((100, 2))
+ >>> num_data[:, 0] = x
+ >>> num_data[:, 1] = y
+ >>> df = frechet_dist(exp_data, num_data)
+
+ """
+ cdef int i, j
+ cdef int n = len(exp_data)
+ cdef int m = len(num_data)
+ c = distance.cdist(exp_data, num_data, metric='minkowski', p=p)
+ #ca = np.ones((n, m))
+ #ca = np.multiply(ca, -1)
+ #ca[0, 0] = c[0, 0]
+ #for i in range(1, n):
+ # ca[i, 0] = max(ca[i-1, 0], c[i, 0])
+ #for j in range(1, m):
+ # ca[0, j] = max(ca[0, j-1], c[0, j])
+ for i in range(1, n):
+ for j in range(1, m):
+ c[i, j] = c_max(
+ c_min(c_min(c[i-1, j], c[i, j-1]) , c[i-1, j-1]),
+ c[i, j])
+ return c[n-1, m-1]
+
+def dtw(exp_data, num_data, metric='euclidean', **kwargs):
+ r"""
+ Compute the Dynamic Time Warping distance.
+
+ This computes a generic Dynamic Time Warping (DTW) distance and follows
+ the algorithm from [1]_. This can use all distance metrics that are
+ available in scipy.spatial.distance.cdist.
+
+ Parameters
+ ----------
+ exp_data : array_like
+ Curve from your experimental data. exp_data is of (M, N) shape, where
+ M is the number of data points, and N is the number of dimmensions
+ num_data : array_like
+ Curve from your numerical data. num_data is of (P, N) shape, where P
+ is the number of data points, and N is the number of dimmensions
+ metric : str or callable, optional
+ The distance metric to use. Default='euclidean'. Refer to the
+ documentation for scipy.spatial.distance.cdist. Some examples:
+ 'braycurtis', 'canberra', 'chebyshev', 'cityblock', 'correlation',
+ 'cosine', 'dice', 'euclidean', 'hamming', 'jaccard', 'kulsinski',
+ 'mahalanobis', 'matching', 'minkowski', 'rogerstanimoto', 'russellrao',
+ 'seuclidean', 'sokalmichener', 'sokalsneath', 'sqeuclidean',
+ 'wminkowski', 'yule'.
+ **kwargs : dict, optional
+ Extra arguments to `metric`: refer to each metric documentation in
+ scipy.spatial.distance.
+ Some examples:
+
+ p : scalar
+ The p-norm to apply for Minkowski, weighted and unweighted.
+ Default: 2.
+
+ w : ndarray
+ The weight vector for metrics that support weights (e.g.,
+ Minkowski).
+
+ V : ndarray
+ The variance vector for standardized Euclidean.
+ Default: var(vstack([XA, XB]), axis=0, ddof=1)
+
+ VI : ndarray
+ The inverse of the covariance matrix for Mahalanobis.
+ Default: inv(cov(vstack([XA, XB].T))).T
+
+ out : ndarray
+ The output array
+ If not None, the distance matrix Y is stored in this array.
+
+ Returns
+ -------
+ r : float
+ DTW distance.
+ d : ndarray (2-D)
+ Cumulative distance matrix
+
+ Notes
+ -----
+ The DTW distance is d[-1, -1].
+
+ This has O(M, P) computational cost.
+
+ The latest scipy.spatial.distance.cdist information can be found at
+ https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html
+
+ Your x locations of data points should be exp_data[:, 0], and the y
+ locations of the data points should be exp_data[:, 1]. Same for num_data.
+
+ This uses the euclidean distance for now. In the future it should be
+ possible to support other metrics.
+
+ DTW is a non-metric distance, which means DTW doesn't hold the triangle
+ inequality.
+ https://en.wikipedia.org/wiki/Triangle_inequality
+
+ References
+ ----------
+ .. [1] Senin, P., 2008. Dynamic time warping algorithm review. Information
+ and Computer Science Department University of Hawaii at Manoa Honolulu,
+ USA, 855, pp.1-23.
+ http://seninp.github.io/assets/pubs/senin_dtw_litreview_2008.pdf
+
+ Examples
+ --------
+ >>> # Generate random experimental data
+ >>> x = np.random.random(100)
+ >>> y = np.random.random(100)
+ >>> exp_data = np.zeros((100, 2))
+ >>> exp_data[:, 0] = x
+ >>> exp_data[:, 1] = y
+ >>> # Generate random numerical data
+ >>> x = np.random.random(100)
+ >>> y = np.random.random(100)
+ >>> num_data = np.zeros((100, 2))
+ >>> num_data[:, 0] = x
+ >>> num_data[:, 1] = y
+ >>> r, d = dtw(exp_data, num_data)
+
+ The euclidean distance is used by default. You can use metric and **kwargs
+ to specify different types of distance metrics. The following example uses
+ the city block or Manhattan distance between points.
+
+ >>> r, d = dtw(exp_data, num_data, metric='cityblock')
+
+ """
+ cdef int i, j
+ cdef int n = len(exp_data)
+ cdef int m = len(num_data)
+ c = distance.cdist(exp_data, num_data, metric=metric, **kwargs)
+
+ d = np.zeros(c.shape)
+ d[0, 0] = c[0, 0]
+ #n, m = c.shape
+ for i in range(1, n):
+ d[i, 0] = d[i-1, 0] + c[i, 0]
+ for j in range(1, m):
+ d[0, j] = d[0, j-1] + c[0, j]
+ for i in range(1, n):
+ for j in range(1, m):
+ d[i, j] = c[i, j] + c_min(c_min(d[i-1, j], d[i, j-1]), d[i-1, j-1])
+ return d[-1, -1], d

similaritymeasures/cy_similaritymeasures/readme.md

@@ -0,0 +1,26 @@
+# cy_similaritymeasures
+
+this folder shall contain stuff to cythonize the similaritymeasures functions
+
+```
+benchmark.py #shall be the script containg a quick performance check, to be executed at last after the cythonized version was correctly built
+build_cython_extension.sh #to invoke the build from the setup.py file a command is required, and for linux such command is written in this file
+cy_similaritymeasures.pyx #the cython code with the cythonized version of the frechet function
+setup.py #a setup.py file used to build the cython code and obtain an extension from it
+```
+
+## extension building
+
+cython code from the `*.pyx` file is not directly invokable from python, instead it should be built and transformed into an extension. To do that there is the script `build_cython_extension.sh`, which invokes `setup.py` to obtain a `.cpp` and `.so` extension file to be invoked from python.
+
+You should have cython on your machine in order to build everything.
+
+## results
+
+to run the benchmark just run the `benchmark.py`. On my machine the results were:
+
+```
+average execution time for non cythonized version: 0.744819
+average execution time for cythonized version: 0.433768
+improvement: 41 %
+```

similaritymeasures/cy_similaritymeasures/setup.py

@@ -0,0 +1,17 @@
+'''
+knowing this may be a bad thing, but i will place here a
+quick setup.py to just
+'''
+
+from setuptools import setup, Extension
+from Cython.Build import cythonize
+
+setup(
+ name = "cy_frechet",
+ ext_modules = cythonize(Extension(
+ "cy_similaritymeasures",
+ sources=["cy_similaritymeasures.pyx"],
+ language="c++"
+ )
+))
+