From c123974f420cebc3c39f103988d594df67903093 Mon Sep 17 00:00:00 2001 From: Alex Leontiev Date: Sat, 20 Jul 2013 15:14:02 +0300 Subject: [PATCH] Eliminated all the calls to std::find() This is done by keeping indexToRow vector, that keeps the information, opposite to those kept by N and B. That is, while N and B help to determine which variable corresponds to given column in column-vector c or row in matrix b, indexToRow helps to determine the corresponding row/column for a given variable. At this point, I'm waiting for comments from pull request reviewer and not working on any upgrades. Comments are appreciated, as usual. --- modules/optim/include/opencv2/optim.hpp | 1 - modules/optim/src/lpsolver.cpp | 68 +++++++++++++++++++-------------- modules/optim/test/test_lpsolver.cpp | 13 ------- 3 files changed, 40 insertions(+), 42 deletions(-) diff --git a/modules/optim/include/opencv2/optim.hpp b/modules/optim/include/opencv2/optim.hpp index 13ed865..0bbedad 100644 --- a/modules/optim/include/opencv2/optim.hpp +++ b/modules/optim/include/opencv2/optim.hpp @@ -47,7 +47,6 @@ namespace cv{namespace optim { - using namespace std; //!the return codes for solveLP() function enum { diff --git a/modules/optim/src/lpsolver.cpp b/modules/optim/src/lpsolver.cpp index e746c39..0690387 100644 --- a/modules/optim/src/lpsolver.cpp +++ b/modules/optim/src/lpsolver.cpp @@ -3,11 +3,8 @@ #include #include -#define ALEX_DEBUG - namespace cv{namespace optim{ using std::vector; -using namespace std; #ifdef ALEX_DEBUG #define dprintf(x) printf x @@ -46,12 +43,14 @@ static void print_simplex_state(const Mat& c,const Mat& b,double v,const std::ve it also initializes N and B and does not make any assumptions about their init values * @return SOLVELP_UNFEASIBLE if problem is unfeasible, 0 if feasible. */ -static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B); -static inline void pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B, int leaving_index,int entering_index); +static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B,vector& indexToRow); +static inline void pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B,int leaving_index, + int entering_index,vector& indexToRow); /**@return SOLVELP_UNBOUNDED means the problem is unbdd, SOLVELP_MULTI means multiple solutions, SOLVELP_SINGLE means one solution. */ -static int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B); +static int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B,vector& indexToRow); static void swap_columns(Mat_& A,int col1,int col2); +#define SWAP(type,a,b) {type tmp=(a);(a)=(b);(b)=tmp;} //return codes:-2 (no_sol - unbdd),-1(no_sol - unfsbl), 0(single_sol), 1(multiple_sol=>least_l2_norm) int solveLP(const Mat& Func, const Mat& Constr, Mat& z){ @@ -75,15 +74,16 @@ int solveLP(const Mat& Func, const Mat& Constr, Mat& z){ Constr.convertTo(bigB.colRange(1,bigB.cols),CV_64FC1); double v=0; vector N,B; + vector indexToRow; - if(initialize_simplex(bigC,bigB,v,N,B)==SOLVELP_UNFEASIBLE){ + if(initialize_simplex(bigC,bigB,v,N,B,indexToRow)==SOLVELP_UNFEASIBLE){ return SOLVELP_UNFEASIBLE; } Mat_ c=bigC.colRange(1,bigC.cols), b=bigB.colRange(1,bigB.cols); int res=0; - if((res=inner_simplex(c,b,v,N,B))==SOLVELP_UNBOUNDED){ + if((res=inner_simplex(c,b,v,N,B,indexToRow))==SOLVELP_UNBOUNDED){ return SOLVELP_UNBOUNDED; } @@ -92,17 +92,17 @@ int solveLP(const Mat& Func, const Mat& Constr, Mat& z){ MatIterator_ it=z.begin(); for(int i=1;i<=c.cols;i++,it++){ std::vector::iterator pos=B.begin(); - if((pos=std::find(B.begin(),B.end(),i))==B.end()){ + if(indexToRow[i](pos-B.begin(),b.cols-1); + *it=b.at(indexToRow[i]-N.size(),b.cols-1); } } return res; } -static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B){ +static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B,vector& indexToRow){ N.resize(c.cols); N[0]=0; for (std::vector::iterator it = N.begin()+1 ; it != N.end(); ++it){ @@ -113,6 +113,11 @@ static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector< for (std::vector::iterator it = B.begin()+1 ; it != B.end(); ++it){ *it=it[-1]+1; } + indexToRow.resize(c.cols+b.rows); + indexToRow[0]=0; + for (std::vector::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){ + *it=it[-1]+1; + } v=0; int k=0; @@ -128,6 +133,9 @@ static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector< if(b(k,b.cols-1)>=0){ N.erase(N.begin()); + for (std::vector::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){ + --(*it); + } return 0; } @@ -141,28 +149,29 @@ static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector< print_simplex_state(c,b,v,N,B); dprintf(("\tWE MAKE PIVOT\n")); - pivot(c,b,v,N,B,k,0); + pivot(c,b,v,N,B,k,0,indexToRow); print_simplex_state(c,b,v,N,B); - inner_simplex(c,b,v,N,B); + inner_simplex(c,b,v,N,B,indexToRow); dprintf(("\tAFTER INNER_SIMPLEX\n")); print_simplex_state(c,b,v,N,B); - vector::iterator iterator=std::find(B.begin(),B.end(),0); - if(iterator!=B.end()){ - int iterator_offset=iterator-B.begin(); + if(indexToRow[0]>=N.size()){ + int iterator_offset=indexToRow[0]-N.size(); if(b(iterator_offset,b.cols-1)>0){ return SOLVELP_UNFEASIBLE; } - pivot(c,b,v,N,B,iterator_offset,0); + pivot(c,b,v,N,B,iterator_offset,0,indexToRow); } + vector::iterator iterator; { - iterator=std::find(N.begin(),N.end(),0); - int iterator_offset=iterator-N.begin(); + int iterator_offset=indexToRow[0]; + iterator=N.begin()+iterator_offset; std::iter_swap(iterator,N.begin()); + SWAP(int,indexToRow[*iterator],indexToRow[0]); swap_columns(c,iterator_offset,0); swap_columns(b,iterator_offset,0); } @@ -174,14 +183,14 @@ static int initialize_simplex(Mat_& c, Mat_& b,double& v,vector< c=0; v=0; for(int I=1;I& c, Mat_& b,double& v,vector< print_simplex_state(c,b,v,N,B); N.erase(N.begin()); + for (std::vector::iterator it = indexToRow.begin()+1 ; it != indexToRow.end(); ++it){ + --(*it); + } return 0; } -static int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B){ +static int inner_simplex(Mat_& c, Mat_& b,double& v,vector& N,vector& B,vector& indexToRow){ int count=0; while(1){ dprintf(("iteration #%d\n",count)); @@ -248,7 +260,7 @@ static int inner_simplex(Mat_& c, Mat_& b,double& v,vector& } dprintf(("the tightest constraint is in row %d with %g\n",l,min)); - pivot(c,b,v,N,B,l,e); + pivot(c,b,v,N,B,l,e,indexToRow); dprintf(("objective, v=%g\n",v)); print_matrix(c); @@ -261,7 +273,8 @@ static int inner_simplex(Mat_& c, Mat_& b,double& v,vector& } } -static inline void pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B, int leaving_index,int entering_index){ +static inline void pivot(Mat_& c,Mat_& b,double& v,vector& N,vector& B, + int leaving_index,int entering_index,vector& indexToRow){ double Coef=b(leaving_index,entering_index); for(int i=0;i& c,Mat_& b,double& v,vector& dprintf(("v was %g\n",v)); v+=Coef*b(leaving_index,b.cols-1); - int tmp=N[entering_index]; - N[entering_index]=B[leaving_index]; - B[leaving_index]=tmp; + SWAP(int,N[entering_index],B[leaving_index]); + SWAP(int,indexToRow[N[entering_index]],indexToRow[B[leaving_index]]); } static inline void swap_columns(Mat_& A,int col1,int col2){ diff --git a/modules/optim/test/test_lpsolver.cpp b/modules/optim/test/test_lpsolver.cpp index 68bc693..5a7f4c4 100644 --- a/modules/optim/test/test_lpsolver.cpp +++ b/modules/optim/test/test_lpsolver.cpp @@ -81,19 +81,6 @@ TEST(Optim_LpSolver, regression_multiple_solutions){ ASSERT_EQ(res,1); ASSERT_EQ(z.dot(A),1); } - - if(false){ - //cormen's example from chapter about initialize_simplex - //online solver told it has inf many solutions, but I'm not sure - A=(cv::Mat_(2,1)<<2,-1); - B=(cv::Mat_(2,3)<<2,-1,2,1,-5,-4); - std::cout<<"here A goes\n"<