Silence std::cerr calls inside of python

add use_cerr flag to backend

add use_cerr to lp.__init__ and passthrough to c++

add use_cerr defs to pyx

fix call to use_cerr to update from ot.backend

def use_cerr in emd.h

emd wrapper cerr

cpp _use_cerr

releases.md

Author: stanleyjs

RELEASES.md

@@ -10,6 +10,7 @@
 
 - Fixed an issue where we could not ask TorchBackend to place a random tensor on GPU
   (Issue #371, PR #373)
+- Fixed an issue where hitting iteration limits would be reported to stderr by std:cerr regardless of Python's stderr stream status.
 
 ## 0.8.2
 

ot/backend.py

@@ -127,6 +127,9 @@
     tf = False
     tf_type = float
 
+def get_cerr():
+    return use_cerr
+use_cerr = False
 
 str_type_error = "All array should be from the same type/backend. Current types are : {}"
 

ot/lp/EMD.h

@@ -28,8 +28,8 @@ enum ProblemType {
     MAX_ITER_REACHED
 };
 
-int EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int maxIter);
-int EMD_wrap_omp(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int maxIter, int numThreads);
+int EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int maxIter, bool use_cerr);
+int EMD_wrap_omp(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int maxIter, int numThreads, bool use_cerr);
 
 
 

ot/lp/EMD_wrapper.cpp

@@ -20,7 +20,8 @@
 
 
 int EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G,
