removed openmp implementation for the simplex algorithm since no real runtime gains observed

Author: nathanneike

ot/lp/EMD.h

@@ -30,8 +30,6 @@ enum ProblemType {
 };
 
 int EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, uint64_t maxIter, double* alpha_init, double* beta_init);
-int EMD_wrap_omp(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, uint64_t maxIter, int numThreads);
-
 int EMD_wrap_sparse(
     int n1,                      
     int n2,                       

ot/lp/EMD_wrapper.cpp

@@ -14,7 +14,6 @@
 
 
 #include "network_simplex_simple.h"
-#include "network_simplex_simple_omp.h"
 #include "sparse_bipartitegraph.h"
 #include "EMD.h"
 #include <cstdint>
@@ -300,87 +299,6 @@ int EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G,
 
 
 
-int EMD_wrap_omp(int n1, int n2, double *X, double *Y, double *D, double *G,
-             double* alpha, double* beta, double *cost, uint64_t maxIter, int numThreads)  {
-    // beware M and C are stored in row major C style!!!
-
-    using namespace lemon_omp;
-    uint64_t n, m, cur;
-
-    typedef FullBipartiteDigraph Digraph;
-    DIGRAPH_TYPEDEFS(Digraph);
-
-    // Get the number of non zero coordinates for r and c
-    n=0;
-    for (int i=0; i<n1; i++) {
-        double val=*(X+i);
-        if (val>0) {
-            n++;
-        }else if(val<0){
-            return INFEASIBLE;
-        }
-    }
-    m=0;
-    for (int i=0; i<n2; i++) {
-        double val=*(Y+i);
-        if (val>0) {
-            m++;
-        }else if(val<0){
-            return INFEASIBLE;
-        }
-    }
-
-    // Define the graph
-
-    std::vector<uint64_t> indI(n), indJ(m);
-    std::vector<double> weights1(n), weights2(m);
-    Digraph di(n, m);
-    const SetupPolicy policy = make_setup_policy(n, m, n1, n2, false);
-    NetworkSimplexSimple<Digraph,double,double, node_id_type> net(
-        di, policy.use_arc_mixing, (int) (n + m), n * m, maxIter, numThreads
-    );
-
-    // Set supply and demand, don't account for 0 values (faster)
-
-    cur=0;
-    for (uint64_t i=0; i<n1; i++) {
-        double val=*(X+i);
-        if (val>0) {
-            weights1[ cur ] = val;
-            indI[cur++]=i;
-        }
-    }
-
-    // Demand is actually negative supply...
-
-    cur=0;
-    for (uint64_t i=0; i<n2; i++) {
-        double val=*(Y+i);
-        if (val>0) {
-            weights2[ cur ] = -val;
-            indJ[cur++]=i;
-        }
-    }
-
-
-    net.supplyMap(&weights1[0], (int) n, &weights2[0], (int) m);
-
-    setup_explicit_arc_costs(net, di, D, n2, indI, indJ, n, m);
-
-    // Solve the problem with the network simplex algorithm
-
-    int ret=net.run();
-    if (ret==(int)net.OPTIMAL || ret==(int)net.MAX_ITER_REACHED) {
-        *cost = 0;
-        extract_compressed_support(
-            net, di, INVALID, D, G, alpha, beta, cost, indI, indJ, n, n2
-        );
-
-    }
-
-    return ret;
-}
-
 // ============================================================================
 // SPARSE VERSION: Accepts edge list instead of dense cost matrix
 // ============================================================================

ot/lp/_barycenter_solvers.py

@@ -213,9 +213,9 @@ def free_support_barycenter(
         Print information along iterations
     log : bool, optional
         record log if True
-    numThreads: int or "max", optional (default=1, i.e. OpenMP is not used)
-        If compiled with OpenMP, chooses the number of threads to parallelize.
-        "max" selects the highest number possible.
+    numThreads: int or "max", optional (default=1)
+        Deprecated compatibility parameter forwarded to EMD. The network
+        simplex solver no longer uses OpenMP, so this parameter is ignored.
 
     Returns
     -------
@@ -349,9 +349,9 @@ def generalized_free_support_barycenter(
         Print information along iterations
     log : bool, optional
         record log if True
-    numThreads: int or "max", optional (default=1, i.e. OpenMP is not used)
-        If compiled with OpenMP, chooses the number of threads to parallelize.
-        "max" selects the highest number possible.
+    numThreads: int or "max", optional (default=1)
+        Deprecated compatibility parameter forwarded to EMD. The network
+        simplex solver no longer uses OpenMP, so this parameter is ignored.
     eps: Stability coefficient for the change of variable matrix inversion
         If the :math:`\mathbf{P}_i^T` matrices don't span :math:`\mathbb{R}^d`, the problem is ill-defined and a matrix
         inversion will fail. In this case one may set eps=1e-8 and get a solution anyway (which may make little sense)

ot/lp/_network_simplex.py

@@ -287,9 +287,9 @@ def emd(
     center_dual: boolean, optional (default=True)
         If True, centers the dual potential using function
         :py:func:`ot.lp.center_ot_dual`.
-    numThreads: int or "max", optional (default=1, i.e. OpenMP is not used)
-        If compiled with OpenMP, chooses the number of threads to parallelize.
-        "max" selects the highest number possible.
+    numThreads: int or "max", optional (default=1)
+        Deprecated compatibility parameter. The network simplex solver no
+        longer uses OpenMP, so this parameter is ignored.
     check_marginals: bool, optional (default=True)
         If True, checks that the marginals mass are equal. If False, skips the
         check.
@@ -590,9 +590,9 @@ def emd2(
     center_dual: boolean, optional (default=True)
         If True, centers the dual potential using function
         :py:func:`ot.lp.center_ot_dual`.
-    numThreads: int or "max", optional (default=1, i.e. OpenMP is not used)
-        If compiled with OpenMP, chooses the number of threads to parallelize.
-        "max" selects the highest number possible.
+    numThreads: int or "max", optional (default=1)
+        Deprecated compatibility parameter. The network simplex solver no
+        longer uses OpenMP, so this parameter is ignored.
     check_marginals: bool, optional (default=True)
         If True, checks that the marginals mass are equal. If False, skips the
         check.

ot/lp/emd_wrap.pyx

@@ -21,7 +21,6 @@ import warnings
 
 cdef extern from "EMD.h":
     int EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, uint64_t maxIter, double* alpha_init, double* beta_init) nogil
-    int EMD_wrap_omp(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, uint64_t maxIter, int numThreads) nogil
     int EMD_wrap_sparse(int n1, int n2, double *X, double *Y, uint64_t n_edges, uint64_t *edge_sources, uint64_t *edge_targets, double *edge_costs, uint64_t *flow_sources_out, uint64_t *flow_targets_out, double *flow_values_out, uint64_t *n_flows_out, uint64_t max_flows_out, double *alpha, double *beta, double *cost, uint64_t maxIter, double* alpha_init, double* beta_init) nogil
     int EMD_wrap_lazy(int n1, int n2, double *X, double *Y, double *coords_a, double *coords_b, int dim, int metric, uint64_t *flow_sources_out, uint64_t *flow_targets_out, double *flow_values_out, uint64_t *n_flows_out, uint64_t max_flows_out, double* alpha, double* beta, double *cost, uint64_t maxIter, double* alpha_init, double* beta_init) nogil
     cdef enum ProblemType: INFEASIBLE, OPTIMAL, UNBOUNDED, MAX_ITER_REACHED
@@ -130,10 +129,7 @@ def emd_c(np.ndarray[double, ndim=1, mode="c"] a, np.ndarray[double, ndim=1, mod
 
     # calling the function
     with nogil:
-        if numThreads == 1:
-            result_code = EMD_wrap(n1, n2, <double*> a.data, <double*> b.data, <double*> M.data, <double*> G.data, <double*> alpha.data, <double*> beta.data, <double*> &cost, max_iter, alpha_init_ptr, beta_init_ptr)
-        else:
-            result_code = EMD_wrap_omp(n1, n2, <double*> a.data, <double*> b.data, <double*> M.data, <double*> G.data, <double*> alpha.data, <double*> beta.data, <double*> &cost, max_iter, numThreads)
+        result_code = EMD_wrap(n1, n2, <double*> a.data, <double*> b.data, <double*> M.data, <double*> G.data, <double*> alpha.data, <double*> beta.data, <double*> &cost, max_iter, alpha_init_ptr, beta_init_ptr)
     return G, cost, alpha, beta, result_code