Port to Python 3 & NumPy API 1.7

.gitignore

@@ -3,3 +3,4 @@ build
 *.pyc
 *.sw*
 *.mex*
+.vscode

README.md

@@ -22,9 +22,9 @@ or by running the setup file manually
     cd pyemd
     python setup.py install
 
-Note the code requires Python 2 (Python 3 is not supported) and depends on the
-`numpy` and `scipy` packages. So have those installed first. The build will
-likely fail if it can't find them. For more information, see:
+Note the code requires Python 3 and depends on the `numpy` and `scipy` packages.
+So have those installed first. The build will likely fail if it can't find them.
+For more information, see:
 
  + [NumPy](http://www.numpy.org/): Library for efficient matrix math in Python
  + [SciPy](http://www.scipy.org/): Library for more MATLAB-like functionality

c_emd/emd.c

@@ -11,8 +11,8 @@
 #include <emd.h>
 
 struct basic_variable **initialize_flow(int n_x, double *weight_x,
-                                       int n_y, double *weight_y,
-                                       double **cost);
+                                        int n_y, double *weight_y,
+                                        double **cost);
 struct basic_variable *init_basic(int row, int col, double flow);
 void insert_basic(struct basic_variable **basis, int size,
                   struct basic_variable *node);
@@ -205,14 +205,14 @@ double emd(int n_x, double *weight_x,
  * using the "Northwest Corner Rule"
  */
 struct basic_variable **initialize_flow(int n_x, double *weight_x,
-                                       int n_y, double *weight_y,
-                                       double **cost){
+                                        int n_y, double *weight_y,
+                                        double **cost){
     struct basic_variable **basis;
     struct basic_variable *basic;
     double *remaining_x;
     double *remaining_y;
     int fx, fy, b, B;
-    
+
     b = 0;
     B = n_x + n_y - 1;
 

c_emd/emd.h

@@ -2,7 +2,7 @@
 /* EMD - compute earth mover's distance */
 /****************************************/
 
-#define EPSILON 1e-12
+#define EPSILON 1e-10
 #define FALSE 0
 #define TRUE 1
 

c_emd/pyemd.c

@@ -2,6 +2,7 @@
 /* Python Wrapper for EMD */
 /**************************/
 #include <Python.h>
+#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #include <numpy/arrayobject.h>
 
 #include <emd.h>
@@ -24,15 +25,15 @@ static PyObject *_emd(PyObject *self, PyObject *args) {
     if (arg1 == NULL || arg2 == NULL || arg3 == NULL || return_flows == NULL) {
         return NULL;
     }
-    
+
     // Convert inputs
     // Assumes array data is contiguous
-    n_x = (int) arg1->dimensions[0];
-    n_y = (int) arg2->dimensions[0];
-    weight_x = (double *)arg1->data;
-    weight_y = (double *)arg2->data;
+    n_x = (int) PyArray_DIMS(arg1)[0];
+    n_y = (int) PyArray_DIMS(arg2)[0];
+    weight_x = (double *)PyArray_DATA(arg1);
+    weight_y = (double *)PyArray_DATA(arg2);
     cost = (double **) malloc((size_t) n_x * sizeof(double *));
-    data_ptr = (double *)arg3->data;
+    data_ptr = (double *)PyArray_DATA(arg3);
     for (i = 0; i < n_x; i++) { cost[i] = data_ptr + i*n_y; }
 
     double **flows_data_ptr = NULL;
@@ -44,7 +45,7 @@ static PyObject *_emd(PyObject *self, PyObject *args) {
     }
 
     distance = emd(n_x, weight_x, n_y, weight_y, cost, flows_data_ptr);
-    
+
     free(cost);
     PyObject *retval = NULL;
     if (flows_data_ptr) {
@@ -77,9 +78,17 @@ static PyMethodDef c_emd_methods[] = {
    { NULL, NULL, 0, NULL }
 };
 
-void initc_emd(void)
+static struct PyModuleDef emdmodule = {
+    PyModuleDef_HEAD_INIT,
+    "c_emd",
+    "Earth Mover's Distance",
+    0,
+    c_emd_methods
+};
+
+PyMODINIT_FUNC
+PyInit_c_emd(void)
 {
-    Py_InitModule3("c_emd", c_emd_methods,
-                   "Earth Mover's Distance");
     import_array();
+    return PyModule_Create(&emdmodule);
 }

emd/__init__.py

@@ -1,4 +1,4 @@
 """
 Earth Mover's Distance in Python
 """
-from emd import emd
+from .emd import emd

emd/emd.py

@@ -23,6 +23,31 @@ def emd(X, Y, X_weights=None, Y_weights=None, distance='euclidean',
     @param return_flows : whether to return the flows between cells in addition
         to the distance (default False)
     """
+
+    if X_weights is not None:
+        if not isinstance(X_weights, np.ndarray):
+            raise ValueError("X_weights has to be a Numpy ndarray, not '{}'.".format(type(X_weights)))
+        if not X_weights.flags['C_CONTIGUOUS']:
+            raise ValueError("X_weights has to be a C-contiguous array.")
+        if not X_weights.dtype == np.float64:
+            raise ValueError("X_weights has to be of dtype double (float64).")
+        if np.any(X_weights < 0):
+            raise ValueError("X_weights has to be non-negative.")
+        if not np.isclose(np.sum(X_weights), 1):
+            raise ValueError("X_weights has to sum to 1.")
+
+    if Y_weights is not None:
+        if not isinstance(Y_weights, np.ndarray):
+            raise ValueError("Y_weights has to be a Numpy ndarray, not '{}'.".format(type(Y_weights)))
+        if not Y_weights.flags['C_CONTIGUOUS']:
+            raise ValueError("Y_weights has to be a C-contiguous array.")
+        if not Y_weights.dtype == np.float64:
+            raise ValueError("Y_weights has to be of dtype double (float64).")
+        if np.any(Y_weights < 0):
+            raise ValueError("Y_weights has to be non-negative.")
+        if not np.isclose(np.sum(Y_weights), 1):
+            raise ValueError("Y_weights has to sum to 1.")
+
     if distance != 'precomputed':
         n = len(X)
         m = len(Y)
@@ -39,6 +64,16 @@ def emd(X, Y, X_weights=None, Y_weights=None, distance='euclidean',
     else:
         if D is None:
             raise ValueError("D must be supplied when distance='precomputed'")
+
+        if not isinstance(D, np.ndarray):
+            raise ValueError("D has to be a Numpy ndarray, not '{}'.".format(type(D)))
+        if not D.flags['C_CONTIGUOUS']:
+            raise ValueError("D has to be a C-contiguous array.")
+        if not D.dtype == np.float64:
+            raise ValueError("D has to be of dtype double (float64).")
+        if np.any(D < 0):
+            raise ValueError("D has to be non-negative.")
+
         n, m = D.shape
         if X_weights is None:
             X_weights = np.ones(n)/n

test/test.py

@@ -12,10 +12,10 @@
               [-1.0, 0.0],
               [0.5, 5.0]])
 
-print emd(X, Y)
+print(emd(X, Y))
 
 dist, flows = emd(X, Y, return_flows=True)
-print flows
+print(flows)
 
 dist, flows = emd(X, Z, return_flows=True)
-print flows
+print(flows)