-                double* alpha, double* beta, double *cost, int maxIter)  {
+                double* alpha, double* beta, double *cost, int maxIter, 
+                bool use_cerr=true)  {
     // beware M and C are stored in row major C style!!!
 
     using namespace lemon;
@@ -54,7 +55,7 @@ int EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G,
     std::vector<int> indI(n), indJ(m);
     std::vector<double> weights1(n), weights2(m);
     Digraph di(n, m);
-    NetworkSimplexSimple<Digraph,double,double, node_id_type> net(di, true, n+m, ((int64_t)n)*((int64_t)m), maxIter);
+    NetworkSimplexSimple<Digraph,double,double, node_id_type> net(di, true, n+m, ((int64_t)n)*((int64_t)m), maxIter, use_cerr);
 
     // Set supply and demand, don't account for 0 values (faster)
 
@@ -122,7 +123,8 @@ 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, int maxIter, int numThreads)  {
+             double* alpha, double* beta, double *cost, int maxIter, int numThreads,
+             bool use_cerr=true)  {
     // beware M and C are stored in row major C style!!!
 
     using namespace lemon_omp;
@@ -156,7 +158,7 @@ int EMD_wrap_omp(int n1, int n2, double *X, double *Y, double *D, double *G,
     std::vector<int> indI(n), indJ(m);
     std::vector<double> weights1(n), weights2(m);
     Digraph di(n, m);
-    NetworkSimplexSimple<Digraph,double,double, node_id_type> net(di, true, n+m, ((int64_t)n)*((int64_t)m), maxIter, numThreads);
+    NetworkSimplexSimple<Digraph,double,double, node_id_type> net(di, true, n+m, ((int64_t)n)*((int64_t)m), maxIter, numThreads, use_cerr);
 
     // Set supply and demand, don't account for 0 values (faster)
 

ot/lp/__init__.py

@@ -24,7 +24,7 @@
 
 from ..utils import dist, list_to_array
 from ..utils import parmap
-from ..backend import get_backend
+from ..backend import get_backend, get_cerr
 
 __all__ = ['emd', 'emd2', 'barycenter', 'free_support_barycenter', 'cvx', ' emd_1d_sorted',
            'emd_1d', 'emd2_1d', 'wasserstein_1d', 'generalized_free_support_barycenter']
@@ -330,7 +330,8 @@ def emd(a, b, M, numItermax=100000, log=False, center_dual=True, numThreads=1):
 
     numThreads = check_number_threads(numThreads)
 
-    G, cost, u, v, result_code = emd_c(a, b, M, numItermax, numThreads)
+    G, cost, u, v, result_code = emd_c(a, b, M, numItermax, numThreads,
+                                        use_cerr=get_cerr())
 
     if center_dual:
         u, v = center_ot_dual(u, v, a, b)
@@ -489,7 +490,8 @@ def emd2(a, b, M, processes=1,
         def f(b):
             bsel = b != 0
 
-            G, cost, u, v, result_code = emd_c(a, b, M, numItermax, numThreads)
+            G, cost, u, v, result_code = emd_c(a, b, M, numItermax, numThreads,
+                                        use_cerr=get_cerr())
 
             if center_dual:
                 u, v = center_ot_dual(u, v, a, b)
@@ -521,7 +523,8 @@ def f(b):
     else:
         def f(b):
             bsel = b != 0
-            G, cost, u, v, result_code = emd_c(a, b, M, numItermax, numThreads)
+            G, cost, u, v, result_code = emd_c(a, b, M, numItermax, numThreads,
+                                        use_cerr=get_cerr())
 
             if center_dual:
                 u, v = center_ot_dual(u, v, a, b)

ot/lp/emd_wrap.pyx

@@ -9,7 +9,7 @@ Cython linker with C solver
 
 import numpy as np
 cimport numpy as np
-
+from libcpp cimport bool as bool_t
 from ..utils import dist
 
 cimport cython
@@ -19,8 +19,8 @@ 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, int maxIter) nogil
-    int EMD_wrap_omp(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int maxIter, int numThreads) nogil
+    int EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int maxIter, bool_t use_cerr) nogil
+    int EMD_wrap_omp(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int maxIter, int numThreads, bool_t use_cerr) nogil
     cdef enum ProblemType: INFEASIBLE, OPTIMAL, UNBOUNDED, MAX_ITER_REACHED
 
 
@@ -39,7 +39,7 @@ def check_result(result_code):
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def emd_c(np.ndarray[double, ndim=1, mode="c"] a, np.ndarray[double, ndim=1, mode="c"]  b, np.ndarray[double, ndim=2, mode="c"]  M, int max_iter, int numThreads):
+def emd_c(np.ndarray[double, ndim=1, mode="c"] a, np.ndarray[double, ndim=1, mode="c"]  b, np.ndarray[double, ndim=2, mode="c"]  M, int max_iter, int numThreads, bool_t use_cerr=False):
     """
         Solves the Earth Movers distance problem and returns the optimal transport matrix
 
@@ -111,9 +111,9 @@ 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)
+            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, use_cerr)
         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_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, use_cerr)
     return G, cost, alpha, beta, result_code
 
 

ot/lp/network_simplex_simple.h

@@ -233,12 +233,13 @@ namespace lemon {
         /// mixed order in the internal data structure.
         /// In special cases, it could lead to better overall performance,
         /// but it is usually slower. Therefore it is disabled by default.
-        NetworkSimplexSimple(const GR& graph, bool arc_mixing, int nbnodes, ArcsType nb_arcs, size_t maxiters) :
+        NetworkSimplexSimple(const GR& graph, bool arc_mixing, int nbnodes, ArcsType nb_arcs, size_t maxiters, bool use_cerr=true) :
         _graph(graph),  //_arc_id(graph),
         _arc_mixing(arc_mixing), _init_nb_nodes(nbnodes), _init_nb_arcs(nb_arcs),
         MAX(std::numeric_limits<Value>::max()),
         INF(std::numeric_limits<Value>::has_infinity ?
-            std::numeric_limits<Value>::infinity() : MAX)
+            std::numeric_limits<Value>::infinity() : MAX),
+        _use_cerr(use_cerr)
         {
             // Reset data structures
             reset();
@@ -292,6 +293,7 @@ namespace lemon {
 
 
     private:
+        bool _use_cerr;
 
         size_t max_iter;
         TEMPLATE_DIGRAPH_TYPEDEFS(GR);
@@ -334,6 +336,7 @@ namespace lemon {
         IntVector _source;  // keep nodes as integers
         IntVector _target;
         bool _arc_mixing;
+
     public:
         // Node and arc data
         CostVector _cost;
@@ -1433,8 +1436,11 @@ namespace lemon {
             while (pivot.findEnteringArc()) {
                 if(max_iter > 0 && ++iter_number>=max_iter&&max_iter>0){
                     char errMess[1000];
-                    sprintf( errMess, "RESULT MIGHT BE INACURATE\nMax number of iteration reached, currently \%d. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher\n",iter_number );
-                    std::cerr << errMess;
+                    
+                    if(_use_cerr){
+                        sprintf( errMess, "RESULT MIGHT BE INACURATE\nMax number of iteration reached, currently \%d. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher\n",iter_number );
+                        std::cerr << errMess;
+                    }
 					retVal = MAX_ITER_REACHED;
                     break;
                 }

ot/lp/network_simplex_simple_omp.h

@@ -244,12 +244,13 @@ namespace lemon_omp {
 		/// mixed order in the internal data structure.
 		/// In special cases, it could lead to better overall performance,
 		/// but it is usually slower. Therefore it is disabled by default.
-		NetworkSimplexSimple(const GR& graph, bool arc_mixing, int nbnodes, ArcsType nb_arcs, size_t maxiters = 0, int numThreads=-1) :
+		NetworkSimplexSimple(const GR& graph, bool arc_mixing, int nbnodes, ArcsType nb_arcs, size_t maxiters = 0, int numThreads=-1, bool use_cerr=true) :
 			_graph(graph),  //_arc_id(graph),
 			_arc_mixing(arc_mixing), _init_nb_nodes(nbnodes), _init_nb_arcs(nb_arcs),
 			MAX(std::numeric_limits<Value>::max()),
 			INF(std::numeric_limits<Value>::has_infinity ?
-				std::numeric_limits<Value>::infinity() : MAX)
+				std::numeric_limits<Value>::infinity() : MAX),
+			_use_cerr(use_cerr)
 		{
 			// Reset data structures
 			reset();
@@ -317,6 +318,7 @@ namespace lemon_omp {
 
 
 	private:
+		bool _use_cerr;
 		size_t max_iter;
 		int num_threads;
 		TEMPLATE_DIGRAPH_TYPEDEFS(GR);
@@ -1611,8 +1613,10 @@ namespace lemon_omp {
 
 				} else {
 					char errMess[1000];
-					sprintf( errMess, "RESULT MIGHT BE INACURATE\nMax number of iteration reached, currently \%d. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher\n",iter_number );
-					std::cerr << errMess;
+					if(_use_cerr){
+						sprintf( errMess, "RESULT MIGHT BE INACURATE\nMax number of iteration reached, currently \%d. Sometimes iterations go on in cycle even though the solution has been reached, to check if it's the case here have a look at the minimal reduced cost. If it is very close to machine precision, you might actually have the correct solution, if not try setting the maximum number of iterations a bit higher\n",iter_number );
+						std::cerr << errMess;
+					}
 					retVal =  MAX_ITER_REACHED;
 					break;
 				}