summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--ot/lp/EMD.h3
-rw-r--r--ot/lp/EMD_wrapper.cpp21
-rw-r--r--ot/lp/emd_wrap.pyx29
-rw-r--r--ot/lp/network_simplex_simple.h46
-rw-r--r--test/test_ot.py52
5 files changed, 102 insertions, 49 deletions
diff --git a/ot/lp/EMD.h b/ot/lp/EMD.h
index 15e9115..bb486de 100644
--- a/ot/lp/EMD.h
+++ b/ot/lp/EMD.h
@@ -26,7 +26,8 @@ typedef unsigned int node_id_type;
enum ProblemType {
INFEASIBLE,
OPTIMAL,
- UNBOUNDED
+ UNBOUNDED,
+ 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 max_iter);
diff --git a/ot/lp/EMD_wrapper.cpp b/ot/lp/EMD_wrapper.cpp
index 8e74462..92663dc 100644
--- a/ot/lp/EMD_wrapper.cpp
+++ b/ot/lp/EMD_wrapper.cpp
@@ -29,14 +29,18 @@ int EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G,
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
@@ -83,16 +87,7 @@ int EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G,
// Solve the problem with the network simplex algorithm
int ret=net.run();
- if (ret!=(int)net.OPTIMAL) {
- if (ret==(int)net.INFEASIBLE) {
- std::cout << "Infeasible problem";
- }
- if (ret==(int)net.UNBOUNDED)
- {
- std::cout << "Unbounded problem";
- }
- } else
- {
+ if (ret==(int)net.OPTIMAL || ret==(int)net.MAX_ITER_REACHED) {
*cost = 0;
Arc a; di.first(a);
for (; a != INVALID; di.next(a)) {
@@ -105,7 +100,7 @@ int EMD_wrap(int n1, int n2, double *X, double *Y, double *D, double *G,
*(beta + indJ[j-n]) = net.potential(j);
}
- };
+ }
return ret;
diff --git a/ot/lp/emd_wrap.pyx b/ot/lp/emd_wrap.pyx
index 7056e0e..9bea154 100644
--- a/ot/lp/emd_wrap.pyx
+++ b/ot/lp/emd_wrap.pyx
@@ -7,6 +7,7 @@ Cython linker with C solver
#
# License: MIT License
+import warnings
import numpy as np
cimport numpy as np
@@ -15,14 +16,14 @@ cimport cython
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 max_iter)
- cdef enum ProblemType: INFEASIBLE, OPTIMAL, UNBOUNDED
+ int EMD_wrap(int n1,int n2, double *X, double *Y,double *D, double *G, double* alpha, double* beta, double *cost, int numItermax)
+ cdef enum ProblemType: INFEASIBLE, OPTIMAL, UNBOUNDED, MAX_ITER_REACHED
@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):
+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 numItermax):
"""
Solves the Earth Movers distance problem and returns the optimal transport matrix
@@ -49,7 +50,7 @@ def emd_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mod
target histogram
M : (ns,nt) ndarray, float64
loss matrix
- max_iter : int
+ numItermax : int
The maximum number of iterations before stopping the optimization
algorithm if it has not converged.
@@ -76,18 +77,20 @@ def emd_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mod
b=np.ones((n2,))/n2
# calling the function
- cdef int resultSolver = 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)
+ cdef int resultSolver = EMD_wrap(n1,n2,<double*> a.data,<double*> b.data,<double*> M.data,<double*> G.data, <double*> alpha.data, <double*> beta.data, <double*> &cost, numItermax)
if resultSolver != OPTIMAL:
if resultSolver == INFEASIBLE:
- print("Problem infeasible. Try to increase numItermax.")
+ warnings.warn("Problem infeasible. Check that a and b are in the simplex")
elif resultSolver == UNBOUNDED:
- print("Problem unbounded")
+ warnings.warn("Problem unbounded")
+ elif resultSolver == MAX_ITER_REACHED:
+ warnings.warn("numItermax reached before optimality. Try to increase numItermax.")
return G, alpha, beta
@cython.boundscheck(False)
@cython.wraparound(False)
-def emd2_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):
+def emd2_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 numItermax):
"""
Solves the Earth Movers distance problem and returns the optimal transport loss
@@ -114,7 +117,7 @@ def emd2_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mo
target histogram
M : (ns,nt) ndarray, float64
loss matrix
- max_iter : int
+ numItermax : int
The maximum number of iterations before stopping the optimization
algorithm if it has not converged.
@@ -140,12 +143,14 @@ def emd2_c( np.ndarray[double, ndim=1, mode="c"] a,np.ndarray[double, ndim=1, mo
if not len(b):
b=np.ones((n2,))/n2
# calling the function
- cdef int resultSolver = 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)
+ cdef int resultSolver = EMD_wrap(n1,n2,<double*> a.data,<double*> b.data,<double*> M.data,<double*> G.data, <double*> alpha.data, <double*> beta.data, <double*> &cost, numItermax)
if resultSolver != OPTIMAL:
if resultSolver == INFEASIBLE:
- print("Problem infeasible. Try to inscrease numItermax.")
+ warnings.warn("Problem infeasible. Check that a and b are in the simplex")
elif resultSolver == UNBOUNDED:
- print("Problem unbounded")
+ warnings.warn("Problem unbounded")
+ elif resultSolver == MAX_ITER_REACHED:
+ warnings.warn("numItermax reached before optimality. Try to increase numItermax.")
return cost, alpha, beta
diff --git a/ot/lp/network_simplex_simple.h b/ot/lp/network_simplex_simple.h
index a7743ee..7c6a4ce 100644
--- a/ot/lp/network_simplex_simple.h
+++ b/ot/lp/network_simplex_simple.h
@@ -34,7 +34,8 @@
#endif
-#define EPSILON 10*2.2204460492503131e-016
+#define EPSILON 2.2204460492503131e-15
+#define _EPSILON 1e-8
#define MAX_DEBUG_ITER 100000
@@ -260,7 +261,9 @@ namespace lemon {
/// The objective function of the problem is unbounded, i.e.
/// there is a directed cycle having negative total cost and
/// infinite upper bound.
- UNBOUNDED
+ UNBOUNDED,
+ /// The maximum number of iteration has been reached
+ MAX_ITER_REACHED
};
/// \brief Constants for selecting the type of the supply constraints.
@@ -683,7 +686,7 @@ namespace lemon {
/// \see resetParams(), reset()
ProblemType run() {
#if DEBUG_LVL>0
- std::cout << "OPTIMAL = " << OPTIMAL << "\nINFEASIBLE = " << INFEASIBLE << "nUNBOUNDED = " << UNBOUNDED << "\n";
+ std::cout << "OPTIMAL = " << OPTIMAL << "\nINFEASIBLE = " << INFEASIBLE << "\nUNBOUNDED = " << UNBOUNDED << "\nMAX_ITER_REACHED" << MAX_ITER_REACHED\n";
#endif
if (!init()) return INFEASIBLE;
@@ -941,15 +944,15 @@ namespace lemon {
// Initialize internal data structures
bool init() {
if (_node_num == 0) return false;
- /*
+
// Check the sum of supply values
_sum_supply = 0;
for (int i = 0; i != _node_num; ++i) {
_sum_supply += _supply[i];
}
- if ( !((_stype == GEQ && _sum_supply <= _epsilon ) ||
- (_stype == LEQ && _sum_supply >= -_epsilon )) ) return false;
- */
+ if ( fabs(_sum_supply) > _EPSILON ) return false;
+
+ _sum_supply = 0;
// Initialize artifical cost
Cost ART_COST;
@@ -1416,13 +1419,11 @@ namespace lemon {
ProblemType start() {
PivotRuleImpl pivot(*this);
double prevCost=-1;
+ ProblemType retVal = OPTIMAL;
// Perform heuristic initial pivots
if (!initialPivots()) return UNBOUNDED;
-#if DEBUG_LVL>0
- int niter=0;
-#endif
int iter_number=0;
//pivot.setDantzig(true);
// Execute the Network Simplex algorithm
@@ -1431,12 +1432,13 @@ namespace lemon {
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;
+ retVal = MAX_ITER_REACHED;
break;
}
#if DEBUG_LVL>0
- if(niter>MAX_DEBUG_ITER)
+ if(iter_number>MAX_DEBUG_ITER)
break;
- if(++niter%1000==0||niter%1000==1){
+ if(iter_number%1000==0||iter_number%1000==1){
double curCost=totalCost();
double sumFlow=0;
double a;
@@ -1445,7 +1447,7 @@ namespace lemon {
for (int i=0; i<_flow.size(); i++) {
sumFlow+=_state[i]*_flow[i];
}
- std::cout << "Sum of the flow " << std::setprecision(20) << sumFlow << "\n" << niter << " iterations, current cost=" << curCost << "\nReduced cost=" << _state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -_pi[_target[in_arc]]) << "\nPrecision = "<< -EPSILON*(a) << "\n";
+ std::cout << "Sum of the flow " << std::setprecision(20) << sumFlow << "\n" << iter_number << " iterations, current cost=" << curCost << "\nReduced cost=" << _state[in_arc] * (_cost[in_arc] + _pi[_source[in_arc]] -_pi[_target[in_arc]]) << "\nPrecision = "<< -EPSILON*(a) << "\n";
std::cout << "Arc in = (" << _node_id(_source[in_arc]) << ", " << _node_id(_target[in_arc]) <<")\n";
std::cout << "Supplies = (" << _supply[_source[in_arc]] << ", " << _supply[_target[in_arc]] << ")\n";
std::cout << _cost[in_arc] << "\n";
@@ -1503,15 +1505,17 @@ namespace lemon {
std::cout << "Sum of the flow " << sumFlow << "\n"<< niter <<" iterations, current cost=" << totalCost() << "\n";
#endif
// Check feasibility
- for (int e = _search_arc_num; e != _all_arc_num; ++e) {
- if (_flow[e] != 0){
- if (abs(_flow[e]) > EPSILON)
- return INFEASIBLE;
- else
- _flow[e]=0;
+ if( retVal == OPTIMAL){
+ for (int e = _search_arc_num; e != _all_arc_num; ++e) {
+ if (_flow[e] != 0){
+ if (abs(_flow[e]) > EPSILON)
+ return INFEASIBLE;
+ else
+ _flow[e]=0;
+ }
}
- }
+ }
// Shift potentials to meet the requirements of the GEQ/LEQ type
// optimality conditions
@@ -1537,7 +1541,7 @@ namespace lemon {
}
}
- return OPTIMAL;
+ return retVal;
}
}; //class NetworkSimplexSimple
diff --git a/test/test_ot.py b/test/test_ot.py
index 6f0f7c9..8a19cf6 100644
--- a/test/test_ot.py
+++ b/test/test_ot.py
@@ -8,6 +8,7 @@ import numpy as np
import ot
from ot.datasets import get_1D_gauss as gauss
+import warnings
def test_doctest():
@@ -100,9 +101,56 @@ def test_emd2_multi():
np.testing.assert_allclose(emd1, emdn)
-def test_dual_variables():
- # %% parameters
+def test_warnings():
+ n = 100 # nb bins
+ m = 100 # nb bins
+
+ mean1 = 30
+ mean2 = 50
+
+ # bin positions
+ x = np.arange(n, dtype=np.float64)
+ y = np.arange(m, dtype=np.float64)
+
+ # Gaussian distributions
+ a = gauss(n, m=mean1, s=5) # m= mean, s= std
+
+ b = gauss(m, m=mean2, s=10)
+ # loss matrix
+ M = ot.dist(x.reshape((-1, 1)), y.reshape((-1, 1))) ** (1. / 2)
+ # M/=M.max()
+
+ # %%
+
+ print('Computing {} EMD '.format(1))
+ G, alpha, beta = ot.emd(a, b, M, dual_variables=True)
+ with warnings.catch_warnings(record=True) as w:
+ # Cause all warnings to always be triggered.
+ warnings.simplefilter("always")
+ # Trigger a warning.
+ print('Computing {} EMD '.format(1))
+ G, alpha, beta = ot.emd(a, b, M, dual_variables=True, numItermax=1)
+ # Verify some things
+ assert "numItermax" in str(w[-1].message)
+ assert len(w) == 1
+ # Trigger a warning.
+ a[0]=100
+ print('Computing {} EMD '.format(2))
+ G, alpha, beta = ot.emd(a, b, M, dual_variables=True)
+ # Verify some things
+ assert "infeasible" in str(w[-1].message)
+ assert len(w) == 2
+ # Trigger a warning.
+ a[0]=-1
+ print('Computing {} EMD '.format(2))
+ G, alpha, beta = ot.emd(a, b, M, dual_variables=True)
+ # Verify some things
+ assert "infeasible" in str(w[-1].message)
+ assert len(w) == 3
+
+
+def test_dual_variables():
n = 5000 # nb bins
m = 6000 # nb bins
generated by cgit v1.2.3 (git 2.39.1) at 2024-06-03 06:36:57 +0000