Merge branch 'master' of https://gricad-gitlab.univ-grenoble-alpes.fr/python-uga/training-hpc

pyfiles/dtw_cort_dist/V3_c_dtw_cort_vect/dtw_cort_dist_mat.py

 #!/usr/bin/env python
 
-import os
-import argparse
 import numpy as np
-from cffi import FFI
 import sysconfig
+from functools import partial
+from cffi import FFI
+from runpy import run_path
+from pathlib import Path
+
+util = run_path(Path(__file__).absolute().parent.parent / "util.py")
 
 ffi = FFI()
 ffi.cdef("double dtw(double* x, double* y, int size_x, int size_y);")
-my_dir = os.path.dirname(os.path.realpath(__file__))
+my_dir = Path(__file__).absolute().parent
 dllib = ffi.dlopen(
-    os.path.join(my_dir, "distances" + sysconfig.get_config_var("EXT_SUFFIX"))
+    str(my_dir / ("cdtw" + sysconfig.get_config_var("EXT_SUFFIX")))
+)
 
 
 def serie_pair_index_generator(number):
@@ -56,35 +60,7 @@ def cort(s1, s2):
     return num / np.sqrt(norm1 * norm2)
 
 
-if __name__ == "__main__":
-    PARSER = argparse.ArgumentParser(
-        description="Computes the distance matrix btw series"
-    )
-    PARSER.add_argument(
-        "input", type=str, help="input file ", nargs="?", default="../data.npy"
-    )
-    PARSER.add_argument(
-        "-n", type=int, default=None, help="number of series (default all 200)"
-    )
-    PARSER.add_argument(
-        "-v", action="store_true", default=None, help="visualize the matrix"
-    )
-
-    ARGS = PARSER.parse_args()
-
-    series = None
-    # load the input data
-    series = np.load(ARGS.input)
-
-    if ARGS.n:
-        nb_series = ARGS.n
-    else:
-        nb_series = series.shape[0]
-
-    from time import time
-
-    t0 = time()
-
+def compute(series, nb_series):
     gen = serie_pair_index_generator(nb_series)
 
     _dist_mat_dtw = np.zeros((nb_series, nb_series), dtype=np.float64)
@@ -97,14 +73,10 @@ if __name__ == "__main__":
         _dist_mat_cort[t1, t2] = dist_cort
         _dist_mat_cort[t2, t1] = dist_cort
 
-    print("\nelapsed time = {:.3f} s".format(time() - t0))
+    return _dist_mat_dtw, _dist_mat_cort
 
-    if ARGS.v:
-        import matplotlib.pyplot as plt
 
-        fig, axes = plt.subplots(2)
-        cax0 = axes[0].imshow(_dist_mat_dtw, cmap=plt.cm.gray)
-        axes[0].set_title("dtw")
-        cax1 = axes[1].imshow(_dist_mat_cort, cmap=plt.cm.gray)
-        axes[1].set_title("cort")
-        plt.show()
+main = partial(util["main"], compute)
+
+if __name__ == "__main__":
+    main()

pyfiles/dtw_cort_dist/V4_c_dtw_c_cort/Makefile

+
+all:
+	python3 setup.py develop
+
+clean:
+	rm -f cdtw.*.so
+	rm -rf build dtw_cort_dist_mat.egg-info

pyfiles/dtw_cort_dist/V4_c_dtw_c_cort/cdtw.c

+
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <float.h>
+#include <math.h>
+#include <float.h> /* for max_float */
+#include "cdtw.h"
+
+
+double min3(const double x, const double y, const double z){  
+  if (x < y)
+    if (x < z) return x;
+  if (y < z) return y;
+  return z;
+}
+
+/*
+ Returns the  minimum-distance warp path 
+ between multivariate time series x and y. 
+*/
+double dtw(double* x, double* y, int size_x, int size_y){
+  int i, j;
+  double* temp;
+  double *costP =  malloc(size_x*sizeof(double));
+  if (! costP) return -1;
+  double *costC = malloc(size_x*sizeof(double));
+  if (! costC){
+    free(costP);
+    return -1; /* TODO:  better handling of the error */
+  }
+  
+  costP[0] =  fabs(x[0] - y[0]);
+  for (i=1; i < size_x; i++){
+    costP[i] = fabs(x[i] - y[0]) + costP[i-1];
+  }
+  for (j=1; j< size_y; j++){
+    costC[0] = costP[0] + fabs(x[0] - y[j]);
+    for (i=1; i < size_x; i++){
+      costC[i] = min3(costC[i-1], costP[i-1], costP[i]) + fabs(x[i] - y[j]);
+    }
+    temp = costP;
+    costP = costC;
+    costC = temp;
+    
+  }
+  /* result is read in costP because switching between costP and costC is done */
+  double ret = costP[size_x-1];
+  free(costP);
+  free(costC);
+  return ret;
+}
+
+double cort(double *x, double *y, int size){
+  int t;
+  double slopex, slopey;
+  double num;
+  double sumSquareXslope, sumSquareYslope;
+
+  num = sumSquareXslope = sumSquareYslope = 0.0;
+  for (t = 0; t < size-1; t++){    
+    slopex = x[t+1] -  x[t]; 
+    slopey = y[t+1] -  y[t];
+    num += slopex * slopey;
+    sumSquareXslope += slopex * slopex;
+    sumSquareYslope += slopey * slopey;      
+  }
+  return num/sqrt(sumSquareXslope*sumSquareYslope);  
+}

