From f8c1c8740f9974dcf4aaf191851d62149dceb91c Mon Sep 17 00:00:00 2001 From: Antoine Rolet Date: Thu, 7 Sep 2017 13:29:46 +0900 Subject: Added MAX_ITER_REACHED flag and warning --- ot/lp/EMD.h | 3 ++- ot/lp/EMD_wrapper.cpp | 21 ++++++++----------- ot/lp/emd_wrap.pyx | 29 +++++++++++++++----------- ot/lp/network_simplex_simple.h | 46 +++++++++++++++++++++++------------------- 4 files changed, 52 insertions(+), 47 deletions(-) (limited to 'ot/lp') 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; i0) { 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, a.data, b.data, M.data, G.data, alpha.data, beta.data, &cost, max_iter) + cdef int resultSolver = EMD_wrap(n1,n2, a.data, b.data, M.data, G.data, alpha.data, beta.data, &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, a.data, b.data, M.data, G.data, alpha.data, beta.data, &cost, max_iter) + cdef int resultSolver = EMD_wrap(n1,n2, a.data, b.data, M.data, G.data, alpha.data, beta.data, &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 -- cgit v1.2.3