pyfiles/dtw_cort_dist/V4_c_dtw_c_cort/cdtw.h

+double dtw(double* x, double* y, int size_x, int size_y);
+double cort(double *x, double *y, int size);

pyfiles/dtw_cort_dist/broken_V4_c_dtw_c_cort/dtw_cort_dist_mat.py → pyfiles/dtw_cort_dist/V4_c_dtw_c_cort/dtw_cort_dist_mat.py

 #!/usr/bin/python3
 
-import os
-import argparse
 import numpy as np
+import sysconfig
+from functools import partial
 from cffi import FFI
+from runpy import run_path
+from pathlib import Path
+
+util = run_path(Path(__file__).absolute().parent.parent / "util.py")
 
 ffi = FFI()
-ffi.cdef("""double dtw(double* x, double* y, int size_x, int size_y);
-            double cort(double *x, double *y, int size);""")
-my_dir = os.path.dirname(os.path.realpath(__file__))
-dllib = ffi.dlopen(os.path.join(my_dir, "build", "lib.linux-x86_64-2.7", "distances.so"))
+ffi.cdef("double dtw(double* x, double* y, int size_x, int size_y);")
+ffi.cdef("double cort(double* x, double* y, int size_x);")
+my_dir = Path(__file__).absolute().parent
+dllib = ffi.dlopen(
+    str(my_dir / ("cdtw" + sysconfig.get_config_var("EXT_SUFFIX")))
+)
 
 
 def serie_pair_index_generator(number):
@@ -44,23 +50,7 @@ def cort(serie_a, serie_b):
     ret = dllib.cort(a_ptr, b_ptr, len(serie_a))
     return ret
 
-
-if __name__ == "__main__":
-    PARSER = argparse.ArgumentParser(description="Computes the distance matrix btw series")
-    PARSER.add_argument("input", type=str, help="input file ")
-    PARSER.add_argument("-n", type=int, default=None, help="number of series (default all 200)")
-    PARSER.add_argument("-v", action="store_true", default=None, help="visualize the matrix")
-
-    ARGS = PARSER.parse_args()
-
-    series = None
-    # load the input data
-    series = np.load(ARGS.input)
-
-    if ARGS.n:
-        nb_series = ARGS.n
-    else:
-        nb_series = series.shape[0]
+def compute(series, nb_series):
     gen = serie_pair_index_generator(nb_series)
 
     _dist_mat_dtw = np.zeros((nb_series, nb_series), dtype=np.float64)
@@ -69,15 +59,14 @@ if __name__ == "__main__":
         dist_dtw = cDTW(series[t1], series[t2])
         _dist_mat_dtw[t1, t2] = dist_dtw
         _dist_mat_dtw[t2, t1] = dist_dtw
-        dist_cort = 0.5*(1-cort(series[t1], series[t2]))
+        dist_cort = 0.5 * (1 - cort(series[t1], series[t2]))
         _dist_mat_cort[t1, t2] = dist_cort
         _dist_mat_cort[t2, t1] = dist_cort
 
-    if ARGS.v:
-        import matplotlib.pyplot as plt
-        fig, axes = plt.subplots(2)
-        cax0 = axes[0].imshow(_dist_mat_dtw, cmap=plt.cm.gray)
-        axes[0].set_title("dtw")
-        cax1 = axes[1].imshow(_dist_mat_cort, cmap=plt.cm.gray)
-        axes[1].set_title("cort")
-        plt.show()
+    return _dist_mat_dtw, _dist_mat_cort
+
+
+main = partial(util["main"], compute)
+
+if __name__ == "__main__":
+    main()

pyfiles/dtw_cort_dist/V4_c_dtw_c_cort/setup.py

+from setuptools import setup, Extension
+
+version = "0.1"
+
+module_distance = Extension(
+    name="cdtw",
+    sources=["cdtw.c"],
+    extra_compile_args=["-O3", "-mavx2"],
+    extra_link_args=["-O3"],
+)
+
+setup(
+    name="dtw_cort_dist_mat",
+    version=version,
+    description="data scientist tool for time series",
+    long_description="data scientist tool for time series",
+    classifiers=[],  # Get strings from http://pypi.python.org/pypi?%3Aaction=list_classifiers
+    author="IKATS project",
+    author_email="Franck.Thollard@imag.fr",
+    url="http://ikats.imag.fr",
+    license="GPL",
+    include_package_data=True,
+    zip_safe=False,
+    install_requires=["cffi", "numpy", "setuptools"],
+    entry_points="",
+    ext_modules=[module_distance],
+)

todo.md

@@ -16,10 +16,6 @@
 
 - [Pierre, Franck] 31_accelerators.ipynb (que retenir des exemples?)
 
-- [Franck] pyfiles/dtw_cort_dist:
-
-  - option -t for testing (no plot) !!!
-
 - [Pierre] Pythran without Transonic, (loop and vec) for Transonic and Numba
 
 - [Pierre, Franck] parallel (presentation + examples)
@@ -33,4 +29,3 @@
 
 - [Franck] Cleanup
   - old_snippet and delete this directory
-  - pyfiles/dtw_cort_dist (fix or delete broken repositories)